-- Copyright (C) 1986  by Xerox Corporation. All rights reserved. 
-- LibmSupportImpl.mesa
-- NFS    22-Apr-86 16:27:27

-- Support functions for libm floating pt. functions, and library functions
-- ldexp(), modf(), fabs(), fmod(), and abs().

DIRECTORY
  DoubleIeee,
  DoubleReal,
  DoubleRealFns,
  DoubleRealOps,
  Inline,
  Libm,
  LibmSupport;
LibmSupportImpl: PROGRAM
  IMPORTS DoubleIeee, DoubleReal, DoubleRealFns, DoubleRealOps, Inline
  EXPORTS Libm, LibmSupport =
  {
  OPEN DoubleReal, DoubleIeee;

  copysign: PUBLIC PROCEDURE [x, y: Double] RETURNS [Double] = {
    a: VeryLongNumber ← LOOPHOLE[x];
    b: VeryLongNumber = LOOPHOLE[y];
    IF DblBitOn[b.hi, DblHiBit] THEN  -- y is negative
      a.hi ← Inline.DBITOR[a.hi, DblHiBit]
    ELSE  -- y is positive
      a.hi ← Inline.DBITAND[a.hi, DblNotHiBit];
    RETURN[LOOPHOLE[a]];
    };

  scalb: PUBLIC PROCEDURE [x: Double, N: INTEGER] RETURNS [Double] = {
    a: DblExt ← UnpackDbl[x];
    a.exp ← a.exp + N;
    IF a.exp >= DblNaNExponent THEN a.det.type ← infinity;
    RETURN[PackDbl[@a]];
    };

  ldexp: PUBLIC PROCEDURE [x: Double, exp: INTEGER] RETURNS [Double] = {
    RETURN[scalb[x, exp]]; };

  logb: PUBLIC PROCEDURE [x: Double] RETURNS [Double] = {
    a: DblExt = UnpackDbl[x];
    SELECT a.det.type FROM
      zero => RETURN[MinusInfinity];
      infinity => RETURN[PlusInfinity];
      nan => RETURN[x];
      ENDCASE => RETURN[DFloat[LONG[a.exp]]];
    };

  finite: PUBLIC PROCEDURE [x: Double] RETURNS [INTEGER] = {
    true: INTEGER = 1;
    false: INTEGER = 0;
    RETURN[
      IF Inline.DBITAND[LOOPHOLE[x, VeryLongNumber].hi, DblExponentMask] =
      DblExponentMask THEN false ELSE true];
    };

  drem: PUBLIC PROCEDURE [x, p: Double] RETURNS [Double] = {
    RETURN[DoubleReal.DFRem[x, p]]; };

  modf: PUBLIC PROCEDURE [d: Double, iptr: LONG POINTER TO Double]
    RETURNS [Double] = {
    fixMode: DoubleRealOps.Mode = [round: rz, traps: DoubleReal.NoExceptions];
    iptr↑ ← DoubleRealOps.DRoundF[d, fixMode];
    RETURN[DoubleReal.DFSub[d,iptr↑]];
  };

  fabs: PUBLIC PROCEDURE [d: Double] RETURNS [Double] = {
    LOOPHOLE[d, VeryLongNumber].hi ← Inline.DBITAND[
      LOOPHOLE[d, VeryLongNumber].hi, DblNotHiBit];
    RETURN[d];
    };

  fmod: PUBLIC PROCEDURE [x, y: Double] RETURNS [Double] = {
    IF y = PlusZero OR y = MinusZero THEN RETURN[x]
      <<Definition of fmod (when y # 0) is not clear.
     The definition used here makes the Rapa validation test pass.>>
    ELSE {
      py: Double = fabs[y];
      rem: Double ← DFRem[fabs[x], py];
      -- If rem < 0 then add py.
      IF DblBitOn[LOOPHOLE[rem, VeryLongNumber].hi, DblHiBit] AND rem # MinusZero
        THEN rem ← DFAdd[rem, py];
      RETURN[copysign[rem, x]];
      };
    };

  abs: PUBLIC PROCEDURE [i: INTEGER] RETURNS [INTEGER] = {RETURN[ABS[i]]; };

  sqrt: PUBLIC PROCEDURE [x: Double] RETURNS [Double] = {
    RETURN[DoubleRealFns.SqRt[x]]; };

  }.