(* Copyright (C) 1992, Digital Equipment Corporation                         *)
(* All rights reserved.                                                      *)
(* See the file COPYRIGHT for a full description.                            *)
(*                                                                           *)
(* Last modified on Tue Jun 16 16:46:27 PDT 1992 by muller                   *)

(* Created by stolfi in Sep 1988                               *)

UNSAFE MODULE RealExtra;

(***********************)
(* IEEE (MIPS) VERSION *)
(***********************)

IMPORT Word;

EXCEPTION InvalidArgument;

CONST 
  LastW    = 16_7f7fffff;           (* Largest finite number *)
  FirstW   = Word.Not(16_00800000); (* = Smallest finite number = 16_ff7fffff *)

  MinPosW  = 16_00000001;           (* Smallest positive number *)
  MaxNegW  = Word.Not(16_7ffffffe); (* Largest negative number = 16_80000001 *)

  MinNormW = 16_00800000;           (* Smallest positive normalized number *)

  NegZeroW = Word.Not(16_7fffffff); (* Minus zero *)

  PInftyW  = 16_7f800000;           (* Plus Infinity *)
  NInftyW  = Word.Not(16_007fffff); (* Minus Infinity = 16_ff800000 *)

CONST
  FrMask = 16_007fffff;            (* Manifest significand (fraction) bits *)
  FrLeng = 23;                     (* Number of manifest significand bits *)
  HidBit = 16_00800000;            (* Hidden units bit of significand (actually lsb of exponent) *)
  ExMask = 16_7f800000;            (* Exponent bits *)
  SFMask = Word.Not(16_7f800000);  (* = 16_807fffff; Sign and fraction bits *)
  HalfEx = 16_3f000000;            (* Exponent of numbers in [0.5 _ 1.0) *)
  SgnBit = Word.Not(16_7fffffff);  (* = 16_80000000; Sign bit *)
  EFMask = 16_7fffffff;            (* Exponent and fraction bits *)
  ExBias = 126;
  MaxExp = 254;                    (* Maximum exponent field of a finite number *)

PROCEDURE PRED(READONLY x: REAL): REAL RAISES {Overflow} =
  VAR w: INTEGER := LOOPHOLE(x, Word.T);
  BEGIN 
    (* This code treats -0.0 as equivalent to 0.0 *)
    (* Uncomment the (*!...!*) to accept and return IEEE infinities *)
    IF Word.And(w, ExMask) = ExMask THEN
      (*!
        IF w = NInftyW THEN
          RAISE Overflow
        ELSIF w = PInftyW THEN
          RETURN Last
        ELSE
          (* Must be NaN *)
          RAISE InvalidArgument
        END
      !*)
      RAISE InvalidArgument
    ELSIF w > 0 THEN
      RETURN LOOPHOLE(w - 1, REAL)
    ELSIF w = 0 OR w = NegZeroW THEN
      RETURN LOOPHOLE(MaxNegW, REAL)
    ELSE (* x < 0 *)
      IF w = FirstW THEN
        (*! RETURN LOOPHOLE(NInftyW, REAL) !*)
        RAISE Overflow
      ELSE
        RETURN LOOPHOLE(w + 1, REAL)
      END
    END;
  END PRED;
  
PROCEDURE SUCC(READONLY x: REAL): REAL RAISES {Overflow} =
  VAR w: INTEGER := LOOPHOLE(x, Word.T);
  BEGIN 
    (* This code treats -0.0 as equivalent to 0.0 *)
    (* Uncomment the (*!...!*) to accept and return IEEE infinities *)
    IF Word.And(w, ExMask) = ExMask THEN
      (*!
        IF w = NInftyW THEN
          RETURN First
        ELSIF w = PInftyW THEN
          RAISE Overflow
        ELSE
          (* Must be NaN *)
          RAISE InvalidArgument
        END
      !*)
      RAISE InvalidArgument
    ELSIF w > 0 THEN
      IF w = LastW THEN
        (*! RETURN LOOPHOLE(PInftyW, REAL) !*)
        RAISE Overflow
      ELSE
        RETURN LOOPHOLE(w + 1, REAL)
      END
    ELSIF w = 0 OR w = NegZeroW THEN
      RETURN LOOPHOLE(MinPosW, REAL)
    ELSE (* x < 0 *)
      RETURN LOOPHOLE(w - 1, REAL)
    END;
  END SUCC;

