(*********************************************************************)
(*                                                                   *)
(* Module ExMathLib0 Copyright © 1995 by Computer Inspirations       *)
(*                                                                   *)
(* Design : Michael Griebling                                        *)
(* Change : Original                                                 *)
(*                                                                   *)
(*********************************************************************)

MODULE ExMathLib0;

IMPORT LR := LongRealConversions, LM := MathIEEEDoubTrans,
       X := ExNumbers;

VAR
  ToRadians : X.ExNumType;
  ToDegrees : X.ExNumType;
  Fact500   : X.ExNumType;
  Fact1000  : X.ExNumType;
  Fact2000  : X.ExNumType;
  Fact3000  : X.ExNumType;


PROCEDURE ExNumToLongReal*(x : X.ExNumType) : LONGREAL;
VAR
  Num : LONGREAL;
  Str : ARRAY 81 OF CHAR;
BEGIN
  (* Convert ExNum into LONGREAL via a string *)
  X.ExNumToStr(x, 0, 0, Str);
  IF LR.StringToReal(Str, Num) THEN
    RETURN Num;
  ELSE
    RETURN 0.0D;
  END;
END ExNumToLongReal;


PROCEDURE LongRealToExNum*(x : LONGREAL; VAR Result : X.ExNumType);
VAR
  Str : ARRAY 81 OF CHAR;
BEGIN
  (* Convert LONGREAL into an ExNum via a string *)
  IF LR.RealToString(x, Str, 1, 52, TRUE) THEN
    X.StrToExNum(Str, Result);
  ELSE
    Result := X.Ex0;
  END;
END LongRealToExNum;


PROCEDURE xtoi*(VAR Result : X.ExNumType; x : X.ExNumType; i : LONGINT);
(* From Knuth, slightly altered : p442, The Art Of Computer Programming, Vol 2 *)
VAR
  Y : X.ExNumType;
  negative : BOOLEAN;
BEGIN
  Y := X.Ex1;
  negative := i < 0;
  i := ABS(i);
  LOOP
    IF ODD(i) THEN X.ExMult(Y, Y, x) END;
    i := i DIV 2;
    IF i = 0 THEN EXIT END;
    X.ExMult(x, x, x);
  END;
  IF negative THEN
    X.ExDiv(Result, X.Ex1, Y);
  ELSE
    Result := Y;
  END;
END xtoi;


PROCEDURE Root *(VAR Result : X.ExNumType;
                    x      : X.ExNumType;
                    i      : LONGINT);
(* Use iterative solution of a general root equation *)
VAR
  y, yp, f, g, t : X.ExNumType;
  iteration : INTEGER;
  root : LONGREAL;
  negate : BOOLEAN;
BEGIN
  IF ((x.Sign = X.negative) & ~ODD(i)) OR (i < 2) THEN
    X.ExStatus := X.IllegalNumber;
    Result := X.Ex0;
  ELSIF X.IsZero(x) THEN
    Result := x;
  ELSE
    (* handle negative roots *)
    IF x.Sign = X.negative THEN X.ExAbs(x); negate := TRUE
    ELSE negate := FALSE
    END;

    (* estimate of the ith root *)
    root := 1.0D / i;
    LongRealToExNum(LM.Pow(root,ExNumToLongReal(x)), yp);
    X.ExNumb(i, 0, 0, f);    (* i *)
    X.ExNumb(i-1, 0, 0, g);  (* i - 1 *)

    (* calculate the root *)
    iteration := 4;
    LOOP
      (* y := 1/i * (yp * (i-1) + x / yp^(i-1)) *)
      xtoi(t, yp, i-1);       (* yp**(i-1) *)
      X.ExMult(y, t, yp);     (* yp**i *)
      X.ExMult(y, y, g);      (* yp**i * (i-1) *)
      X.ExAdd(y, y, x);       (* yp**i * (i-1) + x *)
      X.ExMult(t, t, f);      (* yp**(i-1) * i *)
      X.ExDiv(y, y, t);
      IF (X.ExCompare(y, yp) = X.ExEqual) OR (iteration = 0) THEN EXIT END;
      DEC(iteration);
      yp := y;
    END;

    (* adjust the number's sign *)
    Result := y;
    IF negate THEN X.ExChgSign(Result) END;
  END;
END Root;


PROCEDURE powerof10(VAR Result : X.ExNumType; x : LONGINT);
BEGIN
  X.ExNumb(1, 0, SHORT(x), Result);
END powerof10;


PROCEDURE RadToDegX*(VAR radianAngle : X.ExNumType);
(* Convert a radian measure into degrees *)
BEGIN
  X.ExMult(radianAngle, ToDegrees, radianAngle);
END RadToDegX;


PROCEDURE DegToRadX*(VAR radianAngle : X.ExNumType);
(* Convert a degree measure into radians *)
BEGIN
  X.ExMult(radianAngle, ToRadians, radianAngle);
END DegToRadX;


PROCEDURE sqrtX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  Root(Result, x, 2);
END sqrtX;


PROCEDURE lnX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Log(ExNumToLongReal(x)), Result);
END lnX;


PROCEDURE logX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Log10(ExNumToLongReal(x)), Result);
END logX;


PROCEDURE factorial(VAR prevn, currentn : LONGINT;
                    VAR PrevFact, Result : X.ExNumType);
(* Implements an incremental factorial using a previously calculated value. *)
VAR
  i : LONGINT;
BEGIN
  FOR i := prevn+1 TO currentn DO
    (* PrevFact := PrevFact * i; *)
    X.ExNumb(i, 0, 0, Result);
    X.ExMult(PrevFact, PrevFact, Result);
  END;
  prevn := currentn;
  Result := PrevFact;
END factorial;


