Various algorithms for the continued fraction of a real algebraic number can be found in the literature, e.g. [1], [2], [3]. This page sketches an algorithm that uses only integer arithmetic, and shows some code (in Delphi Pascal) to execute it.
As usual, we suppose that the algebraic number is a simple root of a given polynomial with rational coefficients. We further assume w.l.o.g. that the polynomial is of the form
| ( | n | ) | bnxn + | ( | n | ) | bn–1xn–1 + | ( | n | ) | bn–2xn–2 + ··· + | ( | n | ) | b0 , |
| 0 | 1 | 2 | n |
where the weighted coefficients bi are integers.
As in the algorithm of Rosen & Shallit [2], only integer arithmetic is used, so that rounding errors cannot occur. For algebraic numbers of degree > 2, however, the integers in the working grow without limit as the algorithm proceeds. This seems inevitable. An algorithm that used only a finite amount of memory would have only finitely many states, and its output of partial quotients would therefore either terminate or go into an infinite loop. This can happen only if the root is either rational, or irrational of degree 2.
As in [2], therefore, the user must supply arithmetic routines that can handle integers of arbitrary size. This should not be too difficult, since the algorithm needs only the methods listed at the start of the code section below. Note that the algorithm does not need multiplication or division, except by 2; these operations could be effected by left and right shift respectively.
When setting up the coefficients at the start of the algorithm, the user can also set up the start of the continued fraction as an array of partial quotients. This allows the user to define which root of the polynomial is required, if the polynomial has more than one real root.
How to determine the number of real roots, their multiplicities, and their approximate values is a separate problem, which is not treated here; see [1] or [2] for a discussion.
This section summarizes the working of the sample code (below). Suppose for simplicity that the polynomial is of degree n = 3. The array c[] holds the weighted coefficients of a polynomial
| F(x) = c0 + 3c1x + 3c2x2 + c3x3 , |
which is initialized to the user’s polynomial F0(x) by setting cj = bj for j = 0, …, n. The method StepForward applies n(n – 1)/2 additions with the effect
| (c0, c1, c2, c3) := (c0 + 3c1 + 3c2 + c3, c1 + 2c2 + c3, c2 + c3, c3) . |
It is easy to check that the new values are the weighted coefficients of F(x + 1). In particular, the new c0 equals F(1). By repeatedly applying StepForward we calculate F(2), F(3), etc.
This process is speeded up by calling DoubleStep, which applies n(n – 1)/2 doublings with the effect
| (c0, c1, c2, c3) := (c0, 2c1, 4c2, 8c3) . |
The new values are the weighted coefficients of F(2x), and StepForward now advances the argument of F by 2 each time. By calling DoubleStep again the step is increased to 4, and so on.
The methods StepBackward and HalveStep are the inverses of StepForward and DoubleStep respectively.
Assuming that F has a unique positive root, we can thus apply a binary search to find the partial quotient a0, as the integer for which either F0(a0) = 0 or F0(a0) and F0(a0 + 1) have opposite signs. The case F0(a0) = 0 is detected by c0 = 0 and the algorithm terminates. Otherwise, we call ReverseCoefficients, which has the effect
| (c0, c1, c2, c3) := (c3, c2, c1, c0) |
The new coefficients define a polynomial
| F1(x) = x3F0(a0 + 1/x). |
By applying the above process to F1 we can find a1, the greatest integer for which F(a0 + 1/a1) has the same sign as F(a0 + 1). Similarly we can find a2, a3, etc. as far as desired (or until a rational root is detected by c0 = 0).
The following Pascal code is part of a Delphi program that has been checked against the examples in references [1,2,3,4].
The method Initialize, besides initializing the array c[], allows the user to input the first few partial quotients in order to specify which root of the poynomial is required. The method GetNextPartialQuotient is then called repeatedly; it returns false if the root is rational and all its partial quotients have already been found.
{===============================================================================
This Delphi code uses a unit UnlimitedInt (not shown here) which defines
the types TUnlimitedInt and TOrdinaryInt.
TUnlimitedInt holds signed integers of arbitrary size.
TOrdinaryInt is used for coefficients and partial quotients, and can be an alias
for a signed integer type of the programming language, e.g. longint in Delphi.
The unit UnlimitedInt should provide the following methods:
function Sign( var u : TUnlimitedInt) : integer; // -1, 0, or +1 as usual
procedure Double( var u : TUnlimitedInt); // u := u*2;
procedure Halve( var u: TUnlimitedInt); // u := u/2
procedure Copy( var u, v : TUnlimitedInt); // u := v
procedure Add( var u, v : TUnlimitedInt); // u := u + v
procedure Sub( var u, v : TUnlimitedInt); // u := u - v
procedure Convert( var u : TUnlimitedInt;
i : TOrdinaryInt); // converts i to TUnlimitedInt
===============================================================================}
var
n : integer; // degree of polynomial
a : TOrdinaryInt; // partial quotient
d : TOrdinaryInt; // step used in building partial quotient
c : array of TUnlimitedInt; // weighted coefficients of current polynomial
a_user : array of TOrdinaryInt; // partial quotients set by user (if any)
n_user : integer; // length of array a_user (can be 0)
index : integer; // index of current partial quotient
foundRational : boolean; // set true if root is found to be rational
{-------------------------------------------------------------------------------
Initialize the algorithm.
(1) Set up array of weighted coefficients. Polynomial has degree n and is
(n|0)*b[n]*x^n + (n|1)*b[n-1]*x^(n-1) + (n|2)*b[n-2]*x^(n-2) + ... + (n|n)*b[0]
where (n|r) is the usual binomial coefficient "n choose r".
(2) Optionally set up the first few partial quotients, in order to specify
which root of the polynomial is required.
Partial quotients are in the array a, which is allowed to be empty.
}
procedure Initialize( b : array of TOrdinaryInt;
a : array of TOrdinaryInt);
var
degree, j : integer;
inputOK : boolean;
begin
// Validate input (sketch).
// A full version would provide more helpful error messages.
inputOK := true;
degree := Length(b) - 1;
if (degree < 1) then inputOK := false; // degree must be at least 1
if (b[degree] = 0) then inputOK := false; // leading coefficient can't be 0
j := Length(a);
while (j > 1) and inputOK do begin
dec(j);
if (a[j] <= 0) then inputOK := false; // must have a[j] > 0 if j > 0
end;
if not inputOK then begin
raise Exception.Create( 'Bad input'); // abort the procedure somehow
end;
// If OK, set up the polynomial and partial quotients
n := degree;
SetLength( c, n + 1);
for j := 0 to n do UnlimitedInt.Convert( c[j], b[j]);
n_user := Length(a);
SetLength( a_user, n_user);
for j := 0 to n_user - 1 do a_user[j] := a[j];
// Reset the algorithm
index := 0;
foundRational := false;
end;
{-------------------------------------------------------------------------------
Method to return the next partial quotient, after initializing the algorithm.
The partial quotients returned include those that were set up initially.
Function return is true if partial quotient is valid, false if no more
partial quotients exist (rational root).
}
function GetNextPartialQuotient( var quotient : TOrdinaryInt) : boolean;
var
a_given : TOrdinaryInt; // partial quotient supplied by user
a_work : TOrdinaryInt; // shifted version of a_given
testSign : integer; // used to detect change of sign
begin
if foundRational then begin
// If here, root was found rational on a previous call.
quotient := 0; // return some arbitrary value
result := false; // tell caller no more partial quotients
end
else if (index < n_user) then begin
// If here, we're still returning the user's partial quotients.
// We need to update the working array c[] neverheless.
a_given := a_user[index];
a_work := Abs(a_given);
d := 1;
while (a_work > 0) do begin
if Odd( a_work) then begin
if (a_given > 0) then StepForward() else StepBackward();
end;
a_work := a_work div 2;
if (a_work > 0) then DoubleStep();
end;
while (d > 1) do HalveStep();
quotient := a_given;
result := true; // got a partial quotient
end
else begin
// If here, we need to calculate the next partial quotient.
// Increment partial quotient till c[0] changes sign
a := 0;
d := 1;
testSign := UnlimitedInt.Sign( c[0]);
StepForward();
if (index > 0) then testSign := UnlimitedInt.Sign( c[0]);
while (UnlimitedInt.Sign(c[0]) = testSign) do begin
if (a > 1) then DoubleStep();
StepForward();
end;
// Determine partial quotient by binary search
while (d > 1) do begin
HalveStep();
if (UnlimitedInt.Sign(c[0]) = 0)
or (UnlimitedInt.Sign(c[0]) = testSign) then StepForward() else StepBackward();
end;
if (UnlimitedInt.Sign(c[0]) = -testSign) then StepBackward();
quotient := a;
result := true; // got a partial quotient
// Test for rational root; if so, set flag for future calls
if (UnlimitedInt.Sign(c[0]) = 0) then foundRational := true;
end;
// If not yet found rational root, update the algorithm
if not foundRational then begin
ReverseCoefficients();
inc( index);
end;
end;
{===============================================================================
Local subroutines.
-------------------------------------------------------------------------------}
procedure StepForward();
var
j, k : integer;
begin
a := a + d;
for k := 1 to n do
for j := 0 to n - k do
UnlimitedInt.Add( c[j], c[j+1]);
end;
{------------------------------------------------------------------------------}
procedure StepBackward();
var
j, k : integer;
begin
a := a - d;
for k := 1 to n do
for j := 0 to n - k do
UnlimitedInt.Sub( c[j], c[j+1]);
end;
{------------------------------------------------------------------------------}
procedure DoubleStep();
var
j, k : integer;
begin
d := d*2;
for k := 1 to n do
for j := 1 to k do
UnlimitedInt.Double( c[k]);
end;
{------------------------------------------------------------------------------}
procedure HalveStep();
var
j, k : integer;
begin
d := d div 2;
for k := 1 to n do
for j := 1 to k do
UnlimitedInt.Halve( c[k]);
end;
{------------------------------------------------------------------------------}
procedure ReverseCoefficients();
var
j : integer;
c_temp : UnlimitedInt.TUnlimitedInt;
begin
for j := 0 to (n div 2) do begin
UnlimitedInt.Copy( c_temp, c[j]);
UnlimitedInt.Copy( c[j], c[n-j]);
UnlimitedInt.Copy( c[n-j], c_temp);
end;
end;
| [1] | David G. Cantor, Paul H. Galyean & Horst G. Zimmer, “A continued fraction algorithm for real algebraic numbers”, Mathematics of Computation, 26 (119), July 1972, 785–791. |
| [2] | David Rosen & Jeffrey Shallit, “A continued fraction algorithm for approximating all real polynomial roots”, Mathematics Magazine, 51 (2), March 1978, 112–116. |
| [3] | Enrico Bombieri & Alfred J. van der Poorten, “Continued fractions of algebraic numbers”, in Wieb Bosma & Alf van der Poorten (eds), Computational Algebra and Number Theory, Kluwer, 1995, 137–152. |
| [4] | R.F. Churchhouse & S.T.E. Muir, “Continued fractions, algebraic numbers and modular invariants”, Journal of the Institute of Mathematics and its Applications, 5, 1969, 318–328. |