PROCEDURE SIGN(READONLY x: REAL): [-1 .. +1] =
  VAR w: INTEGER := LOOPHOLE(x, Word.T);
  BEGIN
    (* Uncomment the (*!...!*) to accept IEEE infinities *)
    IF Word.And(w, ExMask) = ExMask THEN
      (*!
        IF w = PInftyW THEN
          RETURN +1
        ELSIF w = NInftyW THEN
          RETURN -1
        ELSE
          (* Must be NaN *)
          RAISE InvalidArgument
        END
      !*)
      RAISE InvalidArgument
    ELSIF w > 0 THEN 
      RETURN +1
    ELSIF w = 0 OR w = NegZeroW THEN
      RETURN 0
    ELSE
      RETURN -1
    END
  END SIGN;

PROCEDURE FRACTION(READONLY x: (*FINITE*) REAL): (*FINITE*) REAL RAISES {} =
  VAR w: INTEGER := LOOPHOLE(x, Word.T);
      abs: INTEGER  := Word.And(w, EFMask);
  BEGIN 
    IF abs = 0 THEN 
      RETURN 0.0
    ELSIF Word.And(w, ExMask) = ExMask THEN
      (* Must be Infinity or NaN *)
      RAISE InvalidArgument
    ELSIF abs >= MinNormW THEN
      RETURN LOOPHOLE(Word.Or(Word.And(w, SFMask), HalfEx), REAL)
    ELSIF abs = MinPosW THEN
      IF w < 0 THEN RETURN -0.5 ELSE RETURN 0.5 END
    ELSE
      (* Needs to normalize by hand: *)
      VAR f := Word.And(w, FrMask);
      BEGIN
        IF Word.And(f, 16_00ffff00) = 0 THEN f := Word.Shift(f, 16) END;
        IF Word.And(f, 16_00ff0000) = 0 THEN f := Word.Shift(f, 08) END;
        IF Word.And(f, 16_00f00000) = 0 THEN f := Word.Shift(f, 04) END;
        IF Word.And(f, 16_00c00000) = 0 THEN f := Word.Shift(f, 02) END;
        IF Word.And(f, 16_00800000) = 0 THEN f := Word.Shift(f, 01) END;
        (* Discard hidden bit: *)
        f := Word.And(f, FrMask);
        (* Put back minus sign: *)
        IF w < 0 THEN f := Word.Or(f, SgnBit) END;
        RETURN LOOPHOLE(Word.Or(f, HalfEx), REAL);
      END
    END
  END FRACTION;
  
PROCEDURE EXPONENT(READONLY x: (*FINITE*) REAL): INTEGER RAISES {} =
  VAR w: INTEGER := LOOPHOLE(x, Word.T);
      abs: INTEGER  := Word.And(w, EFMask);
  BEGIN 
    IF abs = 0 THEN 
      RETURN 0
    ELSIF Word.And(w, ExMask) = ExMask THEN
      (* Must be Infinity or NaN *)
      RAISE InvalidArgument
    ELSIF abs >= MinNormW THEN
      RETURN Word.Shift(Word.And(w, ExMask), -FrLeng) - ExBias
    ELSIF abs = MinPosW THEN
      RETURN 1 - ExBias - FrLeng
    ELSE
      (* Needs to normalize by hand: *)
      VAR f := Word.And(w, FrMask);
          e := 1 - ExBias;
      BEGIN
        IF Word.And(f, 16_00ffff00) = 0 THEN f := Word.Shift(f, 16); e := e-16 END;
        IF Word.And(f, 16_00ff0000) = 0 THEN f := Word.Shift(f, 08); e := e-08 END;
        IF Word.And(f, 16_00f00000) = 0 THEN f := Word.Shift(f, 04); e := e-04 END;
        IF Word.And(f, 16_00c00000) = 0 THEN f := Word.Shift(f, 02); e := e-02 END;
        IF Word.And(f, 16_00800000) = 0 THEN f := Word.Shift(f, 01); e := e-01 END;
        RETURN e
      END
    END
  END EXPONENT;

