Package $approx -- find simple approximations to rationals
-- Export
Export FloatApprox;
//Export CFApprox; now BI
Export CFApproximants;
Export ContFrac;
Export ContFracToRat;


Define About()
  PrintLn "    Author: J Abbott";
  PrintLn "    Date: 24 May 2003";
  PrintLn "    Comment: first version, efficiency not a strong point.";
EndDefine; -- About
-----------------------------------------------------------------------------
-- The exported functions of this package are:
--  FloatApprox    finds an approximation of the form M*10^E within a
--                 maximum specified relative error.
--  CFApprox       finds the "simplest" rational approximation within a
--                 maximum specified relative error.
--  CFApproximants returns a list of all continued fraction approximants
--                 to a specified rational
--  ContFrac       returns a list of the continued fraction "denominators"
--                 for a given rational number.

Define FloatApprox(N, RelativeError)
  If type(N) <> INT And type(N) <> RAT Then error("FloatApprox: first argument must be rational"); EndIf;
  If type(RelativeError) <> INT And type(RelativeError) <> RAT Then
    error("FloatApprox: second argument must be rational");
  EndIf;
  If RelativeError <= 0 Or RelativeError > 0.1 Then RelativeError := 0.1; EndIf;
  If N = 0 Then Return 0; EndIf;
  Ndigits := 1;
  Power := 1/5;
  While Power > RelativeError Do
    Power := Power/10;
    Ndigits := Ndigits+1;
  EndWhile;
  Tmp := MantissaAndExponent(N, Ndigits);
  Return Tmp.Mantissa*10^(Tmp.Exponent - Ndigits + 1);
EndDefine; -- FloatApprox


Define ContFrac(N)
  If type(N) = INT Then Return [N]; EndIf;
  If type(N) <> RAT Then error("ContFrac: argument must be rational"); EndIf;
  If N = 0 Then Return [0]; EndIf; -- We define this result for zero.
  D := den(N);
  N := num(N);
  Ans := [];
  While D <> 0 Do
    Quot := div(N, D);
    append(Ref Ans, Quot);
    N := N - Quot*D;
    Swap := N; N := D; D := Swap;
  EndWhile;
  Return Ans;
EndDefine; -- ContFrac


Define ContFracToRat(L)
  If type(L) <> LIST Or MakeSet(shape(L)) <> [INT] Then
    error("ContFracToRat: argument should be a list of (positive) integer \"quotients\"");
  EndIf;
  N1 := 0; N2 := 1;
  D1 := 1; D2 := 0;
  Foreach Q In L Do
    If D2 <> 0 And Q < 1 Then error("ContFracToRat: negative continued fraction quotient"); EndIf;
    N1 := N1 + Q*N2;
    Swap := N1; N1 := N2; N2 := Swap;
    D1 := D1 + Q*D2;
    Swap := D1; D1 := D2; D2 := Swap;
  EndForeach;
  Return N2/D2;
EndDefine; -- ContFracToRat
  

Define CFApprox(N, RelativeError)
  If type(N) <> INT And type(N) <> RAT Then error("CFApprox: first argument must be rational"); EndIf;
  If type(RelativeError) <> INT And type(RelativeError) <> RAT Then
    error("CFApprox: second argument must be rational");
  EndIf;
  If RelativeError <= 0 Or RelativeError > 0.1 Then RelativeError := 0.1; EndIf;
  If N = 0 Then Return N; EndIf;
  AbsError := abs(N)*RelativeError;
  RecipError := 1+div(den(AbsError), num(AbsError));
  CF := ContFrac(N);
  N1 := 0; N2 := 1;
  D1 := 1; D2 := 0;
  Foreach Denom In CF Do
    If D2*Denom > RecipError Then Break; EndIf;
    N1 := N1 + Denom*N2;
    Swap := N1; N1 := N2; N2 := Swap;
    D1 := D1 + Denom*D2;
    Swap := D1; D1 := D2; D2 := Swap;
  EndForeach;
  Return N2/D2;
EndDefine; -- CFApprox


Define CFApproximants(N)
  If type(N) = INT Then Return [N]; EndIf;
  If N = 0 Then Return [0]; EndIf;
  If type(N) <> RAT Then error("CFApproximants: argument must be rational"); EndIf;
  CF := ContFrac(N);
  N1 := 0; N2 := 1;
  D1 := 1; D2 := 0;
  Ans := [];
  Foreach Denom In CF Do
    N1 := N1 + Denom*N2;
    Swap := N1; N1 := N2; N2 := Swap;
    D1 := D1 + Denom*D2;
    Swap := D1; D1 := D2; D2 := Swap;
    append(Ref Ans, N2/D2);
  EndForeach;
  Return Ans;
EndDefine; -- CFApproximants

EndPackage; -- end of package $approx
