MODULE SMath;
(**
project = "BlackBox"
organization = "www.oberon.ch"
contributors = "Oberon microsystems"
version = "System/Rsrc/About"
copyright = "System/Rsrc/About"
license = "Docu/BB-License"
changes = ""
issues = ""
**)
IMPORT SYSTEM;
VAR eps, e: SHORTREAL;
(* code procedures for 80387 math coprocessor *)
PROCEDURE [code] FLD (x: SHORTREAL);
PROCEDURE [code] TOP (): SHORTREAL;
PROCEDURE [code] FSW (): INTEGER 0DFH, 0E0H;
PROCEDURE [code] FSWs (): SET 0DFH, 0E0H;
PROCEDURE [code] ST0 (): SHORTREAL 0D9H, 0C0H;
PROCEDURE [code] ST1 (): SHORTREAL 0D9H, 0C1H;
PROCEDURE [code] FXCH 0D9H, 0C9H;
PROCEDURE [code] FLDst0 0D9H, 0C0H; (* doublicate st[0] *)
PROCEDURE [code] FSTPst0 0DDH, 0D8H; (* remove st[0] *)
PROCEDURE [code] FSTPst1 0DDH, 0D9H; (* remove st[1] *)
PROCEDURE [code] FSTPDe 0DBH, 05DH, 0F4H; (* FSTPD -12[FP] *) (* COMPILER DEPENDENT *)
PROCEDURE [code] WAIT 09BH;
PROCEDURE [code] FNOP 0D9H, 0D0H;
PROCEDURE [code] FLD0 0D9H, 0EEH;
PROCEDURE [code] FLD1 0D9H, 0E8H;
PROCEDURE [code] FLDPI 0D9H, 0EBH;
PROCEDURE [code] FLDLN2 0D9H, 0EDH;
PROCEDURE [code] FLDLG2 0D9H, 0ECH;
PROCEDURE [code] FLDL2E 0D9H, 0EAH;
PROCEDURE [code] FADD 0DEH, 0C1H;
PROCEDURE [code] FADDst0 0D8H, 0C0H;
PROCEDURE [code] FSUB 0DEH, 0E9H;
PROCEDURE [code] FSUBn 0DCH, 0E9H; (* no pop *)
PROCEDURE [code] FSUBR 0DEH, 0E1H;
PROCEDURE [code] FSUBst1 0D8H, 0E1H;
PROCEDURE [code] FMUL 0DEH, 0C9H;
PROCEDURE [code] FMULst0 0D8H, 0C8H;
PROCEDURE [code] FMULst1st0 0DCH, 0C9H;
PROCEDURE [code] FDIV 0DEH, 0F9H;
PROCEDURE [code] FDIVR 0DEH, 0F1H;
PROCEDURE [code] FDIVRst1 0D8H, 0F9H;
PROCEDURE [code] FCHS 0D9H, 0E0H;
PROCEDURE [code] FCOM 0D8H, 0D1H;
PROCEDURE [code] FSWax 0DFH, 0E0H;
PROCEDURE [code] SAHF 09EH;
PROCEDURE [code] JBE4 076H, 004H;
PROCEDURE [code] JAE4 073H, 004H;
PROCEDURE [code] FRNDINT 0D9H, 0FCH;
PROCEDURE [code] FSCALE 0D9H, 0FDH; (* st[0] * 2^FLOOR(st[1]) *)
PROCEDURE [code] FXTRACT 0D9H, 0F4H; (* exp -> st[1]; mant -> st[0] *)
PROCEDURE [code] FXAM 0D9H, 0E5H;
PROCEDURE [code] FSQRT 0D9H, 0FAH; (* st[0] >= 0 *)
PROCEDURE [code] FSIN 0D9H, 0FEH; (* |st[0]| < 2^63 *)
PROCEDURE [code] FCOS 0D9H, 0FFH; (* |st[0]| < 2^63 *)
PROCEDURE [code] FTAN 0D9H, 0F2H; (* |st[0]| < 2^63 *)
PROCEDURE [code] FATAN 0D9H, 0F3H; (* atan2(st[1], st[0]) *)
PROCEDURE [code] FYL2X 0D9H, 0F1H; (* st[1] * log2(st[0]), st[0] > 0 *)
PROCEDURE [code] FYL2XP1 0D9H, 0F9H; (* st[1] * log2(1 + st[0]), |st[0]| < 1-sqrt(2)/2 *)
PROCEDURE [code] F2XM1 0D9H, 0F0H; (* 2^st[0] - 1, |st[0]| <= 1 *)
PROCEDURE IsNan (x: SHORTREAL): BOOLEAN;
BEGIN
FLD(x); FXAM; FSTPst0; WAIT; RETURN FSWs() * {8, 10} = {8}
END IsNan;
(* sin, cos, tan argument reduction *)
PROCEDURE Reduce;
BEGIN
FXAM; WAIT;
IF ~(8 IN FSWs()) & (ABS(ST0()) > 1.0E18) THEN
(* to be completed *)
FSTPst0; FLD0
END;
END Reduce;
(** SHORTREAL precision **)
PROCEDURE Pi* (): SHORTREAL;
BEGIN
FLDPI; RETURN TOP()
END Pi;
PROCEDURE Eps* (): SHORTREAL;
BEGIN
RETURN eps
END Eps;
PROCEDURE Sqrt* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, argument of Sqrt must not be negative *)
FLD(x); FSQRT; WAIT; RETURN TOP()
END Sqrt;
PROCEDURE Exp* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 2 ^ (x * 1/ln(2)) *)
FLD(x); FLDL2E; FMUL;
IF ABS(ST0()) = INF THEN FLD1
ELSE FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD
END;
FSCALE; FSTPst1; RETURN TOP()
END Exp;
PROCEDURE Ln* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, argument of Ln must not be negative *)
(* ln(2) * ld(x) *)
FLDLN2; FLD(x); FYL2X; WAIT; RETURN TOP()
END Ln;
PROCEDURE Log* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, argument of Log must not be negative *)
(* log(2) * ld(x) *)
FLDLG2; FLD(x); FYL2X; WAIT; RETURN TOP()
END Log;
PROCEDURE Power* (x, y: SHORTREAL): SHORTREAL;
BEGIN
ASSERT(x >= 0, 20);
ASSERT((x # 0.0)OR(y # 0.0), 21);
ASSERT((x # INF)OR(y # 0.0), 22);
ASSERT((x # 1.0)OR(ABS(y) # INF), 23);
(* 2 ^ (y * ld(x)) *)
FLD(y); FLD(x); FYL2X;
IF ABS(ST0()) = INF THEN FLD1
ELSE FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD
END;
FSCALE; FSTPst1; WAIT; RETURN TOP()
END Power;
PROCEDURE IntPower* (x: SHORTREAL; n: INTEGER): SHORTREAL;
BEGIN
FLD1; FLD(x);
IF n = MIN(INTEGER) THEN RETURN IntPower(x, n + 1) / x END;
IF n <= 0 THEN FDIVRst1; (* 1 / x *) n := -n END;
WHILE n > 0 DO
IF ODD(n) THEN FMULst1st0; (* y := y * x *) DEC(n)
ELSE FMULst0; (* x := x * x *) n := n DIV 2
END
END;
FSTPst0; RETURN TOP()
END IntPower;
PROCEDURE Sin* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, ABS(x) # INF *)
FLD(x); Reduce; FSIN; WAIT; RETURN TOP()
END Sin;
PROCEDURE Cos* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, ABS(x) # INF *)
FLD(x); Reduce; FCOS; WAIT; RETURN TOP()
END Cos;
PROCEDURE Tan* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, ABS(x) # INF *)
FLD(x); Reduce; FTAN; FSTPst0; WAIT; RETURN TOP()
END Tan;
PROCEDURE ArcSin* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, -1.0 <= x <= 1.0 *)
(* atan2(x, sqrt(1 - x*x)) *)
FLD(x); FLDst0; FMULst0; FLD1; FSUBR; FSQRT; FNOP; FATAN; WAIT; RETURN TOP()
END ArcSin;
PROCEDURE ArcCos* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, -1.0 <= x <= 1.0 *)
(* atan2(sqrt(1 - x*x), x) *)
FLD(x); FMULst0; FLD1; FSUBR; FSQRT; FLD(x); FATAN; WAIT; RETURN TOP()
END ArcCos;
PROCEDURE ArcTan* (x: SHORTREAL): SHORTREAL;
BEGIN
(* atan2(x, 1) *)
FLD(x); FLD1; FATAN; RETURN TOP()
END ArcTan;
PROCEDURE ArcTan2* (y, x: SHORTREAL): SHORTREAL;
BEGIN
ASSERT((y # 0)OR (x # 0), 20);
ASSERT((ABS(y) # INF)OR(ABS(x)# INF), 21);
FLD(y); FLD(x); FATAN; WAIT; RETURN TOP()
END ArcTan2;
PROCEDURE Sinh* (x: SHORTREAL): SHORTREAL;
BEGIN
(* IF IsNan(x) THEN RETURN x END; *)
(* abs(x) * 1/ln(2) *)
FLD(ABS(x)); FLDL2E; FMUL;
IF ST0() < 0.5 THEN
(* (2^z - 1) + (2^z - 1) / ((2^z - 1) + 1) *)
F2XM1; FLDst0; FLDst0; FLD1; FADD; FDIV; FADD
ELSIF ST0() # INF THEN
(* 2^z - 1 / 2^z *)
FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD; FSCALE;
FSTPst1; FLDst0; FLD1; FDIVR; FSUB
END;
IF x < 0 THEN FCHS END;
RETURN TOP() * 0.5
END Sinh;
PROCEDURE Cosh* (x: SHORTREAL): SHORTREAL;
BEGIN
(* IF IsNan(x) THEN RETURN x END; *)
(* 2^(abs(x) * 1/ln(2)) *)
FLD(ABS(x));
IF ST0() # INF THEN
FLDL2E; FMUL; FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD; FSCALE;
FSTPst1;
(* z + 1/z *)
FLDst0; FLD1; FDIVR; FADD
END;
RETURN TOP() * 0.5
END Cosh;
PROCEDURE Tanh* (x: SHORTREAL): SHORTREAL;
BEGIN
(* IF IsNan(x) THEN RETURN x END; *)
(* abs(x) * 1/ln(2) * 2 *)
FLD(ABS(x)); FLDL2E; FMUL; FADDst0;
IF ST0() < 0.5 THEN
(* (2^z - 1) / (2^z + 1) *)
F2XM1; FLDst0; FLD(2); FADD; FDIV
ELSIF ST0() < 65 THEN
(* 1 - 2 / (2^z + 1) *)
FLDst0; FRNDINT; FXCH; FSUBst1; FNOP; F2XM1; FLD1; FADD; FSCALE; FSTPst1; FLD1; FADD; FLD(2); FDIVR; FLD1; FSUBR
ELSE
FSTPst0; FLD1
END;
IF x < 0 THEN FCHS END;
RETURN TOP()
END Tanh;
PROCEDURE ArcSinh* (x: SHORTREAL): SHORTREAL;
BEGIN
(* IF IsNan(x) THEN RETURN x END; *)
(* x*x *)
FLDLN2; FLD(ABS(x)); FLDst0; FMULst0;
IF ST0() < 0.067 THEN
(* ln(2) * ld(1 + x*x / (sqrt(x*x + 1) + 1) + x) *)
FLDst0; FLD1; FADD; FSQRT; FLD1; FADD; FDIV; FADD; FYL2XP1
ELSE
(* ln(2) * ld(x + sqrt(x*x + 1)) *)
FLD1; FADD; FSQRT; FADD; FYL2X
END;
IF x < 0 THEN FCHS END;
RETURN TOP()
END ArcSinh;
PROCEDURE ArcCosh* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, x >= 1.0 *)
(* IF IsNan(x) THEN RETURN x END; *)
(* ln(2) * ld(x + sqrt(x*x - 1)) *)
FLDLN2; FLD(x); FLDst0; FMULst0; FLD1; FSUB; FSQRT; FADD; FYL2X; WAIT; RETURN TOP()
END ArcCosh;
PROCEDURE ArcTanh* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, -1.0 <= x <= 1.0 *)
(* IF IsNan(x) THEN RETURN x END; *)
(* |x| *)
FLDLN2; FLD(ABS(x));
IF ST0() < 0.12 THEN
(* ln(2) * ld(1 + 2*x / (1 - x)) *)
FLDst0; FLD1; FSUBR; FDIV; FADDst0; FYL2XP1
ELSE
(* ln(2) * ld((1 + x) / (1 - x)) *)
FLDst0; FLD1; FADD; FXCH; FLD1; FSUBR; FDIV; FNOP; FYL2X
END;
IF x < 0 THEN FCHS END;
WAIT;
RETURN TOP() * 0.5
END ArcTanh;
PROCEDURE Floor* (x: SHORTREAL): SHORTREAL;
BEGIN
FLD(x); FLDst0; FRNDINT; FCOM; FSWax; FSTPst1; SAHF; JBE4; FLD1; FSUB; RETURN TOP()
END Floor;
PROCEDURE Ceiling* (x: SHORTREAL): SHORTREAL;
BEGIN
FLD(x); FLDst0; FRNDINT; FCOM; FSWax; FSTPst1; SAHF; JAE4; FLD1; FADD; RETURN TOP()
END Ceiling;
PROCEDURE Round* (x: SHORTREAL): SHORTREAL;
BEGIN
FLD(x);
IF ABS(ST0()) = INF THEN RETURN TOP() END;
FLDst0; FRNDINT; FSUBn; FXCH;
IF TOP() = 0.5 THEN FLD1; FADD END;
RETURN TOP()
END Round;
PROCEDURE Trunc* (x: SHORTREAL): SHORTREAL;
BEGIN
FLD(x); FLDst0; FRNDINT;
IF ST1() >= 0 THEN
FCOM; FSWax; FSTPst1; SAHF; JBE4; FLD1; FSUB
ELSE
FCOM; FSWax; FSTPst1; SAHF; JAE4; FLD1; FADD
END;
RETURN TOP()
END Trunc;
PROCEDURE Frac* (x: SHORTREAL): SHORTREAL;
BEGIN
(* 20, x # INF&x # -INF *)
FLD(x); FLDst0; FRNDINT;
IF ST1() >= 0 THEN
FCOM; FSWax; SAHF; JBE4; FLD1; FSUB
ELSE
FCOM; FSWax; SAHF; JAE4; FLD1; FADD
END;
FSUB; WAIT; RETURN TOP()
END Frac;
PROCEDURE Sign* (x: SHORTREAL): SHORTREAL;
BEGIN
FLD(x); FXAM; WAIT;
CASE FSW() DIV 256 MOD 8 OF
| 0, 2: FSTPst0; RETURN 0.0
| 1, 4, 5: FSTPst0; RETURN 1.0
| 3, 6, 7: FSTPst0; RETURN -1.0
END
END Sign;
PROCEDURE Mantissa* (x: SHORTREAL): SHORTREAL;
BEGIN
FLD(x); FXAM; WAIT;
CASE FSW() DIV 256 MOD 8 OF
| 4, 6: FXTRACT; FSTPst1; RETURN TOP()
| 0, 2: FSTPst0; RETURN 0.0 (* zero *)
| 5: FSTPst0; RETURN 1.0 (* inf *)
| 7: FSTPst0; RETURN -1.0 (* -inf *)
| 1: FSTPst0; RETURN 1.5 (* nan *)
| 3: FSTPst0; RETURN -1.5 (* -nan *)
END
END Mantissa;
PROCEDURE Exponent* (x: SHORTREAL): INTEGER; (* COMPILER DEPENDENT *)
VAR e: INTEGER; (* e is set by FSTPDe! *)
BEGIN
FLD(x); FXAM; WAIT;
CASE FSW() DIV 256 MOD 8 OF
| 4, 6: FXTRACT; FSTPst0; FSTPDe; WAIT; RETURN e
| 0, 2: FSTPst0; RETURN 0 (* zero *)
| 1, 3, 5, 7: FSTPst0; RETURN MAX(INTEGER) (* inf or nan*)
END
END Exponent;
PROCEDURE Real* (m: SHORTREAL; e: INTEGER): SHORTREAL;
VAR s: SET;
BEGIN
IF (m = 0) THEN RETURN 0.0 END;
ASSERT(~IsNan(m) & (1 <= ABS(m)) & (ABS(m) < 2), 20);
IF e = MAX(INTEGER) THEN
SYSTEM.GET(SYSTEM.ADR(m) + 4, s);
SYSTEM.PUT(SYSTEM.ADR(m) + 4, s + {20..30});
RETURN m
ELSE
FLD(e); FLD(m); FSCALE; FSTPst1; RETURN TOP()
END
END Real;
BEGIN
eps := 1.0E+0; e := 2.0E+0;
WHILE e > 1.0E+0 DO eps := eps/2.0E+0; e := 1.0E+0 + eps END; eps := 2.0E+0 * eps;
END SMath.