PROCEDURE NORMALIZE(VAR x: (*FINITE*) REAL; VAR (*OUT*) e: INTEGER) RAISES {} =
  VAR w: INTEGER := LOOPHOLE(x, Word.T);
      abs: INTEGER  := Word.And(w, EFMask);
  BEGIN 
    IF abs = 0 THEN 
      x := 0.0;
      e := 0
    ELSIF Word.And(w, ExMask) = ExMask THEN
      (* Must be Infinity or NaN *)
      RAISE InvalidArgument
    ELSIF abs >= MinNormW THEN
      e := Word.Shift(Word.And(w, ExMask), -FrLeng) - ExBias;
      x := LOOPHOLE(Word.Or(Word.And(w, SFMask), HalfEx), REAL)
    ELSIF abs = MinPosW THEN
      e := 1 - ExBias - FrLeng;
      IF w < 0 THEN x := -0.5 ELSE x := 0.5 END;
    ELSE
      (* Needs to normalize by hand: *)
      VAR f := Word.And(w, FrMask);
      BEGIN
        e := 1 - ExBias;
        IF Word.And(f, 16_00ffff00) = 0 THEN f := Word.Shift(f, 16); e := e-16 END;
        IF Word.And(f, 16_00ff0000) = 0 THEN f := Word.Shift(f, 08); e := e-08 END;
        IF Word.And(f, 16_00f00000) = 0 THEN f := Word.Shift(f, 04); e := e-04 END;
        IF Word.And(f, 16_00c00000) = 0 THEN f := Word.Shift(f, 02); e := e-02 END;
        IF Word.And(f, 16_00800000) = 0 THEN f := Word.Shift(f, 01); e := e-01 END;
        (* Discard hidden bit: *)
        f := Word.And(f, FrMask);
        (* Put back minus sign: *)
        IF w < 0 THEN f := Word.Or(f, SgnBit) END;
        x := LOOPHOLE(Word.Or(f, HalfEx), REAL);
      END
    END
  END NORMALIZE;

CONST
  MaxScale = 256 + FrLeng;  (* Scaling a nonzero number by 2^MaxScale is a sure overflow *)
  
