Package $float  -- convert rationals to floating-point decimal strings

Export DecimalStr;
Export FloatStr;
Export MantissaAndExponent;

-----------------------------------------------------------------------------
--  Function to convert a rational to a decimal string of the form
--  -d.ddddd*10^-ddd; the optional second argument specifies how
-- many digits to print in the mantissa (min 2 digits).  There is no limit
-- on the exponent, but the code gets slow for very large or very
-- small inputs.

Define FloatStr(X, Opt N)
  ErrMesg := "FloatStr: arg must be rational, with optional second integer arg";
  BadArgs := False;
  If IsDefined(N) And type(N) <> INT Then error(ErrMesg); EndIf;
  If not(type(X) IsIn [INT,RAT] Or IsRational(Ref X,X))
    Then error(ErrMesg);
  EndIf;
  If not(IsDefined(N)) Then SigFig := 10;
  Elif N < 2 Then SigFig := 2;
  Else SigFig := N;
  EndIf;
  Sign := "";
  If X < 0 Then Sign := "-"; X := -X; EndIf;
  Rec := $.MantissaAndExponent(X, SigFig);
  If Rec.mantissa <> 0 Then Mant := Sprint(Rec.mantissa);
  Else Mant := sum(NewList(SigFig, "0")); EndIf;
  If Rec.exponent < 0 Then Expo := "*10^("+Sprint(Rec.exponent)+")";
  Else Expo := "*10^"+Sprint(Rec.exponent); EndIf;
  Return Sign+Mant[1]+"."+sum([Mant[I] | I In 2..len(Mant)])+Expo;
EndDefine;


-- This function returns a Record[Mantissa, Exponent].
-- If X = 0 then both components in the record are set to 0.
-- Otherwise Mantissa is an integer with exactly SigFig decimal digits,
-- and Exponent is that integer such that 1 <= abs(X)*10^Exponent < 10

Define MantissaAndExponent(X, SigFig)
  ErrMesg := "MantissaAndExponent: first arg must be a rational, and second an integer";
  If type(SigFig) <> INT Then error(ErrMesg); EndIf;
  If not(type(X) IsIn [INT,RAT] Or IsRational(Ref X,X))
    Then error(ErrMesg);
  EndIf;
  If SigFig < 1 Then SigFig := 1; EndIf;
  If X = 0 Then Return Record[mantissa:=0, exponent:=0]; EndIf;
  Sign := 1;
  If X < 0 Then Sign := -1; X := -X; EndIf;
  N := num(X);
  D := den(X); -- we are guaranteed that D and N are coprime
  If N >= D Then
    Exponent := $float.ILogBase(div(N, D), 10);
    D := D*10^Exponent;
  Else
    Exponent := - $float.ILogBase(div(D, N), 10);
    Power10 := 10^(-Exponent);
    If N = 1 And D = Power10 Then -- X = 1/10^n is a special case
      D := 1;
    Else
      Exponent := Exponent - 1;
      N := 10*N*Power10;
    EndIf;
  EndIf;
  Scale := 10^(SigFig-1);
  -- the next line uses wastefully high precision in most cases...
  Mantissa := div(Scale*N+div(D, 2), D);
  If Mantissa >= 10*Scale Then -- Mantissa was rounded up to next "decade"
    Exponent := Exponent+1;
    Mantissa := div(Mantissa, 10);
  EndIf;
  Mantissa := Sign*Mantissa;
  Return Record[mantissa:=Mantissa, exponent:=Exponent];
EndDefine;