PROCEDURE factorialX*(VAR Result : X.ExNumType; n : LONGINT);
CONST
  MaxFactorial = 3249;
VAR
  fact : LONGINT;
  prev : X.ExNumType;
BEGIN
  IF (n < 0) OR (n > MaxFactorial) THEN
    X.ExStatus := X.IllegalNumber;
    Result := X.Ex0;
    RETURN;
  END;
  IF    n < 500  THEN prev := X.Ex1;      fact := 0
  ELSIF n < 1000 THEN prev := Fact500;  fact := 500
  ELSIF n < 2000 THEN prev := Fact1000; fact := 1000
  ELSIF n < 3000 THEN prev := Fact2000; fact := 2000
  ELSE                prev := Fact3000; fact := 3000
  END;
  factorial(fact, n, prev, Result);
END factorialX;


PROCEDURE expX*(VAR Result : X.ExNumType; x : X.ExNumType);
VAR
  xPower : LONGREAL;
BEGIN
  xPower := ExNumToLongReal(x);
  X.ExFrac(x);
  IF (ABS(xPower) < MAX(LONGINT)) & X.IsZero(x) THEN
    xtoi(Result, X.e, ENTIER(xPower));
  ELSE
    LongRealToExNum(LM.Exp(xPower), Result);
  END;
END expX;


PROCEDURE powerX*(VAR Result : X.ExNumType; x, y : X.ExNumType);
VAR
  yPower : LONGREAL;
BEGIN
  yPower := ExNumToLongReal(y);
  X.ExFrac(y);
  IF (ABS(yPower) < MAX(LONGINT)) & X.IsZero(y) THEN
    xtoi(Result, x, ENTIER(yPower));
  ELSE
    LongRealToExNum(LM.Pow(yPower,ExNumToLongReal(x)),Result);
  END;
END powerX;


PROCEDURE rootX*(VAR Result : X.ExNumType; x, y : X.ExNumType);
VAR
  yRoot : LONGREAL;
BEGIN
  yRoot := ExNumToLongReal(y);
  X.ExFrac(y);
  IF (ABS(yRoot) < MAX(LONGINT)) & X.IsZero(y) THEN
    Root(Result, x, ENTIER(yRoot));
  ELSE
    yRoot := 1.0D / yRoot;
    LongRealToExNum(LM.Pow(yRoot,ExNumToLongReal(x)),Result);
  END;
END rootX;


PROCEDURE sinX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Sin(ExNumToLongReal(x)), Result);
END sinX;


PROCEDURE cosX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Cos(ExNumToLongReal(x)), Result);
END cosX;


PROCEDURE tanX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Tan(ExNumToLongReal(x)), Result);
END tanX;


PROCEDURE arctanX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Atan(ExNumToLongReal(x)), Result);
END arctanX;


PROCEDURE coshX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Cosh(ExNumToLongReal(x)), Result);
END coshX;


PROCEDURE sinhX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Sinh(ExNumToLongReal(x)), Result);
END sinhX;


PROCEDURE tanhX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Tanh(ExNumToLongReal(x)), Result);
END tanhX;


PROCEDURE arccoshX*(VAR Result : X.ExNumType; x : X.ExNumType);
VAR
  Temp : X.ExNumType;
BEGIN
  (* Result = ln(x + sqrt(x*x - 1)) *)
  X.ExMult(Temp, x, x);
  X.ExSub(Temp, Temp, X.Ex1);
  sqrtX(Temp, Temp);
  X.ExAdd(Temp, x, Temp);
  lnX(Result, Temp);
END arccoshX;


PROCEDURE arcsinhX*(VAR Result : X.ExNumType; x : X.ExNumType);
VAR
  Temp : X.ExNumType;
BEGIN
  (* Result = ln(x + sqrt(x*x + 1)) *)
  X.ExMult(Temp, x, x);
  X.ExAdd(Temp, Temp, X.Ex1);
  sqrtX(Temp, Temp);
  X.ExAdd(Temp, x, Temp);
  lnX(Result, Temp);
END arcsinhX;


PROCEDURE arctanhX*(VAR Result : X.ExNumType; x : X.ExNumType);
VAR
  Temp, Temp2 : X.ExNumType;
BEGIN
  (* Result = ln((1 + x) / (1 - x)) / 2 *)
  X.ExAdd(Temp, X.Ex1, x);
  X.ExSub(Temp2, X.Ex1, x);
  X.ExDiv(Temp, Temp, Temp2);
  lnX(Result, Temp);
  X.ExNumb(0, 5, 0, Temp);
  X.ExMult(Result, Result, Temp);
END arctanhX;


PROCEDURE arcsinX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  LongRealToExNum(LM.Asin(ExNumToLongReal(x)), Result);
END arcsinX;


PROCEDURE arccosX*(VAR Result : X.ExNumType; x : X.ExNumType);
BEGIN
  (* Replacement algorithm *)
  LongRealToExNum(LM.Acos(ExNumToLongReal(x)), Result);
END arccosX;


BEGIN
  (* Initialize a few internal conversion constants *)
  X.StrToExNum(
  "5.729577951308232087679815481410517033240547246656420E+1",
  ToDegrees);
  X.StrToExNum(
  "1.745329251994329576923690768488612713442871888541727E-2",
  ToRadians);

  (* Speed up very large factorials *)
  X.StrToExNum(
  "1.220136825991110068701238785423046926253574342803193E+1134",
  Fact500);
  X.StrToExNum(
  "4.023872600770937735437024339230039857193748642107146E+2567",
  Fact1000);
  X.StrToExNum(
  "3.316275092450633241175393380576324038281117208105780E+5735",
  Fact2000);
  X.StrToExNum(
  "4.149359603437854085556867093086612170951119194931810E+9130",
  Fact3000);
END ExMathLib0.