Project-Oberon-BlackBox/BlackBox/Po/Files/Math.Mod.txt

MODULE Math;   (*Standard functions;  NW 12.10.2013*)

  PROCEDURE sqrt*(x: REAL): REAL;
    CONST c1 = 0.70710680;   (* 1/sqrt(2) *)
      c2 = 0.590162067;
      c3 = 1.4142135;  (*sqrt(2)*)
    VAR s: REAL; e: INTEGER;
  BEGIN ASSERT(x >= 0.0);
    IF x > 0.0 THEN
      UNPK(x, e);
      s := c2*(x+c1);
      s := s + (x/s);
      s := 0.25*s + x/s;
      s := 0.5 * (s + x/s);
      IF ODD(e) THEN s := c3*s END ;
      PACK(s, e DIV 2)
    ELSE s := 0.0
    END ;
    RETURN s
  END sqrt;

  PROCEDURE exp*(x: REAL): REAL;
    CONST
      c1 = 1.4426951;  (*1/ln(2) *)
      p0 = 1.513864173E3;
      p1= 2.020170000E1;
      p2 = 2.309432127E-2;
      q0 = 4.368088670E3;
      q1 = 2.331782320E2;
    VAR n: INTEGER; p, y, yy: REAL;
  BEGIN y := c1*x;
    n := FLOOR(y + 0.5); y := y - FLT(n);
    yy := y*y;
    p := ((p2*yy + p1)*yy + p0)*y;
    p := p/((yy + q1)*yy + q0 - p) + 0.5;
    PACK(p, n+1); RETURN p
  END exp;

  PROCEDURE ln*(x: REAL): REAL;
    CONST c1 = 0.70710680;  (* 1/sqrt(2) *)
      c2 =  0.69314720;  (* ln(2) *)
      p0 = -9.01746917E1;
      p1 =  9.34639006E1;
      p2 = -1.83278704E1;
      q0 = -4.50873458E1;
      q1 =  6.76106560E1;
      q2 = -2.07334879E1;
    VAR e: INTEGER; xx, y: REAL;
  BEGIN ASSERT(x > 0.0); UNPK(x, e);
    IF x < c1 THEN x := x*2.0; e := e-1 END ;
    x := (x-1.0)/(x+1.0);
    xx := x;
    y := c2*FLT(e) + x*((p2*xx + p1)*xx + p0) / (((xx + q2)*xx + q1)*xx + q0);
    RETURN y
  END ln;

  PROCEDURE sin*(x: REAL): REAL;
    CONST
      c1 =  6.3661977E-1; (*2/pi*)
      p0 =  7.8539816E-1;
      p1 = -8.0745512E-2;
      p2 =  2.4903946E-3;
      p3 = -3.6576204E-5;
      p4 =  3.1336162E-7;
      p5 = -1.7571493E-9;
      p6 =  6.8771004E-12;
      q0 =  9.9999999E-1;
      q1 = -3.0842514E-1;
      q2 =  1.5854344E-2;
      q3 = -3.2599189E-4;
      q4 =  3.5908591E-6;
      q5 = -2.4609457E-8;
      q6 =  1.1363813E-10;
    VAR n: INTEGER; y, yy, f: REAL;
  BEGIN y := c1*x;
    IF y >= 0.0 THEN n := FLOOR(y + 0.5) ELSE n := FLOOR(y - 0.5) END ;
    y := (y - FLT(n)) * 2.0; yy := y*y;
    IF ODD(n) THEN f := (((((q6*yy + q5)*yy + q4)*yy + q3)*yy + q2)*yy + q1)*yy + q0
    ELSE f := ((((((p6*yy + p5)*yy + p4)*yy + p3)*yy + p2)*yy + p1)*yy + p0)*y
    END ;
    IF ODD(n DIV 2) THEN f := -f END ;
    RETURN f
  END sin;

  PROCEDURE cos*(x: REAL): REAL;
    CONST
      c1 =  6.3661977E-1; (*2/pi*)
      p0 =  7.8539816E-1;
      p1 = -8.0745512E-2;
      p2 =  2.4903946E-3;
      p3 = -3.6576204E-5;
      p4 =  3.1336162E-7;
      p5 = -1.7571493E-9;
      p6 = 6.8771004E-12;
      q0 =  9.9999999E-1;
      q1 = -3.0842514E-1;
      q2 =  1.5854344E-2;
      q3 = -3.2599189E-4;
      q4 =  3.5908591E-6;
      q5 = -2.4609457E-8;
      q6 =  1.1363813E-10;
    VAR n: INTEGER; y, yy, f: REAL;
  BEGIN y := c1*x;
    IF y >= 0.0 THEN n := FLOOR(y + 0.5) ELSE n := FLOOR(y - 0.5) END ;
    y := (y - FLT(n)) * 2.0; yy := y*y;
    IF ~ODD(n) THEN f := (((((q6*yy + q5)*yy + q4)*yy + q3)*yy + q2)*yy + q1)*yy + q0
    ELSE f := ((((((p6*yy + p5)*yy + p4)*yy + p3)*yy + p2)*yy + p1)*yy + p0)*y
    END ;
    IF ODD((n+1) DIV 2) THEN f := -f END ;
    RETURN f
  END cos;
END Math.