search.cpan.org is shutting down
Kevin Ryde > Math-NumSeq > Math::NumSeq::AlgebraicContinued

Math-NumSeq-72.tar.gz

Dependencies

Annotate this POD

Website

# CPAN RT

 Open 0
View/Report Bugs
Module Version: 72

# 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,

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.

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