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

NAME

Math::NumSeq::CbrtContinued -- continued fraction expansion of a cube root

SYNOPSIS

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

DESCRIPTION

In progress ... general algebraic irrational ?

This is terms in the continued fraction expansion of a cube root. For example cbrt(2),

    1, 3, 1, 5, 1, 1, 4, 1, 1, 8, 1, 14, 1, 10, 2, 1, 4, 12, 2, ...

The continued fraction approaches the root by

                      1   
   cbrt(C) = a[0] + ----------- 
                    a[1] +   1
                           -----------
                           a[2] +   1
                                  ----------
                                  a[3] + ...

The first term a[0] is the integer part of the root, leaving a remainder 0 < r < 1 which is expressed as r=1/R with R > 1, so

                     1   
   cbrt(C) = a[0] + ---
                     R

Then a[1] is the integer part of that R, and so on recursively.

The current code uses a fairly generic algebraic root approach and Math::BigInt coefficients (see "FORMULAS" below). The coefficients tend to be come large and the calculation a bit slow.

FUNCTIONS

See "FUNCTIONS" in Math::NumSeq for behaviour common to all sequence classes.

$seq = Math::NumSeq::CbrtContinued->new (cbrt => $s)

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 cubic equation with real root equal to the remainder R at each stage. Initially this is cbrt(C) so

    -x^3 + C = 0

and later in general

    p*x^3 + q*x^2 + r*x + s = 0
    p,q,r,s integers and p < 0

From such a cubic 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 p<0 means the cubic is positive for small x and becomes negative after the root.

Various root finding algorithms could probably be used, but the current code is just a binary search. The integer part of the remainder is often 1, so often it's enough to make a single check to see if x=2 gives cubic<0.

The integer part is subtracted R-a and inverted 1/(R-a) for the continued fraction. This is applied to the cubic equation by a substitution x+a,

    p*x^3 + (3pa+q)*x^2 + (3pa^2+2qa+r)x + (pa^3+qa^2+ra+s)

and then 1/x which is a reversal p,q,r,s -> s,r,q,p, then a term-wise negation to keep p<0. So

    new p = -(p*a^3 + q*a^2 + r*a + s)
    new q = -(3p*a^2 + 2q*a + r)
    new r = -(3p*a + q)
    new s = -p

The values p,q,r,s are integers but may become large, growing by about 1.7 bits per term generated. Presumably growth is related to the average size of the a[i] terms.

This same approach extends to any algebraic number, ie. root of a polynomial. But quadratics can be handled more easily.

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