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

MODULE ExNumbers;

IMPORT io, Cnv := Conversions, S := Strings;

CONST
  MaxExp * = 10000;
  MinExp * = -MaxExp;
  HighBoundsManArray * = 52; (* max possible digits--must be multiple of 4. *)

TYPE
  ExStatusType * = INTEGER;

CONST
  (* ExStatusType values *)
  Okay              *= 0;
  Overflow          *= 1;
  Underflow         *= 2;
  DivideByZero      *= 3;
  TooFewDigits      *= 4;
  TooManyDigits     *= 5;
  IllegalNumber     *= 6;
  UndefinedStorage  *= 7;
  IllegalOperator   *= 8;
  MismatchBraces    *= 9;

TYPE
  ExCompareType = INTEGER;

CONST
  (* ExCompareType values *)
  ExLess    *= 0;
  ExEqual   *= 1;
  ExGreater *= 2;

TYPE
  SignType = SHORTINT;

CONST
  (* SignType values *)
  positive *= 0;
  negative *= 1;

TYPE
  ManType   * = ARRAY (HighBoundsManArray DIV 4)+2 OF INTEGER;
  ExNumType * = RECORD
                  Man  -: ManType;
                  Sign -: SignType;
                  Zero -: BOOLEAN;
                  Exp  -: INTEGER;
                END;

VAR
  ExStatus * : ExStatusType;

  (* Useful constants *)
  e-, ln2-, ln10-, pi-, Ex0-, Ex1-: ExNumType;


CONST
  MaxLengthNumber = 2 * HighBoundsManArray;
  Dec = 10;

VAR
  MaxDigits, MaxQuads : INTEGER;


PROCEDURE SetMaxDigits *(D : INTEGER);
(* Set maximum digits in extended real numbers -- must be
   a multiple of 4 *)
BEGIN
  IF D < 4 THEN
    MaxDigits := 4;
    ExStatus := TooFewDigits;
  ELSIF D > HighBoundsManArray THEN
    MaxDigits := HighBoundsManArray;
    ExStatus := TooManyDigits;
  ELSE
    MaxDigits := D DIV 4;   (* Force a multiple of 4 *)
    IF D MOD 4 > 0 THEN INC(MaxDigits) END;
    MaxDigits := MaxDigits * 4;
  END;
  MaxQuads := MaxDigits DIV 4;
END SetMaxDigits;


PROCEDURE ExTimes10 *(VAR A : ExNumType);
(* A := A * 10 -- much faster than ExMult *)
BEGIN
  INC(A.Exp);
  IF A.Exp > MaxExp THEN
    ExStatus := Overflow;
  END;
END ExTimes10;


PROCEDURE ExDiv10 *(VAR A : ExNumType);
(* A := A / 10 -- much faster than ExDiv *)
BEGIN
  DEC(A.Exp);
  IF A.Exp < MinExp THEN
    ExStatus := Underflow;
  END;
END ExDiv10;


PROCEDURE IsZero *(A : ExNumType) : BOOLEAN;
VAR
  i : INTEGER;
  Zero : BOOLEAN;
BEGIN
  (* check for zero *)
  i := 0;
  Zero := TRUE;
  WHILE (i <= MaxQuads) AND Zero DO
    IF A.Man[i] # 0 THEN
      Zero := FALSE;
    END;
    INC(i);
  END;
  RETURN Zero;
END IsZero;


PROCEDURE ExShiftRight(Carry : INTEGER; VAR A : ExNumType);
(* shift all mantissa digits in A to the right one place.
   The most significant digit is replaced with the Carry. *)
VAR
  i : INTEGER;
BEGIN
  (* shift right *)
  FOR i := MaxQuads TO 1 BY -1 DO
    A.Man[i] := A.Man[i] DIV 10 + (A.Man[i-1] MOD 10) * 1000;
  END;

  (* put Carry in most significant position *)
  A.Man[0] := A.Man[0] DIV 10 + 1000 * Carry;
END ExShiftRight;


