-- RealFnsImpl.mesa
-- Last modified: Stewart on August 27, 1982 1:07 pm

DIRECTORY
  Ieee USING [CVExtended, ExponentBias, SingleReal, Unpack],
  Real USING [Extended, Fix, RealError, RealException],
  Inline USING [LowHalf],
  RealFns;

RealFnsImpl: CEDAR PROGRAM IMPORTS Ieee, Inline, Real EXPORTS RealFns =
  BEGIN OPEN Ieee;

  Precision: CARDINAL _ 5;
  LoopCount: INTEGER _ 3;

  Ln2: REAL = 0.693147181;
  LogBase2ofE: REAL = 1/Ln2;

  --compute e^n using continued fractions
  --See ACM Algorithm 229, Morelock, James C., "Elementary functions
  --	by continued fractions"
  Exp: PUBLIC PROC [x: REAL] RETURNS [REAL] = {
    r, f, kludge: REAL;
    k: LONG INTEGER;
    i, ikludge: INTEGER;
    IF x = 0 THEN RETURN[1];

    x _ x*LogBase2ofE; --exponent of 2 (instead of e)
    k _ Real.Fix[x]; --integer part
    x _ x - k; --fractional part

    x _ x*Ln2; --now result = 2^k * exp(x) and 0<=x<Ln2
    r _ x*x;

    --f_4*LoopCount+2;
    i _ 4*LoopCount + 2;
    f _ i;

    FOR i DECREASING IN [1..LoopCount] DO
      ikludge _ 4*i - 2; kludge _ ikludge; f _ kludge + r/f; ENDLOOP;

    --unscaled result is (f+x)/(f-x)
    x _ (f + x)/(f - x); --scale: multipy by 2^k by adding k to the exponent
    x _ Pow2[x, LOOPHOLE[Inline.LowHalf[LOOPHOLE[k]]]];
    RETURN[x];
    };

  Pow2: PROC [x: REAL, exp: INTEGER] RETURNS [REAL] = INLINE {
    fl: SingleReal _ LOOPHOLE[1.0, SingleReal];
    IF exp NOT IN [-126..127] THEN ERROR;
    fl.exp _ exp + ExponentBias;
    --scale: multipy by 2^exp by adding exp to the exponent
    x _ x*LOOPHOLE[fl, REAL];
    RETURN[x];
    };

  Log: PUBLIC PROC [base, arg: REAL] RETURNS [REAL] = {
    --logxy = Ln(y)/Ln(x)
    x, y: REAL;
    x _ Ln[base];
    y _ Ln[arg];
    RETURN[y/x];
    };

  Ln: PUBLIC PROC [x: REAL] RETURNS [REAL] = {
    fl: SingleReal;
    i: CARDINAL;
    m: INTEGER;
    kludge: REAL;
    OneMinusX, OneMinusXPowerRep, PreviousVal: REAL;
    IF x < 0 THEN {
      [] _ SIGNAL Real.RealException[flags: [invalidOperation: TRUE], vp: NEW[Real.Extended _ Ieee.CVExtended[Ieee.Unpack[x]]]];
      ERROR Real.RealError[];
      };

    --if x = z * 2^m then Ln(x) = m * Ln2 + Ln(z).  Here, 1/2 <= z < 1
    fl _ LOOPHOLE[x, SingleReal];
    m _ fl.exp - (ExponentBias - 1);
    fl.exp _ ExponentBias - 1; --now 1/2 <= z < 1
    x _ LOOPHOLE[fl, REAL];

    --special case for limit condition (x a power of 2)
    IF x = 0.5 THEN {m _ m - 1; kludge _ m; x _ Ln2*kludge; RETURN[x]; };

    OneMinusX _ 1 - x;
    x _ m*Ln2;
    OneMinusXPowerRep _ 1;
    FOR i IN [1..100] DO
      PreviousVal _ x;
      OneMinusXPowerRep _ OneMinusXPowerRep*OneMinusX;
      x _ x - OneMinusXPowerRep/i;
      IF x = PreviousVal THEN EXIT;
      ENDLOOP;

    RETURN[x];
    };


  SqRt: PUBLIC PROC [x: REAL] RETURNS [REAL] = {
    epsilon: REAL = Pow2[1, -23];
    y, x1: REAL;
    fl: SingleReal;


    IF x <= 0 THEN RETURN[0];

    fl _ LOOPHOLE[x, SingleReal];
    fl.exp _ (fl.exp - ExponentBias)/2 + ExponentBias;
    y _ LOOPHOLE[fl, REAL]; --initial guess = half original exponent
    DO
      x1 _ y*y;
      y _ y - (x1 - x)/(2*y); --slope of the curve here is 1/2
      IF ABS[(x1 - x)/x] <= epsilon THEN EXIT; -- may need ABS[]

      ENDLOOP;

    RETURN[y];
    };

  Root: PUBLIC PROC [index, arg: REAL] RETURNS [REAL] = {
    --yth root of x = e^(Ln(x)/y)
    RETURN[Exp[Ln[arg]/index]];
    };


  Power: PUBLIC PROC [base, exponent: REAL] RETURNS [REAL] = {
    --x^y = e^(y*Ln(x))
    x: REAL;
    IF base = 0 THEN IF exponent = 0 THEN RETURN[1] ELSE RETURN[0];
    x _ Ln[base];
    RETURN[Exp[exponent*x]];
    };

  Cos: PUBLIC PROC [radians: REAL] RETURNS [cos: REAL] = {
    -- This function is good to 7.33 decimal places and is taken from:
    -- Computer Approximations, John F. Hart et.al. pp118(top 2nd col).
    negsign: REAL _ 1;
    nangle: REAL _ ABS[radians - Real.Fix[radians*recpi]*twoPI];
    SELECT nangle FROM
      IN [PIovr2..PI3ovr2) => BEGIN negsign _ -1; nangle _ ABS[PI - nangle]; END;
      IN [PI3ovr2..twoPI) => nangle _ twoPI - nangle;
      ENDCASE;
    nangle _ nangle*nangle;
    cos _ negsign*(cp0 + nangle*(cp1 + nangle*(cp2 + nangle*(cp3 + nangle*cp4))));
    };

  CosDeg: PUBLIC PROC [degrees: REAL] RETURNS [cos: REAL] = {
    radians: REAL _ degrees*degtorad; cos _ Cos[radians]; };

  Sin: PUBLIC PROC [radians: REAL] RETURNS [sin: REAL] = {
    radians _ PIovr2 - radians; sin _ Cos[radians]; };

  SinDeg: PUBLIC PROC [degrees: REAL] RETURNS [sin: REAL] = {
    radians: REAL _ degrees*degtorad; sin _ Sin[radians]; };

  Tan: PUBLIC PROC [radians: REAL] RETURNS [tan: REAL] = {
    tan _ Sin[radians]/Cos[radians]; };

  TanDeg: PUBLIC PROC [degrees: REAL] RETURNS [tan: REAL] = {
    radians: REAL _ degrees*degtorad; tan _ Sin[radians]/Cos[radians]; };

  ArcTan: PUBLIC PROC [y, x: REAL] RETURNS [radians: REAL] = {
    -- This function is good to 8.7 decimal places and is taken from:
    -- Computer Approximations, John F. Hart et.al. pp129(top 2nd col).
    s, v, t, t2, c, q: REAL;
    IF ABS[x] <= ABS[smallnumber] THEN
      IF ABS[y] <= ABS[smallnumber] THEN RETURN[0]
      ELSE IF y < 0 THEN RETURN[-PIovr2] ELSE RETURN[PIovr2];
    IF x < 0 THEN IF y < 0 THEN {q _ -PI; s _ 1; } ELSE {q _ PI; s _ -1; }
    ELSE IF y < 0 THEN {q _ 0; s _ -1; } ELSE {q _ 0; s _ 1; };
    v _ ABS[y/x];
    SELECT v FROM
      IN [0..tanPI16) => {t _ v; c _ 0; };
      IN [tanPI16..tan3PI16) => {t _ (x2) - ((x22)/(x2 + v)); c _ PIovr8; };
      IN [tan3PI16..tan5PI16) => {t _ (x4) - ((x44)/(x4 + v)); c _ PIovr4; };
      IN [tan5PI16..tan7PI16) => {t _ (x6) - ((x66)/(x6 + v)); c _ PI3ovr8; };
      >= tan7PI16 => {t _ -1/v; c _ PIovr2; };
      ENDCASE;
    t2 _ t*t;
    radians _ s*(t*(p0 + t2*(p1 + t2*(p2 + t2*p3))) + c) + q;
    };

  ArcTanDeg: PUBLIC PROC [y, x: REAL] RETURNS [degrees: REAL] = {
    radians: REAL _ ArcTan[y, x]; degrees _ radians*radtodeg; };

  AlmostZero: PUBLIC PROC [x: REAL, distance: INTEGER [-126..127]]
    RETURNS [BOOLEAN] = {
    fl: Ieee.SingleReal _ LOOPHOLE[x];
    RETURN[fl.exp < (distance + Ieee.ExponentBias)];
    };

  AlmostEqual: PUBLIC PROC [y, x: REAL, distance: INTEGER [-126..0]]
    RETURNS [BOOLEAN] = {
    fl: Ieee.SingleReal _ LOOPHOLE[(x - y)];
    oldexp: INTEGER _ MAX[
      LOOPHOLE[x, Ieee.SingleReal].exp, LOOPHOLE[y, Ieee.SingleReal].exp];
    RETURN[fl.exp < (distance + oldexp)];
    };

  smallnumber: REAL = 0.00000000005;
  PI: REAL = 3.14159265;
  twoPI: REAL = PI*2;
  PI3ovr2: REAL = 3*PI/2;
  PIovr2: REAL = 1.570796327;
  PI3ovr8: REAL = 1.1780972;
  PIovr4: REAL = .785398163;
  PIovr8: REAL = .392699;
  cp0: REAL = .999999953;
  cp1: REAL = -.499999053;
  cp2: REAL = .0416635847;
  cp3: REAL = -.001385370426;
  cp4: REAL = .0000231539317;
  radtodeg: REAL = 180/PI;
  degtorad: REAL = PI/180;
  recpi: REAL = 1/twoPI;

  tanPI16: REAL = .1989123673;
  x2: REAL = 1/(.414213562);
  x22: REAL = x2*x2 + 1;
  tan3PI16: REAL = .668178638;
  x4: REAL = 1/(1.0);
  x44: REAL = x4*x4 + 1;
  tan5PI16: REAL = 1.496605763;
  x6: REAL = 1/(2.41421356);
  x66: REAL = x6*x6 + 1;
  tan7PI16: REAL = 5.02733949;

  p0: REAL = .999999998;
  p1: REAL = -.333331726;
  p2: REAL = .1997952738;
  p3: REAL = -.134450639;

  END.
September 16, 1980  3:13 PM, Stewart; Bug in Exp
September 28, 1980  8:28 PM, Stewart; Add AlmostEqual and AlmostZero, format
November 7, 1980  3:33 PM, Stewart; fix AlmostZero
August 27, 1982 1:07 pm, Stewart; CEDAR