and 1 contributors

# 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 + -------------
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 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.