PROCEDURE ExShiftLeft(VAR A : ExNumType) : INTEGER;
(* shift all mantissa digits in A to the left one place.
   The digit shifted out of the number is returned.
   The least significant digit is replaced with zero. *)
VAR
  i, d : INTEGER;
BEGIN
  (* shift left *)
  d := A.Man[0] DIV 1000;
  FOR i := 0 TO MaxQuads DO
    A.Man[i] := (A.Man[i] MOD 1000) * 10 + A.Man[i+1] DIV 1000;
  END;

  (* put zero in least significant position *)
  A.Man[MaxQuads] := (A.Man[MaxQuads] MOD 1000) * 10;
  RETURN d;
END ExShiftLeft;


PROCEDURE ExChgSign *(VAR A : ExNumType);
(* A := -A *)
BEGIN
  IF A.Sign = positive THEN
    A.Sign := negative;
  ELSE
    A.Sign := positive;
  END;
END ExChgSign;


PROCEDURE ExAbs *(VAR A : ExNumType);
(* A := ABS(A) *)
BEGIN
  A.Sign := positive;
END ExAbs;


PROCEDURE ExNorm *(VAR A : ExNumType);
(* Normalise A *)
VAR d : INTEGER;
BEGIN
  (* normalise *)
  IF IsZero(A) THEN
    (* normalize zero *)
    A.Sign := positive;
    A.Exp := 0;
  ELSE
    (* shift mantissa to left until most significant digit is
       non-zero, increment exponent with each shift *)
    WHILE A.Man[0] DIV 1000 = 0 DO
      d := ExShiftLeft(A);
      ExDiv10(A);
    END;
  END;
END ExNorm;


PROCEDURE GetMaxDigits *() : INTEGER;
(* Get the current number of digits in extended real numbers *)
BEGIN
  RETURN MaxDigits;
END GetMaxDigits;


PROCEDURE GetExpMant *(x : ExNumType; VAR exp : INTEGER;
                       VAR mant : ExNumType);
(* Returned `mant' number will be between -10.0 and 10.0 *)
BEGIN
  exp := x.Exp;
  mant := x;
  mant.Exp := 0;
END GetExpMant;


PROCEDURE PutDigit(VAR A : INTEGER; Digit, Index : INTEGER);
BEGIN
  IF Index = 0 THEN
    A := A MOD 1000 + Digit * 1000;
  ELSIF Index = 1 THEN
    A := A DIV 1000 * 1000 + A MOD 100 + Digit * 100;
  ELSIF Index = 2 THEN
    A := A DIV 100 * 100 + A MOD 10 + Digit * 10;
  ELSE
    A := A DIV 10 * 10 + Digit;
  END;
END PutDigit;


PROCEDURE ExTrunc *(VAR A : ExNumType);
(* Truncate A so no decimal places are kept. *)
VAR
  i : INTEGER;
BEGIN
  IF A.Exp+1 < 0 THEN A := Ex0; RETURN END;
  FOR i := A.Exp+1 TO MaxDigits-1 DO
    (* zero these digits *)
    PutDigit(A.Man[i DIV 4], 0, i MOD 4);
  END;
END ExTrunc;


PROCEDURE ExFrac *(VAR A : ExNumType);
(* Keep only the fraction portion of A. *)
VAR
  i : INTEGER;
BEGIN
  FOR i := 0 TO A.Exp DO (* zero these digits *)
    PutDigit(A.Man[i DIV 4], 0, i MOD 4);
  END;
  ExNorm(A);             (* normalize the fraction *)
END ExFrac;


PROCEDURE ExToLongInt *(A : ExNumType) : LONGINT;
(* Convert the extended real number `A' into a INTEGER --
   saturating if necessary. *)
CONST
  MaxDigits = 10;
VAR
  Cnt : INTEGER;
  Int : LONGINT;
  Digit : INTEGER;
  Negative : BOOLEAN;
