/usr/share/gap/lib/contfrac.gi is in gap-libs 4r7p5-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 | #############################################################################
##
#W contfrac.gi Stefan Kohl
##
#Y Copyright (C) 2004 The GAP Group
##
##
## 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> )\n",
"for a polynomial P with integer coefficients and a ",
"positive integer <n>.\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.\n"); 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; fi; # Root is rational.
if d = 2 and n = 0 and Pi in pols then break; fi; # One period is done.
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> )\n",
"for a polynomial P with integer coefficients and a ",
"positive integer <n>.\n");
fi;
a := ContinuedFractionExpansionOfRoot(P,n);
M := Product(a,a_i->[[a_i,1],[1,0]]);
return M[1][1]/M[2][1];
end );
#############################################################################
##
#E contfrac.gi . . . . . . . . . . . . . . . . . . . . . . . . . . ends here
|