Define ILogBase(N, Base)
  If not(type(N) IsIn [INT, RAT]) Or type(Base) <> INT Then
    error("ILogBase: first arg must be integer or rational, second arg must be integer");
  EndIf;
  If N = 0 Then error("ILogBase: cannot compute logarithm of 0"); EndIf;
  If Base < 0 Then Base := -Base; EndIf;
  If Base < 2 Then error("ILogBase: base must be at least 2"); EndIf;
  If N < 0 Then N := -N; EndIf;
  If type(N) = RAT Then
    If num(N) >= den(N) Then
      N := div(num(N), den(N));
    Else
      Tmp := ILogBase(1/N, Base);
      If num(N) = 1 And den(N) <> Base^Tmp Then Tmp := Tmp+1; EndIf;
      Return -Tmp;
    EndIf;
  EndIf;
  CurrPower := Base;
  Powers := [];
  Exponent := 1;
  While N >= CurrPower Do
    append(ref Powers, CurrPower);
    CurrPower := CurrPower^2;
    Exponent := 2*Exponent;
  EndWhile;
  reverse(ref Powers);
  LogN := 0;
  Foreach Power In Powers Do
    Exponent := Exponent/2;
    If N >= Power Then
      N := div(N, Power);
      LogN := LogN + Exponent;
    EndIf;
  EndForeach;
  Return LogN;
EndDefine; -- ILogBase


-----------------------------------------------------------------------------

--  Function to convert a rational to a decimal string of the form
--  ddddd.ddd; the optional second argument specifies how many decimal
--  digits to print (default 3 digits).  There is no limit on the
--  precision of the integer part.

Define DecimalStr(X, Opt SigFig)
  If not(IsDefined(SigFig)) Then SigFig := 3; EndIf;
  ErrMesg := "DecimalStr: arg must be rational, with optional second integer arg";
  If not(type(X) IsIn [INT,RAT,RINGELEM]) Or type(SigFig) <> INT Then
    error(ErrMesg);
  EndIf;
  If type(X) = RINGELEM And not(IsRational(Ref X,X)) Then
    error(ErrMesg);
  EndIf;
  Rec := $float.IntegerAndDecimal(X, SigFig);
  If Rec.Sign = -1 Then
    Return "-" + Sprint(Rec.Integer) + Rec.Decimal;
  EndIf;
  Return Sprint(Rec.Integer) + Rec.Decimal;
EndDefine; -- DecimalStr


-- This function returns a Record[Integer, Decimal].
-- If X = 0 then both components in the record are set to 0.
-- Otherwise Integer is the integer part given as an integer,
-- and Decimal is the truncated decimal part given as a string which
--  . If X = Integer : is the empty string
--  . if X = Integer+Decimal : contains the dot and <=SigFig digits
--  . contains the dot and =SigFig digits otherwise

Define IntegerAndDecimal(X, SigFig)
  ErrMesg := "IntegerAndDecimal: first arg must be a rational, and second an integer";
  If type(SigFig) <> INT Then error(ErrMesg); EndIf;
  If not(type(X) IsIn [INT,RAT] Or IsRational(Ref X,X))
    Then error(ErrMesg);
  EndIf;
  Rec := Record[Sign:=0, Integer:=0, Decimal:=""];
  If X = 0 Then Return Rec; EndIf;
  If X > 0 Then Rec.Sign := 1; Else Rec.Sign := -1; X := -X; EndIf;
  N := num(X);
  D := den(X);
  Rec.Integer := div(N,D);
  N := mod(N, D);
  If N = 0 Or SigFig <=0 Then Return Rec; EndIf;
  Rec.Decimal := ".";
  For Digits := 1 To SigFig Do
    N := N*10;
    Rec.Decimal := Rec.Decimal + Sprint(div(N, D));
    N := mod(N, D);
    If N = 0 Then Return Rec; EndIf;
  EndFor;
  Return Rec;
EndDefine; -- IntegerAndDecimal

-----------------------------------------------------------------------------

-- This is a very simplistic implementation!
Define atan(N, Digits) -- computes 10^Digits * atan(1/N)
  ans := 0;
  summand := div(10^Digits, N);
  I := 1;
  sign := 1;
  While summand > 0 Do
    ans := ans + div(summand, I)*sign;
    sign := -sign;
    I := I+2;
    summand := div(summand, N*N);
  EndWhile;
  Return ans;
EndDefine; -- Atan

Define Pi(N)
  -- Better to use 4*(2*atan(3/79) + 3*atan(1/7))
  -- Or even 4*(5*atan(3/79) + 7*atan(29/278))
  numer := 4*(4*atan(5, N+5) - atan(239, N+5));
  Return numer/10^(N+5);
EndDefine; -- Pi


EndPackage; -- end of package float.