BEGIN
  Negative := FALSE;
  IF A.Sign = negative THEN
    Negative := TRUE;
    ExAbs(A);
  END;
  IF A.Exp < 0 THEN
    Int := 0;
  ELSIF A.Exp >= MaxDigits THEN
    Int := MAX(LONGINT);
  ELSE
    Int := 0;
    FOR Cnt := 0 TO A.Exp DO
      Digit := ExShiftLeft(A);
      IF Cnt = MaxDigits-1 THEN
        IF Int > MAX(LONGINT) DIV 10 THEN
          RETURN Int;
        END;
        IF (Int = MAX(LONGINT) DIV 10) & (Digit > 6) THEN
          Digit := 6;
        END;
      END;
      Int := Int * 10 + Digit;
    END;
  END;
  IF Negative THEN
    RETURN -Int;
  ELSE
    RETURN Int;
  END;
END ExToLongInt;


PROCEDURE ExCompare *(A, B : ExNumType) : ExCompareType;
(* Compares the two extended real numbers. *)
VAR
  Done : BOOLEAN;
  i : INTEGER;
BEGIN
  IF A.Sign # B.Sign THEN
    (* A and B have different signs *)
    IF A.Sign = positive THEN
      (* A and B have different signs and A is positive so A>B *)
      RETURN ExGreater;
    ELSE
      (* A and B have different signs and A is negative so A<B *)
      RETURN ExLess;
    END;
  ELSE
    (* A and B have the same sign *)
    IF (A.Exp # B.Exp) & NOT IsZero(B) & NOT IsZero(A) THEN
      IF A.Exp > B.Exp THEN
        (* A exponent > B exponent *)
        IF A.Sign = positive THEN
          RETURN ExGreater;
        ELSE
          RETURN ExLess;
        END;
      ELSE
        (* A exponent <= B exponent *)
        IF A.Sign = positive THEN
          RETURN ExLess;
        ELSE
          RETURN ExGreater;
        END;
      END;
    ELSE
      (* A & B have same sign and A exponent = B exponent *)
      Done := FALSE;
      i := 0;

      (* compare each digit until a difference is found or
         we reach the end *)
      WHILE (i <= MaxQuads) AND NOT Done DO
        IF A.Man[i] # B.Man[i] THEN
          Done := TRUE;
        ELSE
          INC(i);
        END;
      END;
      IF i > MaxQuads THEN
        (* end reached and all digits match *)
        RETURN ExEqual;
      ELSE
        (* compare different digits *)
        IF A.Man[i] < B.Man[i] THEN
          IF A.Sign = positive THEN
            RETURN ExLess;
          ELSE
            RETURN ExGreater;
          END;
        ELSE
          IF A.Sign = positive THEN
            RETURN ExGreater;
          ELSE
            RETURN ExLess;
          END;
        END;
      END;
    END;
  END;
END ExCompare;


PROCEDURE ExMin *(VAR A : ExNumType; B, C : ExNumType);
(* Return the smaller of B and C in A *)
BEGIN
  IF ExCompare(B, C) = ExGreater THEN
    A := C;
  ELSE
    A := B;
  END;
END ExMin;


PROCEDURE ExMax *(VAR A : ExNumType; B, C : ExNumType);
(* Return the larger of B and C in A *)
BEGIN
  IF ExCompare(B, C) = ExLess THEN
    A := C;
  ELSE
    A := B;
  END;
END ExMax;


PROCEDURE ExAddUtility(VAR A : ExNumType; B, C : ExNumType);
(* A := ABS(B) + ABS(C) *)
VAR
  i, j, joff, carry, quad, total : INTEGER;
  Exl1, Ex2 : ExNumType;
BEGIN
  IF IsZero(B) THEN
    A := C;
  ELSIF IsZero(C) THEN
    A := B;
  ELSE
    IF B.Exp > C.Exp THEN
      Exl1 := B;
      Ex2 := C;
    ELSE
      Exl1 := C;
      Ex2 := B;
    END;
    A := Ex0;
    A.Exp := Exl1.Exp;
    carry := 0;

    (* shift smallest number until quad-aligned relative to
       larger number *)
    j := (Exl1.Exp - Ex2.Exp) MOD 4;
    FOR i := j TO 1 BY -1 DO
      ExShiftRight(0, Ex2);
      INC(Ex2.Exp);
    END;
    joff := (Ex2.Exp - Exl1.Exp) DIV 4;

    (* add the two numbers together *)
    FOR i := MaxQuads TO 0 BY -1 DO
      (* j = index to Ex2 *)
      j := i + joff;

      (* check that j falls within array bounds *)
      IF (j >= 0) AND (j <= MaxQuads) THEN
        (* get quad digit from Ex2 *)
        quad := Ex2.Man[j];
      ELSE
        (* j is outside array bounds, use 0 for quad digit *)
        quad := 0;
      END;

      (* perform addition with carry *)
      total := Exl1.Man[i] + quad + carry;

      (* check for carry *)
      IF total >= 10000 THEN
        DEC(total, 10000);
        carry := 1;
      ELSE
        carry := 0;
      END;
      A.Man[i] := total;
    END;

    (* handle final carry *)
    IF carry = 1 THEN
      (* shift carry into top of mantissa *)
      ExShiftRight(carry, A);

      (* multiply by ten to update exponent *)
      ExTimes10(A);
    END;
  END;

  (* set ExStatus *)
  IF A.Exp > MaxExp THEN
    ExStatus := Overflow;
  END;
END ExAddUtility;


PROCEDURE ExSubUtility(VAR A : ExNumType; B, C : ExNumType);
(* A := ABS(B) - ABS(C) *)
VAR
  PositiveResult : BOOLEAN;
  i, j, joff, borrow, quad, result : INTEGER;
  Exl1, Ex2 : ExNumType;
BEGIN
  ExAbs(B);
  ExAbs(C);
  IF IsZero(B) THEN
    A := C;
  ELSIF IsZero(C) THEN
    A := B;
  ELSE
    IF B.Exp > C.Exp THEN
      Exl1 := B;
      Ex2 := C;
    ELSE
      Exl1 := C;
      Ex2 := B;
    END;
    PositiveResult := ExCompare(Exl1, Ex2) = ExGreater;
    A := Ex0;
    A.Exp := Exl1.Exp;
    borrow := 0;

    (* shift smallest number until quad-aligned relative to
       larger number *)
    j := (Exl1.Exp - Ex2.Exp) MOD 4;
    FOR i := j TO 1 BY -1 DO
      ExShiftRight(0, Ex2);
      INC(Ex2.Exp);
    END;
    joff := (Ex2.Exp - Exl1.Exp) DIV 4;

    (* subtract the two numbers *)
    FOR i := MaxQuads TO 0 BY -1 DO
      (* j = index to Ex2 *)
      j := i + joff;

      (* check that j falls within array bounds *)
      IF (j >= 0) AND (j <= MaxQuads) THEN
        (* get quad from Ex2 *)
        quad := Ex2.Man[j];
      ELSE
        (* j is outside array bounds, use 0 for quad *)
        quad := 0;
      END;

      (* perform subtraction with borrow *)
      IF PositiveResult THEN
        result := Exl1.Man[i] - quad - borrow;
      ELSE
        result := quad - Exl1.Man[i] - borrow;
      END;

      (* check for borrow *)
      IF result < 0 THEN
        INC(result, 10000);
        borrow := 1;
      ELSE
        borrow := 0;
      END;
      A.Man[i] := result;
    END;
  END;

  (* normalise *)
  ExNorm(A);

  (* adjust sign *)
  IF ExCompare(B, C) = ExLess THEN
    ExChgSign(A);
  END;
END ExSubUtility;


PROCEDURE ExAdd *(VAR A : ExNumType; B, C : ExNumType);
(* A = B + C *)
BEGIN
  IF B.Sign = C.Sign THEN
    (* B and C have the same sign -- just add *)
    ExAddUtility(A, B, C);
    IF B.Sign = negative THEN
      ExChgSign(A);
    END;
  ELSE
    (* B and C have different signs *)
    IF B.Sign = positive THEN
      ExSubUtility(A, B, C);
    ELSE
      ExSubUtility(A, C, B);
    END;
  END;
END ExAdd;


PROCEDURE ExSub *(VAR A : ExNumType; B, C : ExNumType);
(* A = B - C *)
BEGIN
  ExChgSign(C);   (* A = B + (-C) *)
  ExAdd(A, B, C);
END ExSub;


PROCEDURE ExRound *(VAR A : ExNumType; D : INTEGER);
(* A := Round(A) *)
VAR
  cindex, index, digit, i : INTEGER;
  Exl : ExNumType;
BEGIN
  IF D <= MaxDigits-1 THEN
    index := (D+1) DIV 4;
    digit := A.Man[index];
    cindex := (D + 1) MOD 4;
    IF cindex = 0 THEN
      digit := digit DIV 1000;
    ELSIF cindex = 1 THEN
      digit := digit DIV 100;
    ELSIF cindex = 2 THEN
      digit := digit DIV 10;
    END;
    IF digit MOD 10 >= 5 THEN
      (* round up *)
      Exl := Ex1;
      Exl.Exp := A.Exp - D;
      IF A.Sign = negative THEN
        ExChgSign(Exl);
      END;
      ExAdd(A, A, Exl);
    END;

    (* make remaining digits zero *)
    IF cindex = 0 THEN
      A.Man[index] := 0;
    ELSIF cindex = 1 THEN
      A.Man[index] := A.Man[index] DIV 1000 * 1000;
    ELSIF cindex = 2 THEN
      A.Man[index] := A.Man[index] DIV 100 * 100;
    ELSIF cindex = 3 THEN
      A.Man[index] := A.Man[index] DIV 10 * 10;
    END;
    FOR i := index+1 TO MaxQuads DO
      A.Man[i] := 0;
    END;
  END;
END ExRound;


PROCEDURE ExMult *(VAR A : ExNumType; B, C : ExNumType);
(* Return B * C *)
VAR
  i, j, carry : INTEGER;
  product : LONGINT;
  Exl : ExNumType;
BEGIN
  IF (ExCompare(B,Ex0) = ExEqual) OR (ExCompare(C,Ex0) = ExEqual) THEN
    (* multiplication by zero *)
    A := Ex0;
  ELSIF ExCompare(C,Ex1) = ExEqual THEN
    A := B;
  ELSIF ExCompare(B,Ex1) = ExEqual THEN
    A := C;
  ELSE
    (* real multiplication *)
    A := Ex0;
    FOR i := MaxQuads TO 0 BY -1 DO
      Exl := Ex0;
      Exl.Exp := B.Exp + C.Exp - i * 4 - 3;
      carry := 0;
      FOR j := MaxQuads TO 0 BY -1 DO
        product := LONG(B.Man[j]) * LONG(C.Man[i]) + LONG(carry);
        Exl.Man[j] := SHORT(product MOD 10000);
        carry := SHORT(product DIV 10000);
      END;

      (* check for final carry *)
      WHILE carry > 0 DO
        ExShiftRight(carry MOD 10, Exl);
        ExTimes10(Exl);
        carry := carry DIV 10;
      END;

      (* perform summation *)
      ExAddUtility(A, A, Exl);
    END;

    (* adjust product sign *)
    IF B.Sign # C.Sign THEN
      ExChgSign(A);
    END;
  END;
END ExMult;


PROCEDURE ExDiv *(VAR A : ExNumType; B, C : ExNumType);
(* A := B / C *)
VAR
  i, j : INTEGER;
  quotient : LONGINT;
  Exl1, Ex2 : ExNumType;
BEGIN
  IF IsZero(C) THEN
    (* attempt to divide by zero *)
    ExStatus := DivideByZero;
  ELSIF IsZero(B) THEN
    (* dividend = 0 *)
    A := Ex0;
  ELSIF ExCompare(C,Ex1) = ExEqual THEN
    (* divisor = 1 *)
    A := B;
  ELSE
    (* real division *)
    A := Ex0;
    A.Exp := B.Exp - C.Exp;

    (* adjust quotient sign *)
    IF B.Sign # C.Sign THEN
      ExChgSign(A);
    END;

    (* let Exl1 = ABS(B) / magnitude of B *)
    Exl1 := B;
    ExAbs(Exl1);
    Exl1.Exp := 0;

    (* let Ex2 = ABS(C) / magnitude of C *)
    Ex2 := C;
    ExAbs(Ex2);
    Ex2.Exp := 0;

    (* actual division *)
    FOR i := 0 TO MaxDigits-1 DO
      quotient := 0;
      WHILE ExCompare(Exl1, Ex2) >= ExEqual DO
        INC(quotient);
        ExSubUtility(Exl1, Exl1, Ex2);
      END;
      A.Man[i DIV 4] := A.Man[i DIV 4] * 10 + SHORT(quotient);
      ExDiv10(Ex2);
    END;

    (* normalize quotient *)
    ExNorm(A);
  END;
END ExDiv;


(* $CopyArrays- *)
PROCEDURE StrToExNum *(Str : ARRAY OF CHAR; VAR A : ExNumType);
(* Convert the string `Str' into an extended real number in A. *)
VAR
  Exp, NumbIndex, InCnt, EndCnt : INTEGER;
  ZeroFlag, NegativeExponent, LeftSide, InExponent : BOOLEAN;
  Done, NegExponent : BOOLEAN;
  ActiveChar : CHAR;

  PROCEDURE SetDigit(VAR Numb : INTEGER);
  BEGIN
    Numb := Numb * 10 + ORD(Str[InCnt]) - ORD('0');
  END SetDigit;

BEGIN
  (* initialize a few counters and stuff *)
  A := Ex0;
  InCnt := 0;             (* character counter *)
  Exp := 0;               (* working exponent *)
  LeftSide := TRUE;
  InExponent := FALSE;
  ZeroFlag := TRUE;
  NegativeExponent := FALSE;
  EndCnt := SHORT(S.Length(Str));
  NumbIndex := 0;

  (* set the sign of `A' to a negative -- if needed *)
  WHILE (InCnt < EndCnt) & (Str[InCnt] = ' ') DO INC(InCnt) END;
  IF Str[InCnt] = '-' THEN
    A.Sign := negative;
    INC(InCnt);
  END;
  WHILE InCnt < EndCnt DO
    ActiveChar := Str[InCnt];
    IF (ActiveChar >= '0') & (ActiveChar <= '9') THEN
      IF InExponent THEN
        SetDigit(Exp);
      ELSE
        IF NumbIndex < MaxDigits THEN  (* enter a digit *)
          SetDigit(A.Man[NumbIndex DIV 4]);
        END;
        IF ZeroFlag & (Str[InCnt] # '0') THEN
          ZeroFlag := FALSE;
        END;
        IF NOT ZeroFlag THEN
          INC(NumbIndex);
          IF LeftSide THEN INC(A.Exp) END;
        ELSE
          IF NOT LeftSide & (A.Exp <= 0) THEN DEC(A.Exp) END;
        END;
      END;
    ELSIF ActiveChar = '.' THEN
      IF ~LeftSide THEN ExStatus := IllegalNumber END;
      LeftSide := FALSE;
    ELSIF ActiveChar = 'E' THEN
      InExponent := TRUE;
      IF Str[InCnt+1] = '-' THEN
        NegativeExponent := TRUE;
        INC(InCnt);
      ELSIF Str[InCnt+1] = '+' THEN
        INC(InCnt);
      END;
    ELSIF ActiveChar = ' ' THEN
      (* do nothing if blanks are encountered *)
    ELSE
      ExStatus := IllegalNumber;
    END; (* IF *)
    INC(InCnt);
  END;

  (* fix up the last quad digits *)
  WHILE (NumbIndex DIV 4 <= MaxQuads) & (NumbIndex MOD 4 > 0) DO
    A.Man[NumbIndex DIV 4] := A.Man[NumbIndex DIV 4] * 10;
    INC(NumbIndex);
  END;

  (* Do some final fixes to the exponent *)
  IF NegativeExponent THEN
    DEC(A.Exp, Exp);
  ELSE
    INC(A.Exp, Exp);
  END;
  DEC(A.Exp);

  (* Ensure valid zero value *)
  IF IsZero(A) THEN A := Ex0 END;
END StrToExNum;


PROCEDURE GetDigit(VAR ExpStr : ARRAY OF CHAR; VAR StrCnt : INTEGER;
                   A : ExNumType; VAR ManIndex : INTEGER) : CHAR;
VAR Quad : LONGINT;
    Ok : BOOLEAN;
BEGIN
  (* Passing all parameters due to a bug in Oberon-2 when this
     was a local procedure *)
  INC(StrCnt);
  IF StrCnt = 4 THEN (* get a quad of digits *)
    Quad := A.Man[ManIndex];
    Ok := Cnv.IntToStr(Quad,ExpStr,Dec,5,'0');
    S.Delete(ExpStr, 0, 1);   (* remove leading digit *)
    INC(ManIndex);
    StrCnt := 0;
  END;
  RETURN ExpStr[StrCnt];
END GetDigit;


PROCEDURE ExNumToStr *(A : ExNumType; Decimal, ExpWidth : INTEGER;
                       VAR Str : ARRAY OF CHAR);
(* Convert the extended real number into a string `S'. *)
VAR
  pos, ManIndex, StrCnt, InCnt, Aexp, MaxExpWidth : INTEGER;
  ExpStr : ARRAY 41 OF CHAR;
  FixPoint, Ok : BOOLEAN;

  PROCEDURE ConcatChar(ch : CHAR);
  BEGIN
    Str[pos] := ch;
    INC(pos);
  END ConcatChar;

BEGIN
  (* initialize a few parameters *)
  pos := 0;
  StrCnt := 3;
  ManIndex := 0;
  ExpStr := '';

  (* force scientific notation for numbers too small or too large *)
  Aexp := ABS(A.Exp);
  MaxExpWidth := ExpWidth;
  IF ((ExpWidth = 0) AND (Aexp > MaxDigits)) OR (ExpWidth > 0) THEN
    (* force scientific notation *)
    IF Aexp > 9999 THEN ExpWidth := 5
    ELSIF Aexp > 999 THEN ExpWidth := 4
    ELSIF Aexp > 99 THEN ExpWidth := 3
    ELSIF Aexp > 9 THEN ExpWidth := 2
    ELSE ExpWidth := 1
    END;
  END;
  IF MaxExpWidth < ExpWidth THEN MaxExpWidth := ExpWidth END;

  (* add the negative sign to the number *)
  IF A.Sign = negative THEN ConcatChar('-') END;

  (* ensure we don't exceed the maximum digits *)
  FixPoint := Decimal # 0;
  IF (Decimal > MaxDigits) OR NOT FixPoint THEN
    Decimal := MaxDigits-1;
  END;

  (* convert the number into scientific notation *)
  IF MaxExpWidth > 0 THEN
    ExRound(A, Decimal);    (* round to appropriate decimal places *)
    ConcatChar(GetDigit(ExpStr, StrCnt, A, ManIndex)); (* leading digit *)
    ConcatChar('.');        (* decimal point *)
    FOR InCnt := 1 TO Decimal DO
      ConcatChar(GetDigit(ExpStr, StrCnt, A, ManIndex));     (* add following digits *)
    END;

    (* add the exponent *)
    ConcatChar('E');
    IF A.Exp >= 0 THEN ConcatChar('+') ELSE ConcatChar('-') END;
    ConcatChar(0X);                        (* terminate the string *)

    Ok := Cnv.IntToStr(Aexp,ExpStr,Dec,SHORT(MaxExpWidth),'0');
    S.Append(Str, ExpStr);
  ELSE
    (* format a non-scientific number *)
    ExRound(A, Decimal+A.Exp); (* round to decimal places *)
    IF A.Exp < 0 THEN
      ConcatChar('0');         (* leading digit *)
      ConcatChar('.');         (* decimal point *)
      FOR InCnt := 2 TO ABS(A.Exp) DO   (* pad with leading zeros *)
        ConcatChar('0');
      END;
      INC(Decimal, A.Exp+1);
    END;
    InCnt := 0;
    REPEAT
      ConcatChar(GetDigit(ExpStr, StrCnt, A, ManIndex));
      IF InCnt > A.Exp THEN
        DEC(Decimal);
      ELSIF InCnt = A.Exp THEN
        ConcatChar('.');
      END;
      INC(InCnt);
    UNTIL (InCnt = MaxDigits) OR (Decimal = 0);
    ConcatChar(0X);

    (* remove any trailing zeros and unneeded digits *)
    InCnt := pos - 2;
    WHILE (InCnt > 1) & (Str[InCnt] = '0') & NOT FixPoint DO
      Str[InCnt] := 0X;
      DEC(InCnt);
    END;
  END;
END ExNumToStr;


PROCEDURE WriteExNum *(A : ExNumType;
                       Width, Decimal, ExpWidth : INTEGER);
(* Write out A to the current output stream in a field of
   `Width' characters, with `Decimal' decimal places, and
   `ExpWidth' exponent width. *)
VAR
  Str : ARRAY MaxLengthNumber+1 OF CHAR;
  i, len : INTEGER;
BEGIN
  ExNumToStr(A, Decimal, ExpWidth, Str);
  len := SHORT(S.Length(Str));
  IF Width >= len THEN
    FOR i := 1 TO Width-len DO io.Write(" ") END;
  END;
  io.WriteString(Str);
END WriteExNum;


PROCEDURE ExNumb *(LeftMan : LONGINT; RightMan : LONGINT;
                   ExpShift : INTEGER; VAR A : ExNumType);
(* create an extended real number which has LeftMan to the left
   of the decimal point and RightMan to the right. The ExpShift
   quantity can shift the decimal point to the right for negative
   values; to the left for positive values. *)
VAR
  i : INTEGER;
BEGIN
  A := Ex0;
  IF LeftMan < 0 THEN
    A.Sign := negative;
    LeftMan := -LeftMan;
  END;
  WHILE RightMan # 0 DO
    ExShiftRight(SHORT(RightMan MOD 10), A);(* shift right 1 position *)
    RightMan := RightMan DIV 10;
  END;
  WHILE LeftMan # 0 DO
    ExShiftRight(SHORT(LeftMan MOD 10), A); (* shift right 1 position *)
    ExTimes10(A);                    (* adjust the exponent *)
    LeftMan := LeftMan DIV 10;
  END;
  ExDiv10(A);                        (* final exponent adjust *)
  INC(A.Exp, ExpShift);              (* shift the decimal point *)
  IF A.Exp > MaxExp THEN             (* signal any errors *)
    ExStatus := Overflow;
  ELSIF A.Exp < MinExp THEN
    ExStatus := Underflow;
  END;
END ExNumb;


BEGIN
  (* create extended number 0 *)
  Ex0.Sign := positive;
  FOR MaxDigits := 0 TO LEN(Ex0.Man)-1 DO
    Ex0.Man[MaxDigits] := 0;
  END;
  Ex0.Exp := 0;

  (* default to max number of digits *)
  SetMaxDigits(HighBoundsManArray);

  (* create some extended number constants *)
  ExNumb(1, 0, 0, Ex1);     (* 1.0 *)

  StrToExNum(
  "3.14159265358979323846264338327950288419716939937511", pi);
  StrToExNum(
  "2.71828182845904523536028747135266249775724709369996", e);
  StrToExNum(
  "0.69314718055994530941723212145817656807550013436026", ln2);
  StrToExNum(
  "2.30258509299404568401799145468436420760110148862877", ln10);
END ExNumbers.