DEADSOFTWARE

936cf15d137e8be0a72bc13cd78e03c21166297f
[bbcp.git] / Trurl-based / System / Mod / Math.txt
1 MODULE Math;
3 (* THIS IS TEXT COPY OF BlackBox 1.6-rc6 System/Mod/Math.odc *)
4 (* DO NOT EDIT *)
6 IMPORT SYSTEM;
8 VAR eps, e: REAL;
11 (* code procedures for 80387 math coprocessor *)
13 PROCEDURE [code] FLD (x: REAL);
14 PROCEDURE [code] TOP (): REAL;
15 PROCEDURE [code] FSW (): INTEGER 0DFH, 0E0H;
16 PROCEDURE [code] FSWs (): SET 0DFH, 0E0H;
17 PROCEDURE [code] ST0 (): REAL 0D9H, 0C0H;
18 PROCEDURE [code] ST1 (): REAL 0D9H, 0C1H;
20 PROCEDURE [code] FXCH 0D9H, 0C9H;
21 PROCEDURE [code] FLDst0 0D9H, 0C0H; (* doublicate st[0] *)
22 PROCEDURE [code] FSTPst0 0DDH, 0D8H; (* remove st[0] *)
23 PROCEDURE [code] FSTPst1 0DDH, 0D9H; (* remove st[1] *)
24 PROCEDURE [code] FSTPDe 0DBH, 05DH, 0F4H; (* FSTPD -12[FP] *) (* COMPILER DEPENDENT *)
25 PROCEDURE [code] WAIT 09BH;
26 PROCEDURE [code] FNOP 0D9H, 0D0H;
28 PROCEDURE [code] FLD0 0D9H, 0EEH;
29 PROCEDURE [code] FLD1 0D9H, 0E8H;
30 PROCEDURE [code] FLDPI 0D9H, 0EBH;
31 PROCEDURE [code] FLDLN2 0D9H, 0EDH;
32 PROCEDURE [code] FLDLG2 0D9H, 0ECH;
33 PROCEDURE [code] FLDL2E 0D9H, 0EAH;
35 PROCEDURE [code] FADD 0DEH, 0C1H;
36 PROCEDURE [code] FADDst0 0D8H, 0C0H;
37 PROCEDURE [code] FSUB 0DEH, 0E9H;
38 PROCEDURE [code] FSUBn 0DCH, 0E9H; (* no pop *)
39 PROCEDURE [code] FSUBR 0DEH, 0E1H;
40 PROCEDURE [code] FSUBst1 0D8H, 0E1H;
41 PROCEDURE [code] FMUL 0DEH, 0C9H;
42 PROCEDURE [code] FMULst0 0D8H, 0C8H;
43 PROCEDURE [code] FMULst1st0 0DCH, 0C9H;
44 PROCEDURE [code] FDIV 0DEH, 0F9H;
45 PROCEDURE [code] FDIVR 0DEH, 0F1H;
46 PROCEDURE [code] FDIVRst1 0D8H, 0F9H;
47 PROCEDURE [code] FCHS 0D9H, 0E0H;
49 PROCEDURE [code] FCOM 0D8H, 0D1H;
50 PROCEDURE [code] FSWax 0DFH, 0E0H;
51 PROCEDURE [code] SAHF 09EH;
52 PROCEDURE [code] JBE4 076H, 004H;
53 PROCEDURE [code] JAE4 073H, 004H;
55 PROCEDURE [code] FRNDINT 0D9H, 0FCH;
56 PROCEDURE [code] FSCALE 0D9H, 0FDH; (* st[0] * 2^FLOOR(st[1]) *)
57 PROCEDURE [code] FXTRACT 0D9H, 0F4H; (* exp -> st[1]; mant -> st[0] *)
58 PROCEDURE [code] FXAM 0D9H, 0E5H;
60 PROCEDURE [code] FSQRT 0D9H, 0FAH; (* st[0] >= 0 *)
61 PROCEDURE [code] FSIN 0D9H, 0FEH; (* |st[0]| < 2^63 *)
62 PROCEDURE [code] FCOS 0D9H, 0FFH; (* |st[0]| < 2^63 *)
63 PROCEDURE [code] FTAN 0D9H, 0F2H; (* |st[0]| < 2^63 *)
64 PROCEDURE [code] FATAN 0D9H, 0F3H; (* atan2(st[1], st[0]) *)
65 PROCEDURE [code] FYL2X 0D9H, 0F1H; (* st[1] * log2(st[0]), st[0] > 0 *)
66 PROCEDURE [code] FYL2XP1 0D9H, 0F9H; (* st[1] * log2(1 + st[0]), |st[0]| < 1-sqrt(2)/2 *)
67 PROCEDURE [code] F2XM1 0D9H, 0F0H; (* 2^st[0] - 1, |st[0]| <= 1 *)
70 PROCEDURE IsNan (x: REAL): BOOLEAN;
71 BEGIN
72 FLD(x); FXAM; FSTPst0; WAIT; RETURN FSWs() * {8, 10} = {8}
73 END IsNan;
76 (* sin, cos, tan argument reduction *)
78 PROCEDURE Reduce;
79 BEGIN
80 FXAM; WAIT;
81 IF ~(8 IN FSWs()) & (ABS(ST0()) > 1.0E18) THEN
82 (* to be completed *)
83 FSTPst0; FLD0
84 END;
85 END Reduce;
88 (** REAL precision **)
90 PROCEDURE Pi* (): REAL;
91 BEGIN
92 FLDPI; RETURN TOP()
93 END Pi;
95 PROCEDURE Eps* (): REAL;
96 BEGIN
97 RETURN eps
98 END Eps;
101 PROCEDURE Sqrt* (x: REAL): REAL;
102 BEGIN
103 (* 20, argument of Sqrt must not be negative *)
104 FLD(x); FSQRT; WAIT; RETURN TOP()
105 END Sqrt;
108 PROCEDURE Exp* (x: REAL): REAL;
109 BEGIN
110 (* 2 ^ (x * 1/ln(2)) *)
111 FLD(x); FLDL2E; FMUL;
112 IF ABS(ST0()) = INF THEN FLD1
113 ELSE FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD
114 END;
115 FSCALE; FSTPst1; RETURN TOP()
116 END Exp;
118 PROCEDURE Ln* (x: REAL): REAL;
119 BEGIN
120 (* 20, argument of Ln must not be negative *)
121 (* ln(2) * ld(x) *)
122 FLDLN2; FLD(x); FYL2X; WAIT; RETURN TOP()
123 END Ln;
125 PROCEDURE Log* (x: REAL): REAL;
126 BEGIN
127 (* 20, argument of Log must not be negative *)
128 (* log(2) * ld(x) *)
129 FLDLG2; FLD(x); FYL2X; WAIT; RETURN TOP()
130 END Log;
132 PROCEDURE Power* (x, y: REAL): REAL;
133 BEGIN
134 ASSERT(x >= 0, 20);
135 ASSERT((x # 0.0) OR (y # 0.0), 21);
136 ASSERT((x # INF) OR (y # 0.0), 22);
137 ASSERT((x # 1.0) OR (ABS(y) # INF), 23);
138 (* 2 ^ (y * ld(x)) *)
139 FLD(y); FLD(x); FYL2X;
140 IF ABS(ST0()) = INF THEN FLD1
141 ELSE FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD
142 END;
143 FSCALE; FSTPst1; WAIT; RETURN TOP()
144 END Power;
146 PROCEDURE IntPower* (x: REAL; n: INTEGER): REAL;
147 BEGIN
148 FLD1; FLD(x);
149 IF n = MIN(INTEGER) THEN RETURN IntPower(x, n + 1) / x END;
150 IF n <= 0 THEN FDIVRst1; (* 1 / x *) n := -n END;
151 WHILE n > 0 DO
152 IF ODD(n) THEN FMULst1st0; (* y := y * x *) DEC(n)
153 ELSE FMULst0; (* x := x * x *) n := n DIV 2
154 END
155 END;
156 FSTPst0; RETURN TOP()
157 END IntPower;
160 PROCEDURE Sin* (x: REAL): REAL;
161 BEGIN
162 (* 20, ABS(x) # INF *)
163 FLD(x); Reduce; FSIN; WAIT; RETURN TOP()
164 END Sin;
166 PROCEDURE Cos* (x: REAL): REAL;
167 BEGIN
168 (* 20, ABS(x) # INF *)
169 FLD(x); Reduce; FCOS; WAIT; RETURN TOP()
170 END Cos;
172 PROCEDURE Tan* (x: REAL): REAL;
173 BEGIN
174 (* 20, ABS(x) # INF *)
175 FLD(x); Reduce; FTAN; FSTPst0; WAIT; RETURN TOP()
176 END Tan;
178 PROCEDURE ArcSin* (x: REAL): REAL;
179 BEGIN
180 (* 20, -1.0 <= x <= 1.0 *)
181 (* atan2(x, sqrt(1 - x*x)) *)
182 FLD(x); FLDst0; FMULst0; FLD1; FSUBR; FSQRT; FNOP; FATAN; WAIT; RETURN TOP()
183 END ArcSin;
185 PROCEDURE ArcCos* (x: REAL): REAL;
186 BEGIN
187 (* 20, -1.0 <= x <= 1.0 *)
188 (* atan2(sqrt(1 - x*x), x) *)
189 FLD(x); FMULst0; FLD1; FSUBR; FSQRT; FLD(x); FATAN; WAIT; RETURN TOP()
190 END ArcCos;
192 PROCEDURE ArcTan* (x: REAL): REAL;
193 BEGIN
194 (* atan2(x, 1) *)
195 FLD(x); FLD1; FATAN; RETURN TOP()
196 END ArcTan;
198 PROCEDURE ArcTan2* (y, x: REAL): REAL;
199 BEGIN
200 ASSERT((y # 0) OR (x # 0), 20);
201 ASSERT((ABS(y) # INF) OR (ABS(x) # INF), 21);
202 FLD(y); FLD(x); FATAN; WAIT; RETURN TOP()
203 END ArcTan2;
206 PROCEDURE Sinh* (x: REAL): REAL;
207 BEGIN
208 (* IF IsNan(x) THEN RETURN x END; *)
209 (* abs(x) * 1/ln(2) *)
210 FLD(ABS(x)); FLDL2E; FMUL;
211 IF ST0() < 0.5 THEN
212 (* (2^z - 1) + (2^z - 1) / ((2^z - 1) + 1) *)
213 F2XM1; FLDst0; FLDst0; FLD1; FADD; FDIV; FADD
214 ELSIF ST0() # INF THEN
215 (* 2^z - 1 / 2^z *)
216 FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD; FSCALE;
217 FSTPst1; FLDst0; FLD1; FDIVR; FSUB
218 END;
219 IF x < 0 THEN FCHS END;
220 RETURN TOP() * 0.5
221 END Sinh;
223 PROCEDURE Cosh* (x: REAL): REAL;
224 BEGIN
225 (* IF IsNan(x) THEN RETURN x END; *)
226 (* 2^(abs(x) * 1/ln(2)) *)
227 FLD(ABS(x));
228 IF ST0() # INF THEN
229 FLDL2E; FMUL; FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD; FSCALE;
230 FSTPst1;
231 (* z + 1/z *)
232 FLDst0; FLD1; FDIVR; FADD
233 END;
234 RETURN TOP() * 0.5
235 END Cosh;
237 PROCEDURE Tanh* (x: REAL): REAL;
238 BEGIN
239 (* IF IsNan(x) THEN RETURN x END; *)
240 (* abs(x) * 1/ln(2) * 2 *)
241 FLD(ABS(x)); FLDL2E; FMUL; FADDst0;
242 IF ST0() < 0.5 THEN
243 (* (2^z - 1) / (2^z + 1) *)
244 F2XM1; FLDst0; FLD(2); FADD; FDIV
245 ELSIF ST0() < 65 THEN
246 (* 1 - 2 / (2^z + 1) *)
247 FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD; FSCALE;
248 FSTPst1; FLD1; FADD; FLD(2); FDIVR; FLD1; FSUBR
249 ELSE
250 FSTPst0; FLD1
251 END;
252 IF x < 0 THEN FCHS END;
253 RETURN TOP()
254 END Tanh;
256 PROCEDURE ArcSinh* (x: REAL): REAL;
257 BEGIN
258 (* IF IsNan(x) THEN RETURN x END; *)
259 (* x*x *)
260 FLDLN2; FLD(ABS(x)); FLDst0; FMULst0;
261 IF ST0() < 0.067 THEN
262 (* ln(2) * ld(1 + x*x / (sqrt(x*x + 1) + 1) + x) *)
263 FLDst0; FLD1; FADD; FSQRT; FLD1; FADD; FDIV; FADD; FYL2XP1
264 ELSE
265 (* ln(2) * ld(x + sqrt(x*x + 1)) *)
266 FLD1; FADD; FSQRT; FADD; FYL2X
267 END;
268 IF x < 0 THEN FCHS END;
269 RETURN TOP()
270 END ArcSinh;
272 PROCEDURE ArcCosh* (x: REAL): REAL;
273 BEGIN
274 (* 20, x >= 1.0 *)
275 (* IF IsNan(x) THEN RETURN x END; *)
276 (* ln(2) * ld(x + sqrt(x*x - 1)) *)
277 FLDLN2; FLD(x); FLDst0; FMULst0; FLD1; FSUB; FSQRT; FADD; FYL2X; WAIT; RETURN TOP()
278 END ArcCosh;
280 PROCEDURE ArcTanh* (x: REAL): REAL;
281 BEGIN
282 (* 20, -1.0 <= x <= 1.0 *)
283 (* IF IsNan(x) THEN RETURN x END; *)
284 (* |x| *)
285 FLDLN2; FLD(ABS(x));
286 IF ST0() < 0.12 THEN
287 (* ln(2) * ld(1 + 2*x / (1 - x)) *)
288 FLDst0; FLD1; FSUBR; FDIV; FADDst0; FYL2XP1
289 ELSE
290 (* ln(2) * ld((1 + x) / (1 - x)) *)
291 FLDst0; FLD1; FADD; FXCH; FLD1; FSUBR; FDIV; FNOP; FYL2X
292 END;
293 IF x < 0 THEN FCHS END;
294 WAIT;
295 RETURN TOP() * 0.5
296 END ArcTanh;
299 PROCEDURE Floor* (x: REAL): REAL;
300 BEGIN
301 FLD(x); FLDst0; FRNDINT; FCOM; FSWax; FSTPst1; SAHF; JBE4; FLD1; FSUB; RETURN TOP()
302 END Floor;
304 PROCEDURE Ceiling* (x: REAL): REAL;
305 BEGIN
306 FLD(x); FLDst0; FRNDINT; FCOM; FSWax; FSTPst1; SAHF; JAE4; FLD1; FADD; RETURN TOP()
307 END Ceiling;
309 PROCEDURE Round* (x: REAL): REAL;
310 BEGIN
311 FLD(x);
312 IF ABS(ST0()) = INF THEN RETURN TOP() END;
313 FLDst0; FRNDINT; FSUBn; FXCH;
314 IF TOP() = 0.5 THEN FLD1; FADD END;
315 RETURN TOP()
316 END Round;
318 PROCEDURE Trunc* (x: REAL): REAL;
319 BEGIN
320 FLD(x); FLDst0; FRNDINT;
321 IF ST1() >= 0 THEN
322 FCOM; FSWax; FSTPst1; SAHF; JBE4; FLD1; FSUB
323 ELSE
324 FCOM; FSWax; FSTPst1; SAHF; JAE4; FLD1; FADD
325 END;
326 RETURN TOP()
327 END Trunc;
329 PROCEDURE Frac* (x: REAL): REAL;
330 BEGIN
331 (* 20, x # INF & x # -INF *)
332 FLD(x); FLDst0; FRNDINT;
333 IF ST1() >= 0 THEN
334 FCOM; FSWax; SAHF; JBE4; FLD1; FSUB
335 ELSE
336 FCOM; FSWax; SAHF; JAE4; FLD1; FADD
337 END;
338 FSUB; WAIT; RETURN TOP()
339 END Frac;
342 PROCEDURE Sign* (x: REAL): REAL;
343 BEGIN
344 FLD(x); FXAM; WAIT;
345 CASE FSW() DIV 256 MOD 8 OF
346 | 0, 2: FSTPst0; RETURN 0.0
347 | 1, 4, 5: FSTPst0; RETURN 1.0
348 | 3, 6, 7: FSTPst0; RETURN -1.0
349 END
350 END Sign;
352 PROCEDURE Mantissa* (x: REAL): REAL;
353 BEGIN
354 FLD(x); FXAM; WAIT;
355 CASE FSW() DIV 256 MOD 8 OF
356 | 4, 6: FXTRACT; FSTPst1; RETURN TOP()
357 | 0, 2: FSTPst0; RETURN 0.0 (* zero *)
358 | 5: FSTPst0; RETURN 1.0 (* inf *)
359 | 7: FSTPst0; RETURN -1.0 (* -inf *)
360 | 1: FSTPst0; RETURN 1.5 (* nan *)
361 | 3: FSTPst0; RETURN -1.5 (* -nan *)
362 END
363 END Mantissa;
365 PROCEDURE Exponent* (x: REAL): INTEGER; (* COMPILER DEPENDENT *)
366 VAR e: INTEGER; (* e is set by FSTPDe! *)
367 BEGIN
368 FLD(x); FXAM; WAIT;
369 CASE FSW() DIV 256 MOD 8 OF
370 | 4, 6: FXTRACT; FSTPst0; FSTPDe; WAIT; RETURN e
371 | 0, 2: FSTPst0; RETURN 0 (* zero *)
372 | 1, 3, 5, 7: FSTPst0; RETURN MAX(INTEGER) (* inf or nan*)
373 END
374 END Exponent;
376 PROCEDURE Real* (m: REAL; e: INTEGER): REAL;
377 VAR s: SET;
378 BEGIN
379 IF (m = 0) THEN RETURN 0.0 END;
380 ASSERT(~IsNan(m) & (1 <= ABS(m)) & (ABS(m) < 2), 20);
381 IF e = MAX(INTEGER) THEN
382 SYSTEM.GET(SYSTEM.ADR(m) + 4, s);
383 SYSTEM.PUT(SYSTEM.ADR(m) + 4, s + {20..30});
384 RETURN m
385 ELSE
386 FLD(e); FLD(m); FSCALE; FSTPst1; RETURN TOP()
387 END
388 END Real;
390 BEGIN
391 eps := 1.0E+0; e := 2.0E+0;
392 WHILE e > 1.0E+0 DO eps := eps/2.0E+0; e := 1.0E+0 + eps END; eps := 2.0E+0 * eps;
393 END Math.
397 PROCEDURE Log* (x: REAL): REAL;
398 BEGIN
399 RETURN Ln(x)/ln10
400 END Log;
402 PROCEDURE Power* (x, y: REAL): REAL;
403 BEGIN
404 RETURN Exp(y * Ln(x))
405 END Power;
407 PROCEDURE IntPower* (x: REAL; n: LONGINT): REAL;
408 VAR y: REAL;
409 BEGIN y := 1.0E+0;
410 IF n < 0 THEN x := 1.0E+0/x; n := -n END;
411 WHILE n > 0 DO
412 IF ODD(n) THEN y := y*x; DEC(n)
413 ELSE x := x * x; n := n DIV 2
414 END
415 END;
416 RETURN y
417 END IntPower;
419 PROCEDURE Tan* (x: REAL): REAL;
420 BEGIN
421 RETURN Sin(x)/Cos(x)
422 END Tan;
424 PROCEDURE ArcSin* (x: REAL): REAL;
425 BEGIN
426 RETURN 2.0E+0 * ArcTan(x/(1.0E+0 + Sqrt(1.0E+0 - x*x)))
427 END ArcSin;
429 PROCEDURE ArcCos* (x: REAL): REAL;
430 BEGIN (* pi/2 - arcsin(x) *)
431 RETURN Pi()/2.0E+0 - 2.0E+0 * ArcTan(x/(1.0E+0 + Sqrt(1.0E+0 - x*x)))
432 (*
433 IF x = -1 THEN RETURN Pi()
434 ELSE RETURN 2 * ArcTan(Sqrt((1 - x) / (1 + x)))
435 END
436 *) END ArcCos;
438 PROCEDURE ArcTan2* (y, x: REAL): REAL;
439 BEGIN
440 IF x = 0.0 THEN
441 RETURN Sign(y) * Pi() / 2.0
442 ELSIF y = 0.0 THEN
443 RETURN (1.0 - Sign(x)) * Pi() / 2.0
444 ELSE
445 RETURN ArcTan(y/x) + (1.0 - Sign(x)) * Sign(y) * Pi() / 2.0
446 END
447 END ArcTan2;
449 PROCEDURE Sinh* (x: REAL): REAL;
450 BEGIN
451 IF ABS(x) < -lneps THEN RETURN (Exp(x)-Exp(-x))/2.0E+0
452 ELSE RETURN Sign(x)*Exp(ABS(x))/2.0E+0
453 END
454 END Sinh;
456 PROCEDURE Cosh* (x: REAL): REAL;
457 BEGIN
458 IF ABS(x) < -lneps THEN RETURN (Exp(x)+Exp(-x))/2.0E+0
459 ELSE RETURN Exp(ABS(x))/2.0E+0
460 END
461 END Cosh;
463 PROCEDURE Tanh* (x: REAL): REAL;
464 VAR e1, e2: REAL;
465 BEGIN
466 IF ABS(x) < -lneps THEN
467 e1 := Exp(x); e2 := 1.0E+0/e1;
468 RETURN (e1-e2)/(e1+e2)
469 ELSE
470 RETURN Sign(x)
471 END
472 END Tanh;
474 PROCEDURE ArcSinh* (x: REAL): REAL;
475 BEGIN
476 IF x >= 0.0E+0 THEN RETURN Ln(x + Sqrt(x*x + 1.0E+0))
477 ELSE RETURN - Ln(-x + Sqrt(x*x + 1.0E+0))
478 END
479 END ArcSinh;
481 PROCEDURE ArcCosh* (x: REAL): REAL;
482 BEGIN
483 RETURN Ln(x + Sqrt(x*x - 1.0E+0))
484 END ArcCosh;
486 PROCEDURE ArcTanh* (x: REAL): REAL;
487 BEGIN
488 RETURN Ln((1.0E+0 + x)/(1.0E+0 - x))/2.0E+0
489 (* Variants:
490 (Ln(1+x)-Ln(1-x))/2.0E+0
491 -Ln((1-x)/Sqrt(1-x*x))
492 arcsinh(x/sqrt(1-x*x))
493 *)
494 END ArcTanh;
496 PROCEDURE Floor* (x: REAL): REAL;
497 BEGIN
498 IF ABS(x) >= 1.0E16 THEN RETURN x
499 ELSE RETURN ENTIER(x)
500 END
501 END Floor;
503 PROCEDURE Ceiling* (x: REAL): REAL;
504 BEGIN
505 IF ABS(x) >= 1.0E16 THEN RETURN x
506 ELSE RETURN -ENTIER(-x)
507 END
508 END Ceiling;
510 PROCEDURE Round* (x: REAL): REAL;
511 BEGIN
512 IF ABS(x) >= 1.0E16 THEN RETURN x
513 ELSE RETURN ENTIER(x + 0.5)
514 END
515 END Round;
517 PROCEDURE Trunc* (x: REAL): REAL;
518 BEGIN
519 IF ABS(x) >= 1.0E16 THEN RETURN x
520 ELSIF x >= 0 THEN RETURN ENTIER(x)
521 ELSE RETURN -ENTIER(-x)
522 END
523 END Trunc;
525 PROCEDURE Frac* (x: REAL): REAL;
526 BEGIN
527 IF ABS(x) >= 1.0E16 THEN RETURN 0.0
528 ELSIF x >= 0 THEN RETURN x - ENTIER(x)
529 ELSE RETURN x + ENTIER(-x)
530 END
531 END Frac;