PROCEDURE SCALE(VAR x: (*FINITE*) REAL; exp: INTEGER) RAISES {(*! Overflow !*)} =
  VAR w: INTEGER := LOOPHOLE(x, Word.T);
      abs: INTEGER := Word.And(w, EFMask);
  BEGIN
    (* Uncomment the (*!...!*) to raise Overflow instead of clipping to [First _ Last] *)
    IF abs = 0 THEN
      x := 0.0;
      RETURN
    ELSIF exp = 0 THEN
      RETURN
    ELSIF exp > MaxScale THEN
      (*! RAISE Overflow !*)
      IF w > 0 THEN x := Last ELSE x := First END;
      RETURN
    ELSIF exp < -MaxScale THEN
      x := 0.0;
      RETURN
    ELSE
      VAR ebold: INTEGER := Word.Shift(Word.And(w, ExMask), -FrLeng); (* Old exp + bias *)
          ebnew: INTEGER := ebold + exp;                                  (* New exp + bias *)
      BEGIN
        IF ebold > MaxExp THEN
          (* Must be infinity or NaN *)
          RAISE InvalidArgument
        ELSIF ebnew < -FrLeng THEN
          x := 0.0
        ELSIF ebold > 0 THEN 
          (* Input was normalized *)
          IF ebnew > 0 THEN
            IF ebnew > MaxExp THEN
              (*! RAISE Overflow !*)
              IF w > 0 THEN x := Last ELSE x := First END
            ELSE
              x := LOOPHOLE(Word.Or(Word.And(w, SFMask), Word.Shift(ebnew, FrLeng)), REAL)
            END
          ELSE
            (* Drats! input was normalized, but result is not. *)
            WITH f = Word.Shift(Word.Or(Word.And(w, FrMask), HidBit), ebnew-1) DO
              IF w < 0 THEN 
                x := LOOPHOLE(Word.Or(f, SgnBit), REAL)
              ELSE
                x := LOOPHOLE(f, REAL)
              END
            END;
          END
        ELSE
          (* Drats! Input was denormalized. *)
          VAR f := Word.And(w, FrMask);
          BEGIN
            IF ebnew < 0 THEN
              (* Need only shift fraction to the right *)
              f := Word.Shift(f, exp);
              IF f = 0 THEN
                x := 0.0
              ELSE
                IF w < 0 THEN f := Word.Or(f, SgnBit) END;
                x := LOOPHOLE(f, REAL)
              END
            ELSE
              <* ASSERT ebnew > 0 *>
              <* ASSERT ebnew = exp *>
              (* Need to shift fraction to the left, taking care not to overshoot *)
              IF f = 1 THEN
                (* Special handling for the common case of MinPosW, MaxNegW: *)
                IF ebnew < FrLeng THEN
                  f := Word.Shift(f, ebnew);
                ELSE
                  ebnew := ebnew - FrLeng + 1;
                  IF ebnew > MaxExp THEN
                    (*! RAISE Overflow !*)
                    IF w > 0 THEN x := Last ELSE x := First END;
                    RETURN
                  ELSE
                    f := Word.Shift(ebnew, FrLeng)
                  END
                END;
                IF w < 0 THEN f := Word.Or(f, SgnBit) END;
                x := LOOPHOLE(f, REAL);
              ELSE
                (* Other denormalized numbers: *)
                IF ebnew >= 16 AND Word.And(f, 16_00ffff00) = 0 THEN f := Word.Shift(f, 16); DEC(ebnew, 16) END;
                IF ebnew >= 8  AND Word.And(f, 16_00ff0000) = 0 THEN f := Word.Shift(f, 08); DEC(ebnew, 08) END;
                IF ebnew >= 4  AND Word.And(f, 16_00f00000) = 0 THEN f := Word.Shift(f, 04); DEC(ebnew, 04) END;
                IF ebnew >= 2  AND Word.And(f, 16_00c00000) = 0 THEN f := Word.Shift(f, 02); DEC(ebnew, 02) END;
                IF ebnew >= 1  AND Word.And(f, 16_00800000) = 0 THEN f := Word.Shift(f, 01); DEC(ebnew, 01) END;
                IF Word.And(f, HidBit) = 0 THEN
                  (* Result is still denormalized, exponent must be zero: *)
                  <* ASSERT ebnew = 0 *>
                ELSIF ebnew >= MaxExp THEN
                  (*! RAISE Overflow !*)
                  IF w > 0 THEN x := Last ELSE x := First END;
                  RETURN
                ELSE
                  (* Result is normalized, get rid of hidden bit and correct exponent: *)
                  INC(ebnew);
                  f := Word.Or(Word.And(f, FrMask), Word.Shift(ebnew, FrLeng))
                END;
                IF w < 0 THEN f := Word.Or(f, SgnBit) END;
                x := LOOPHOLE(f, REAL);
              END;
            END
          END
        END
      END
    END
  END SCALE;
  
BEGIN
  First   := LOOPHOLE(FirstW, REAL);
  Last    := LOOPHOLE(LastW, REAL);
  MinPos  := LOOPHOLE(MinPosW, REAL);
  MinNorm := LOOPHOLE(MinNormW, REAL);
END RealExtra.