The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.
 / 
Math-NumSeq-54
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, ...

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 since the coefficients become large, representing an every more precise approximation to the target value.

Expression

The expression parameter currently only accepts a couple of forms for a cube root or general specified 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 cofficient 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 4throot 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.

SEE ALSO

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

HOME PAGE

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

LICENSE

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