Math::NumSeq::AlgebraicContinued(3pm)

1Math::NumSeq::AlgebraicUCsoenrtiCnounetdr(i3b)uted PerlMDaotchu:m:eNnutmaSteiqo:n:AlgebraicContinued(3)
2
3
4

NAME

6       Math::NumSeq::AlgebraicContinued -- continued fraction expansion of
7       algebraic numbers
8

SYNOPSIS

10        use Math::NumSeq::AlgebraicContinued;
11        my $seq = Math::NumSeq::AlgebraicContinued->new (expression => 'cbrt 2');
12        my ($i, $value) = $seq->next;
13

DESCRIPTION

15       This is terms in the continued fraction expansion of an algebraic
16       number such as a cube root or Nth root.  For example cbrt(2),
17
18           1, 3, 1, 5, 1, 1, 4, 1, 1, 8, 1, 14, 1, 10, 2, 1, 4, 12, 2, ...
19           starting i=0
20
21       A continued fraction approaches the root by a form
22
23                        1
24           C = a[0] + -------------
25                      a[1] +   1
26                             -------------
27                             a[2] +   1
28                                    ----------
29                                    a[3] + ...
30
31       The first term a[0] is the integer part of C, leaving a remainder
32       0 < r < 1 which is expressed as r=1/R with R > 1, so
33
34                      1
35          C = a[0] + ---
36                      R
37
38       Then a[1] is the integer part of that R, and so on repeatedly.
39
40       The current code uses a generic approach manipulating a polynomial with
41       "Math::BigInt" coefficients (see "FORMULAS" below).  It tends to be a
42       little slow because the coefficients become large, representing an ever
43       more precise approximation to the target value.
44
45   Expression
46       The "expression" parameter currently only accepts a couple of forms for
47       a cube root or Nth root.
48
49           cbrt 123
50           7throot 123
51
52       The intention would be to perhaps take some simple fractions or
53       products if they can be turned into a polynomial easily.  Or take an
54       initial polynomial directly.
55

FUNCTIONS

57       See "FUNCTIONS" in Math::NumSeq for behaviour common to all sequence
58       classes.
59
60       "$seq = Math::NumSeq::AlgebraicContinued->new (expression => $str)"
61           Create and return a new sequence object.
62
63       "$i = $seq->i_start ()"
64           Return 0, the first term in the sequence being at i=0.
65

FORMULAS

67   Next
68       The continued fraction can be developed by maintaining a polynomial
69       with single real root equal to the remainder R at each stage.  (As per
70       for example Knuth volume 2 Seminumerical Algorithms section 4.5.3
71       exercise 13, following Lagrange.)
72
73       As an example, a cube root cbrt(C) begins
74
75           -x^3 + C = 0
76
77       and later has a set of coefficients p,q,r,s
78
79           p*x^3 + q*x^2 + r*x + s = 0
80           p,q,r,s integers and p < 0
81
82       From such an equation the integer part of the root can be found by
83       looking for the biggest integer x with
84
85           p*x^3 + q*x^2 + r*x + s < 0
86
87       Choosing the signs so the high coefficient "p<0" means the polynomial
88       is positive for small x and becomes negative above the root.
89
90       Various root finding algorithms could probably be used, but the current
91       code is a binary search.
92
93       The integer part is subtracted R-c and then inverted 1/(R-c) for the
94       continued fraction.  This is applied to the cubic equation first by a
95       substitution x+c,
96
97           p*x^3 + (3pc+q)*x^2 + (3pc^2+2qc+r)x + (pc^3+qc^2+rc+s)
98
99       and then 1/x which is a reversal p,q,r,s -> s,r,q,p, and a term-wise
100       negation to keep p<0.  So
101
102           new p = -(p*c^3 + q*c^2 + r*c + s)
103           new q = -(3p*c^2 + 2q*c + r)
104           new r = -(3p*c + q)
105           new s = -p
106
107       The values p,q,r,s are integers but may become large.  For a cube root
108       they seem to grow by about 1.7 bits per term.  Presumably growth is
109       related to the average size of the a[i] terms.
110
111       For a general polynomial the substitution x+c becomes a set of binomial
112       factors for the coefficients.
113
114       For a square root or other quadratic equation q*x^2+rx+s the continued
115       fraction terms repeat and can be calculated more efficiently than this
116       general approach (see Math::NumSeq::SqrtContinued).
117
118       The binary search or similar root finding algorithm for the integer
119       part is important.  The integer part is often 1, and in that case a
120       single check to see if x=2 gives poly<0 suffices.  But a term can be
121       quite large so a linear search 1,2,3,4,etc is undesirable.  An example
122       with large terms can be found in Sloane's OEIS,
123
124           <http://oeis.org/A093876> continued fraction of 4th root of 9.1,
125           ie. (91/10)^(1/4)
126
127       The first few terms include 75656 and 262344, before settling down to
128       more usual size terms it seems.
129

HOME PAGE

134       <http://user42.tuxfamily.org/math-numseq/index.html>
135

LICENSE

137       Copyright 2012, 2013, 2014, 2016, 2019, 2020 Kevin Ryde
138
139       Math-NumSeq is free software; you can redistribute it and/or modify it
140       under the terms of the GNU General Public License as published by the
141       Free Software Foundation; either version 3, or (at your option) any
142       later version.
143
144       Math-NumSeq is distributed in the hope that it will be useful, but
145       WITHOUT ANY WARRANTY; without even the implied warranty of
146       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
147       General Public License for more details.
148
149       You should have received a copy of the GNU General Public License along
150       with Math-NumSeq.  If not, see <http://www.gnu.org/licenses/>.
151
152
153
154perl v5.36.0                      2022-07-22Math::NumSeq::AlgebraicContinued(3)