The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.
 / 
Math-NumSeq-72
River stage zero No dependents
3 ++
/ Math::NumSeq::AlgebraicContinued

NAME

Math::NumSeq::AlgebraicContinued -- continued fraction expansion of algebraic numbers

SYNOPSIS

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

DESCRIPTION

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[0] + -------------
               a[1] +   1
                      -------------
                      a[2] +   1
                             ----------
                             a[3] + ...

The first term a[0] 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[0] + ---
               R

Then a[1] 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

The 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.

FUNCTIONS

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.

FORMULAS

Next

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,

The first few terms include 75656 and 262344, before settling down to more usual size terms it seems.

SEE ALSO

Math::NumSeq, Math::NumSeq::SqrtContinued

HOME PAGE

http://user42.tuxfamily.org/math-numseq/index.html

LICENSE

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/>.