#############################################################################
##
## This file is part of GAP, a system for computational discrete algebra.
## This file's authors include Stefan Kohl.
##
## Copyright of GAP belongs to its developers, whose names are too numerous
## to list here. Please refer to the COPYRIGHT file for details.
##
## SPDX-License-Identifier: GPL-2.0-or-later
##
## This file contains code for computing (with) continued fraction
## expansions of real numbers.
##
#############################################################################
##
#F ContinuedFractionExpansionOfRoot( <P>, <n> )
##
InstallGlobalFunction( ContinuedFractionExpansionOfRoot,
function ( P, n )
local a, ai, Pi, Pi_1, pols, d, x0, step, bincoeff, i, j, k;
if not IsUnivariatePolynomial(P)
or not (IsPosInt(n) or DegreeOfLaurentPolynomial(P) = 2 and n = 0)
or not ForAll(CoefficientsOfLaurentPolynomial(P)[1],IsInt)
or LeadingCoefficient(P) < 0
then
Error("usage: ContinuedFractionExpansionOfRoot( <P>, <n> ) ",
"for a polynomial P with integer coefficients and a ",
"positive integer <n>");
fi;
P := CoefficientsOfLaurentPolynomial(P);
P := Concatenation(ListWithIdenticalEntries(P[2],0),P[1]);
d := Length(P) - 1;
bincoeff := List([0..d],n->List([0..d],k->Binomial(n,k)));
if ValuePol(P,0) >= 0 then
Error("the value of <P> at x = 0 has to be negative");
fi;
a := []; Pi := ShallowCopy(P); pols := []; i := 1;
while i <= n or n = 0 do
if d = 2 and n = 0 then
Add(pols,Pi);
fi;
x0 := 1; step := 1;
while ValuePol(Pi,x0) < 0 do
x0 := x0 + step;
step := 2 * step;
od;
step := step/4;
while step >= 1 do
if ValuePol(Pi,x0) > 0 then
x0 := x0 - step;
else
x0 := x0 + step;
fi;
step := step/2;
od;
if ValuePol(Pi,x0) > 0 then
ai := x0 - 1;
else
ai := x0;
fi;
a[i] := ai;
Pi_1 := ShallowCopy(Pi); Pi := ListWithIdenticalEntries(d+1,0);
for j in [1..d+1] do
for k in [1..d-j+2] do
Pi[k] := Pi[k] + Pi_1[d-j+2]*bincoeff[d-j+2][k]*ai^(d+2-j-k);
od;
od;
Pi := -Reversed(Pi);
if Pi[d+1] = 0 then
break; # Root is rational.
fi;
if d = 2 and n = 0 and Pi in pols then
break; # One period is done.
fi;
i := i + 1;
od;
return a;
end );
#############################################################################
##
#F ContinuedFractionApproximationOfRoot( <P>, <n> )
##
InstallGlobalFunction( ContinuedFractionApproximationOfRoot,
function ( P, n )
local a, M;
if not IsUnivariatePolynomial(P)
or not (IsPosInt(n) or DegreeOfLaurentPolynomial(P) = 2 and n = 0)
or not ForAll(CoefficientsOfLaurentPolynomial(P)[1],IsInt)
or LeadingCoefficient(P) < 0
then
Error("usage: ContinuedFractionApproximationOfRoot( <P>, <n> ) ",
"for a polynomial P with integer coefficients and a ",
"positive integer <n>");
fi;
a := ContinuedFractionExpansionOfRoot(P,n);
M := Product(a,a_i->[[a_i,1],[1,0]]);
return M[1][1]/M[2][1];
end );
