lib/ack/libm2/Mathlib.mod

   1 (*
   2   (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3   See the copyright notice in the ACK home directory, in the file "Copyright".
   4 *)
   5
   6 (*$R-*)
   7 IMPLEMENTATION MODULE Mathlib;
   8 (*
   9   Module:       Mathematical functions
  10   Author:       Ceriel J.H. Jacobs
  11   Version:      $Header$
  12 *)
  13
  14   FROM  EM IMPORT       FIF, FEF;
  15   FROM  Traps IMPORT    Message;
  16
  17   CONST
  18         OneRadianInDegrees      = 57.295779513082320876798155D;
  19         OneDegreeInRadians      =  0.017453292519943295769237D;
  20         OneOverSqrt2            = 0.70710678118654752440084436210484904D;
  21
  22   (* basic functions *)
  23
  24   PROCEDURE pow(x: REAL; i: INTEGER): REAL;
  25   BEGIN
  26         RETURN SHORT(longpow(LONG(x), i));
  27   END pow;
  28
  29   PROCEDURE longpow(x: LONGREAL; i: INTEGER): LONGREAL;
  30     VAR val: LONGREAL;
  31         ri: LONGREAL;
  32   BEGIN
  33         ri := FLOATD(i);
  34         IF x < 0.0D THEN
  35                 val := longexp(longln(-x) * ri);
  36                 IF ODD(i) THEN RETURN -val;
  37                 ELSE RETURN val;
  38                 END;
  39         ELSIF x = 0.0D THEN
  40                 RETURN 0.0D;
  41         ELSE
  42                 RETURN longexp(longln(x) * ri);
  43         END;
  44   END longpow;
  45
  46   PROCEDURE sqrt(x: REAL): REAL;
  47   BEGIN
  48         RETURN SHORT(longsqrt(LONG(x)));
  49   END sqrt;
  50
  51   PROCEDURE longsqrt(x: LONGREAL): LONGREAL;
  52     VAR
  53         temp: LONGREAL;
  54         exp, i: INTEGER;
  55   BEGIN
  56         IF x <= 0.0D THEN
  57                 IF x < 0.0D THEN
  58                         Message("sqrt: negative argument");
  59                         HALT
  60                 END;
  61                 RETURN 0.0D;
  62         END;
  63         temp := FEF(x,exp);
  64         (*
  65          * NOTE
  66          * this wont work on 1's comp
  67          *)
  68         IF ODD(exp) THEN
  69                 temp := 2.0D * temp;
  70                 DEC(exp);
  71         END;
  72         temp := 0.5D*(1.0D + temp);
  73
  74         WHILE exp > 28 DO
  75                 temp := temp * 16384.0D;
  76                 exp := exp - 28;
  77         END;
  78         WHILE exp < -28 DO
  79                 temp := temp / 16384.0D;
  80                 exp := exp + 28;
  81         END;
  82         WHILE exp >= 2 DO
  83                 temp := temp * 2.0D;
  84                 exp := exp - 2;
  85         END;
  86         WHILE exp <= -2 DO
  87                 temp := temp / 2.0D;
  88                 exp := exp + 2;
  89         END;
  90         FOR i := 0 TO 5 DO
  91                 temp := 0.5D*(temp + x/temp);
  92         END;
  93         RETURN temp;
  94   END longsqrt;
  95
  96   PROCEDURE ldexp(x:LONGREAL; n: INTEGER): LONGREAL;
  97   BEGIN
  98         WHILE n >= 16 DO
  99                 x := x * 65536.0D;
 100                 n := n - 16;
 101         END;
 102         WHILE n > 0 DO
 103                 x := x * 2.0D;
 104                 DEC(n);
 105         END;
 106         WHILE n <= -16 DO
 107                 x := x / 65536.0D;
 108                 n := n + 16;
 109         END;
 110         WHILE n < 0 DO
 111                 x := x / 2.0D;
 112                 INC(n);
 113         END;
 114         RETURN x;
 115   END ldexp;
 116
 117   PROCEDURE exp(x: REAL): REAL;
 118   BEGIN
 119         RETURN SHORT(longexp(LONG(x)));
 120   END exp;
 121
 122   PROCEDURE longexp(x: LONGREAL): LONGREAL;
 123   (*    Algorithm and coefficients from:
 124                 "Software manual for the elementary functions"
 125                 by W.J. Cody and W. Waite, Prentice-Hall, 1980
 126   *)
 127     CONST
 128         p0 = 0.25000000000000000000D+00;
 129         p1 = 0.75753180159422776666D-02;
 130         p2 = 0.31555192765684646356D-04;
 131         q0 = 0.50000000000000000000D+00;
 132         q1 = 0.56817302698551221787D-01;
 133         q2 = 0.63121894374398503557D-03;
 134         q3 = 0.75104028399870046114D-06;
 135
 136     VAR
 137         neg: BOOLEAN;
 138         n: INTEGER;
 139         xn, g, x1, x2: LONGREAL;
 140   BEGIN
 141         neg := x < 0.0D;
 142         IF neg THEN
 143                 x := -x;
 144         END;
 145         n := TRUNC(x/longln2 + 0.5D);
 146         xn := FLOATD(n);
 147         x1 := FLOATD(TRUNCD(x));
 148         x2 := x - x1;
 149         g := ((x1 - xn * 0.693359375D)+x2) - xn * (-2.1219444005469058277D-4);
 150         IF neg THEN
 151                 g := -g;
 152                 n := -n;
 153         END;
 154         xn := g*g;
 155         x := g*((p2*xn+p1)*xn+p0);
 156         INC(n);
 157         RETURN ldexp(0.5D + x/((((q3*xn+q2)*xn+q1)*xn+q0) - x), n);
 158   END longexp;
 159
 160   PROCEDURE ln(x: REAL): REAL;  (* natural log *)
 161   BEGIN
 162         RETURN SHORT(longln(LONG(x)));
 163   END ln;
 164
 165   PROCEDURE longln(x: LONGREAL): LONGREAL;      (* natural log *)
 166   (*    Algorithm and coefficients from:
 167                 "Software manual for the elementary functions"
 168                 by W.J. Cody and W. Waite, Prentice-Hall, 1980
 169    *)
 170     CONST
 171         p0 = -0.64124943423745581147D+02;
 172         p1 =  0.16383943563021534222D+02;
 173         p2 = -0.78956112887491257267D+00;
 174         q0 = -0.76949932108494879777D+03;
 175         q1 =  0.31203222091924532844D+03;
 176         q2 = -0.35667977739034646171D+02;
 177         q3 =  1.0D;
 178     VAR
 179         exp: INTEGER;
 180         z, znum, zden, w: LONGREAL;
 181
 182   BEGIN
 183         IF x <= 0.0D THEN
 184                 Message("ln: argument <= 0");
 185                 HALT
 186         END;
 187         x := FEF(x, exp);
 188         IF x > OneOverSqrt2 THEN
 189                 znum := (x - 0.5D) - 0.5D;
 190                 zden := x * 0.5D + 0.5D;
 191         ELSE
 192                 znum := x - 0.5D;
 193                 zden := znum * 0.5D + 0.5D;
 194                 DEC(exp);
 195         END;
 196         z := znum / zden;
 197         w := z * z;
 198         x := z + z * w * (((p2*w+p1)*w+p0)/(((q3*w+q2)*w+q1)*w+q0));
 199         z := FLOATD(exp);
 200         x := x + z * (-2.121944400546905827679D-4);
 201         RETURN x + z * 0.693359375D;
 202   END longln;
 203
 204   PROCEDURE log(x: REAL): REAL; (* log with base 10 *)
 205   BEGIN
 206         RETURN SHORT(longlog(LONG(x)));
 207   END log;
 208
 209   PROCEDURE longlog(x: LONGREAL): LONGREAL;     (* log with base 10 *)
 210   BEGIN
 211         RETURN longln(x)/longln10;
 212   END longlog;
 213
 214   (* trigonometric functions; arguments in radians *)
 215
 216   PROCEDURE sin(x: REAL): REAL;
 217   BEGIN
 218         RETURN SHORT(longsin(LONG(x)));
 219   END sin;
 220
 221   PROCEDURE sinus(x: LONGREAL; cosflag: BOOLEAN) : LONGREAL;
 222   (*    Algorithm and coefficients from:
 223                 "Software manual for the elementary functions"
 224                 by W.J. Cody and W. Waite, Prentice-Hall, 1980
 225   *)
 226     CONST
 227         r0 = -0.16666666666666665052D+00;
 228         r1 =  0.83333333333331650314D-02;
 229         r2 = -0.19841269841201840457D-03;
 230         r3 =  0.27557319210152756119D-05;
 231         r4 = -0.25052106798274584544D-07;
 232         r5 =  0.16058936490371589114D-09;
 233         r6 = -0.76429178068910467734D-12;
 234         r7 =  0.27204790957888846175D-14;
 235         A1 =  3.1416015625D;
 236         A2 = -8.908910206761537356617D-6;
 237     VAR
 238         x1, x2, y : LONGREAL;
 239         neg : BOOLEAN;
 240   BEGIN
 241         IF x < 0.0D THEN
 242                 neg := TRUE;
 243                 x := -x
 244         ELSE    neg := FALSE
 245         END;
 246         IF cosflag THEN
 247                 neg := FALSE;
 248                 y := longhalfpi + x
 249         ELSE
 250                 y := x
 251         END;
 252         y := y / longpi + 0.5D;
 253
 254         IF FIF(y, 1.0D, y) < 0.0D THEN ; END;
 255         IF FIF(y, 0.5D, x1) # 0.0D THEN neg := NOT neg END;
 256         IF cosflag THEN y := y - 0.5D END;
 257         x2 := FIF(x, 1.0, x1);
 258         x := x1 - y * A1;
 259         x := x + x2;
 260         x := x - y * A2;
 261
 262         IF x < 0.0D THEN
 263                 neg := NOT neg;
 264                 x := -x
 265         END;
 266         y := x * x;
 267         x := x + x * y * (((((((r7*y+r6)*y+r5)*y+r4)*y+r3)*y+r2)*y+r1)*y+r0);
 268         IF neg THEN RETURN -x END;
 269         RETURN x;
 270   END sinus;
 271
 272   PROCEDURE longsin(x: LONGREAL): LONGREAL;
 273   BEGIN
 274         RETURN sinus(x, FALSE);
 275   END longsin;
 276
 277   PROCEDURE cos(x: REAL): REAL;
 278   BEGIN
 279         RETURN SHORT(longcos(LONG(x)));
 280   END cos;
 281
 282   PROCEDURE longcos(x: LONGREAL): LONGREAL;
 283   BEGIN
 284         IF x < 0.0D THEN x := -x; END;
 285         RETURN sinus(x, TRUE);
 286   END longcos;
 287
 288   PROCEDURE tan(x: REAL): REAL;
 289   BEGIN
 290         RETURN SHORT(longtan(LONG(x)));
 291   END tan;
 292
 293   PROCEDURE longtan(x: LONGREAL): LONGREAL;
 294   (*    Algorithm and coefficients from:
 295                 "Software manual for the elementary functions"
 296                 by W.J. Cody and W. Waite, Prentice-Hall, 1980
 297   *)
 298
 299     CONST
 300         p1 = -0.13338350006421960681D+00;
 301         p2 =  0.34248878235890589960D-02;
 302         p3 = -0.17861707342254426711D-04;
 303
 304         q0 =  1.0D;
 305         q1 = -0.46671683339755294240D+00;
 306         q2 =  0.25663832289440112864D-01;
 307         q3 = -0.31181531907010027307D-03;
 308         q4 =  0.49819433993786512270D-06;
 309
 310         A1 =  1.57080078125D;
 311         A2 = -4.454455103380768678308D-06;
 312
 313     VAR y, x1, x2: LONGREAL;
 314         negative: BOOLEAN;
 315         invert: BOOLEAN;
 316   BEGIN
 317         negative := x < 0.0D;
 318         y := x / longhalfpi + 0.5D;
 319
 320         (*      Use extended precision to calculate reduced argument.
 321                 Here we used 12 bits of the mantissa for a1.
 322                 Also split x in integer part x1 and fraction part x2.
 323         *)
 324         IF FIF(y, 1.0D, y) < 0.0D THEN ; END;
 325         invert := FIF(y, 0.5D, x1) # 0.0D;
 326         x2 := FIF(x, 1.0D, x1);
 327         x := x1 - y * A1;
 328         x := x + x2;
 329         x := x - y * A2;
 330
 331         y := x * x;
 332         x := x + x * y * ((p3*y+p2)*y+p1);
 333         y := (((q4*y+q3)*y+q2)*y+q1)*y+q0;
 334         IF negative THEN x := -x END;
 335         IF invert THEN RETURN -y/x END;
 336         RETURN x/y;
 337   END longtan;
 338
 339   PROCEDURE arcsin(x: REAL): REAL;
 340   BEGIN
 341         RETURN SHORT(longarcsin(LONG(x)));
 342   END arcsin;
 343
 344   PROCEDURE arcsincos(x: LONGREAL; cosfl: BOOLEAN): LONGREAL;
 345     CONST
 346         p0 = -0.27368494524164255994D+02;
 347         p1 =  0.57208227877891731407D+02;
 348         p2 = -0.39688862997540877339D+02;
 349         p3 =  0.10152522233806463645D+02;
 350         p4 = -0.69674573447350646411D+00;
 351
 352         q0 = -0.16421096714498560795D+03;
 353         q1 =  0.41714430248260412556D+03;
 354         q2 = -0.38186303361750149284D+03;
 355         q3 =  0.15095270841030604719D+03;
 356         q4 = -0.23823859153670238830D+02;
 357         q5 =  1.0D;
 358     VAR
 359         negative : BOOLEAN;
 360         big: BOOLEAN;
 361         g: LONGREAL;
 362   BEGIN
 363         negative := x < 0.0D;
 364         IF negative THEN x := -x; END;
 365         IF x > 0.5D THEN
 366                 big := TRUE;
 367                 IF x > 1.0D THEN
 368                         Message("arcsin or arccos: argument > 1");
 369                         HALT
 370                 END;
 371                 g := 0.5D - 0.5D * x;
 372                 x := -longsqrt(g);
 373                 x := x + x;
 374         ELSE
 375                 big := FALSE;
 376                 g := x * x;
 377         END;
 378         x := x + x * g *
 379           ((((p4*g+p3)*g+p2)*g+p1)*g+p0)/(((((q5*g+q4)*g+q3)*g+q2)*g+q1)*g+q0);
 380         IF cosfl AND NOT negative THEN x := -x END;
 381         IF cosfl = NOT big THEN
 382                 x := (x + longquartpi) + longquartpi;
 383         ELSIF cosfl AND negative AND big THEN
 384                 x := (x + longhalfpi) + longhalfpi;
 385         END;
 386         IF negative AND NOT cosfl THEN x := -x END;
 387         RETURN x;
 388   END arcsincos;
 389
 390   PROCEDURE longarcsin(x: LONGREAL): LONGREAL;
 391   BEGIN
 392         RETURN arcsincos(x, FALSE);
 393   END longarcsin;
 394
 395   PROCEDURE arccos(x: REAL): REAL;
 396   BEGIN
 397         RETURN SHORT(longarccos(LONG(x)));
 398   END arccos;
 399
 400   PROCEDURE longarccos(x: LONGREAL): LONGREAL;
 401   BEGIN
 402         RETURN arcsincos(x, TRUE);
 403   END longarccos;
 404
 405   PROCEDURE arctan(x: REAL): REAL;
 406   BEGIN
 407         RETURN SHORT(longarctan(LONG(x)));
 408   END arctan;
 409
 410   VAR A: ARRAY[0..3] OF LONGREAL;
 411       arctaninit: BOOLEAN;
 412
 413   PROCEDURE longarctan(x: LONGREAL): LONGREAL;
 414   (*    Algorithm and coefficients from:
 415                 "Software manual for the elementary functions"
 416                 by W.J. Cody and W. Waite, Prentice-Hall, 1980
 417   *)
 418     CONST
 419         p0 = -0.13688768894191926929D+02;
 420         p1 = -0.20505855195861651981D+02;
 421         p2 = -0.84946240351320683534D+01;
 422         p3 = -0.83758299368150059274D+00;
 423         q0 =  0.41066306682575781263D+02;
 424         q1 =  0.86157349597130242515D+02;
 425         q2 =  0.59578436142597344465D+02;
 426         q3 =  0.15024001160028576121D+02;
 427         q4 =  1.0D;
 428     VAR
 429         g: LONGREAL;
 430         neg: BOOLEAN;
 431         n: INTEGER;
 432   BEGIN
 433         IF NOT arctaninit THEN
 434                 arctaninit := TRUE;
 435                 A[0] := 0.0D;
 436                 A[1] := 0.52359877559829887307710723554658381D; (* p1/6 *)
 437                 A[2] := longhalfpi;
 438                 A[3] := 1.04719755119659774615421446109316763D; (* pi/3 *)
 439         END;
 440         neg := FALSE;
 441         IF x < 0.0D THEN
 442                 neg := TRUE;
 443                 x := -x;
 444         END;
 445         IF x > 1.0D THEN
 446                 x := 1.0D/x;
 447                 n := 2
 448         ELSE
 449                 n := 0
 450         END;
 451         IF x > 0.26794919243112270647D (* 2-sqrt(3) *) THEN
 452                 INC(n);
 453                 x := (((0.73205080756887729353D*x-0.5D)-0.5D)+x)/
 454                         (1.73205080756887729353D + x);
 455         END;
 456         g := x*x;
 457         x := x + x * g * (((p3*g+p2)*g+p1)*g+p0) / ((((q4*g+q3)*g+q2)*g+q1)*g+q0);
 458         IF n > 1 THEN x := -x END;
 459         x := x + A[n];
 460         IF neg THEN RETURN -x; END;
 461         RETURN x;
 462   END longarctan;
 463
 464   (* hyperbolic functions *)
 465   (* The C math library has better implementations for some of these, but
 466      they depend on some properties of the floating point implementation,
 467      and, for now, we don't want that in the Modula-2 system.
 468   *)
 469
 470   PROCEDURE sinh(x: REAL): REAL;
 471   BEGIN
 472         RETURN SHORT(longsinh(LONG(x)));
 473   END sinh;
 474
 475   PROCEDURE longsinh(x: LONGREAL): LONGREAL;
 476     VAR expx: LONGREAL;
 477   BEGIN
 478         expx := longexp(x);
 479         RETURN (expx - 1.0D/expx)/2.0D;
 480   END longsinh;
 481
 482   PROCEDURE cosh(x: REAL): REAL;
 483   BEGIN
 484         RETURN SHORT(longcosh(LONG(x)));
 485   END cosh;
 486
 487   PROCEDURE longcosh(x: LONGREAL): LONGREAL;
 488     VAR expx: LONGREAL;
 489   BEGIN
 490         expx := longexp(x);
 491         RETURN (expx + 1.0D/expx)/2.0D;
 492   END longcosh;
 493
 494   PROCEDURE tanh(x: REAL): REAL;
 495   BEGIN
 496         RETURN SHORT(longtanh(LONG(x)));
 497   END tanh;
 498
 499   PROCEDURE longtanh(x: LONGREAL): LONGREAL;
 500     VAR expx: LONGREAL;
 501   BEGIN
 502         expx := longexp(x);
 503         RETURN (expx - 1.0D/expx) / (expx + 1.0D/expx);
 504   END longtanh;
 505
 506   PROCEDURE arcsinh(x: REAL): REAL;
 507   BEGIN
 508         RETURN SHORT(longarcsinh(LONG(x)));
 509   END arcsinh;
 510
 511   PROCEDURE longarcsinh(x: LONGREAL): LONGREAL;
 512     VAR neg: BOOLEAN;
 513   BEGIN
 514         neg := FALSE;
 515         IF x < 0.0D THEN
 516                 neg := TRUE;
 517                 x := -x;
 518         END;
 519         x := longln(x + longsqrt(x*x+1.0D));
 520         IF neg THEN RETURN -x; END;
 521         RETURN x;
 522   END longarcsinh;
 523
 524   PROCEDURE arccosh(x: REAL): REAL;
 525   BEGIN
 526         RETURN SHORT(longarccosh(LONG(x)));
 527   END arccosh;
 528
 529   PROCEDURE longarccosh(x: LONGREAL): LONGREAL;
 530   BEGIN
 531         IF x < 1.0D THEN
 532                 Message("arccosh: argument < 1");
 533                 HALT
 534         END;
 535         RETURN longln(x + longsqrt(x*x - 1.0D));
 536   END longarccosh;
 537
 538   PROCEDURE arctanh(x: REAL): REAL;
 539   BEGIN
 540         RETURN SHORT(longarctanh(LONG(x)));
 541   END arctanh;
 542
 543   PROCEDURE longarctanh(x: LONGREAL): LONGREAL;
 544   BEGIN
 545         IF (x <= -1.0D) OR (x >= 1.0D) THEN
 546                 Message("arctanh: ABS(argument) >= 1");
 547                 HALT
 548         END;
 549         RETURN longln((1.0D + x)/(1.0D - x)) / 2.0D;
 550   END longarctanh;
 551
 552   (* conversions *)
 553
 554   PROCEDURE RadianToDegree(x: REAL): REAL;
 555   BEGIN
 556         RETURN SHORT(longRadianToDegree(LONG(x)));
 557   END RadianToDegree;
 558
 559   PROCEDURE longRadianToDegree(x: LONGREAL): LONGREAL;
 560   BEGIN
 561         RETURN x * OneRadianInDegrees;
 562   END longRadianToDegree;
 563
 564   PROCEDURE DegreeToRadian(x: REAL): REAL;
 565   BEGIN
 566         RETURN SHORT(longDegreeToRadian(LONG(x)));
 567   END DegreeToRadian;
 568
 569   PROCEDURE longDegreeToRadian(x: LONGREAL): LONGREAL;
 570   BEGIN
 571         RETURN x * OneDegreeInRadians;
 572   END longDegreeToRadian;
 573
 574 BEGIN
 575         arctaninit := FALSE;
 576 END Mathlib.