#############################################################################
##
#W  contfrac.gi                                                   Stefan Kohl
##
#Y  Copyright (C) 2004 The GAP Group
##
#H  @(#)$Id: contfrac.gi,v 4.5.2.1 2005/05/03 13:44:09 stefan Exp $
##
##  This file contains code for computing (with) continued fraction
##  expansions of real numbers.
##
Revision.contfrac_gi :=
    "@(#)$Id: contfrac.gi,v 4.5.2.1 2005/05/03 13:44:09 stefan Exp $";

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