source: trunk/minix/lib/ack/libm2/Mathlib.mod@ 20

Last change on this file since 20 was 9, checked in by Mattia Monga, 14 years ago

Minix 3.1.2a
File size: 12.3 KB
RevLine 
[9]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-*)
7IMPLEMENTATION MODULE Mathlib;
8(*
9 Module: Mathematical functions
10 Author: Ceriel J.H. Jacobs
11 Version: $Header: /cvsup/minix/src/lib/ack/libm2/Mathlib.mod,v 1.1 2005/10/10 15:27:46 beng Exp $
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
574BEGIN
575 arctaninit := FALSE;
576END Mathlib.
Note: See TracBrowser for help on using the repository browser.