Math::NumSeq::AlgebraicContinued -- continued fraction expansion of algebraic numbers
use Math::NumSeq::AlgebraicContinued; my $seq = Math::NumSeq::AlgebraicContinued->new (expression => 'cbrt 2'); my ($i, $value) = $seq->next;
This is terms in the continued fraction expansion of an algebraic number such as a cube root or Nth root. For example cbrt(2),
1, 3, 1, 5, 1, 1, 4, 1, 1, 8, 1, 14, 1, 10, 2, 1, 4, 12, 2, ... starting i=0
A continued fraction approaches the root by a form
1 C = a + ------------- a + 1 ------------- a + 1 ---------- a + ...
The first term a is the integer part of C, leaving a remainder 0 < r < 1 which is expressed as r=1/R with R > 1, so
1 C = a + --- R
Then a is the integer part of that R, and so on repeatedly.
The current code uses a generic approach manipulating a polynomial with
Math::BigInt coefficients (see "FORMULAS" below). It tends to be a little slow because the coefficients become large, representing an ever more precise approximation to the target value.
expression parameter currently only accepts a couple of forms for a cube root or Nth root.
cbrt 123 7throot 123
The intention would be to perhaps take some simple fractions or products if they can be turned into a polynomial easily. Or take an initial polynomial directly.
See "FUNCTIONS" in Math::NumSeq for behaviour common to all sequence classes.
$seq = Math::NumSeq::AlgebraicContinued->new (expression => $str)
Create and return a new sequence object.
$i = $seq->i_start ()
Return 0, the first term in the sequence being at i=0.
The continued fraction can be developed by maintaining a polynomial with single real root equal to the remainder R at each stage. (As per for example Knuth volume 2 Seminumerical Algorithms section 4.5.3 exercise 13, following Lagrange.)
As an example, a cube root cbrt(C) begins
-x^3 + C = 0
and later has a set of coefficients p,q,r,s
p*x^3 + q*x^2 + r*x + s = 0 p,q,r,s integers and p < 0
From such an equation the integer part of the root can be found by looking for the biggest integer x with
p*x^3 + q*x^2 + r*x + s < 0
Choosing the signs so the high coefficient
p<0 means the polynomial is positive for small x and becomes negative above the root.
Various root finding algorithms could probably be used, but the current code is a binary search.
The integer part is subtracted R-c and then inverted 1/(R-c) for the continued fraction. This is applied to the cubic equation first by a substitution x+c,
p*x^3 + (3pc+q)*x^2 + (3pc^2+2qc+r)x + (pc^3+qc^2+rc+s)
and then 1/x which is a reversal p,q,r,s -> s,r,q,p, and a term-wise negation to keep p<0. So
new p = -(p*c^3 + q*c^2 + r*c + s) new q = -(3p*c^2 + 2q*c + r) new r = -(3p*c + q) new s = -p
The values p,q,r,s are integers but may become large. For a cube root they seem to grow by about 1.7 bits per term. Presumably growth is related to the average size of the a[i] terms.
For a general polynomial the substitution x+c becomes a set of binomial factors for the coefficients.
For a square root or other quadratic equation q*x^2+rx+s the continued fraction terms repeat and can be calculated more efficiently than this general approach (see Math::NumSeq::SqrtContinued).
The binary search or similar root finding algorithm for the integer part is important. The integer part is often 1, and in that case a single check to see if x=2 gives poly<0 suffices. But a term can be quite large so a linear search 1,2,3,4,etc is undesirable. An example with large terms can be found in Sloane's OEIS,
http://oeis.org/A093876 continued fraction of 4th root of 9.1, ie. (91/10)^(1/4)
The first few terms include 75656 and 262344, before settling down to more usual size terms it seems.
Copyright 2012, 2013, 2014 Kevin Ryde
Math-NumSeq is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3, or (at your option) any later version.
Math-NumSeq is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with Math-NumSeq. If not, see <http://www.gnu.org/licenses/>.