(*******************************************************************************
Name:		MathLib.mod
Version:	Amiga 1.02
Purpose:	Standard MathLib
Created:	27.10.86  18:00  jr
Modified:	21.07.87  cn/red Version 3.0
Modified:	12.02.89  ms Konstantendeklarationen für die Fehlermeldungen
Modified:	27.03.89  ms Deutsche Version
Modified:	08.06.90  bp Resident Version
Modified:	02.01.91  bp Ausnutzung der FPU und MathIEESingTrans
*******************************************************************************)
IMPLEMENTATION MODULE MathLib0;
(*$
   StackChk:=FALSE
   RangeChk:=FALSE
   OverflowChk:=FALSE
   LargeVars:=FALSE
   Volatile:=FALSE
   StackParms:=FALSE
   EntryClear:=FALSE
*)

FROM SYSTEM IMPORT
 ADR, ASSEMBLE, BITSET, CAST;
IMPORT Arts;

(*$ IF m68881 AND NOT m68040 *) (* 68040 muß selbst sehen.. *)
IMPORT FPUExc;

(*$ EntryExitCode:=FALSE *)
PROCEDURE sin	(x: REAL): REAL;
BEGIN
  ASSEMBLE(
	FSIN.S	x-4(A7),FP0 (* Parallelität ausnutzen! *)
	MOVE.L	(A7)+,A0  (* Dies kostet GAR NICHTS! *)
	ADDQ.L	#4,A7
	FMOVE.S	FP0,D0
	JMP	(A0)
  END);
END sin;

(*$ EntryExitCode:=FALSE *)
PROCEDURE cos	(x: REAL): REAL;
BEGIN
  ASSEMBLE(
	FCOS.S	x-4(A7),FP0 (* Parallelität ausnutzen! *)
	MOVE.L	(A7)+,A0
	ADDQ.L	#4,A7
	FMOVE.S	FP0,D0
	JMP	(A0)
  END);
END cos;

(*$ EntryExitCode:=FALSE *)
PROCEDURE arctan(x: REAL): REAL;
BEGIN
  ASSEMBLE(
	FATAN.S	x-4(A7),FP0 (* Parallelität ausnutzen! *)
	MOVE.L	(A7)+,A0
	ADDQ.L	#4,A7
	FMOVE.S	FP0,D0
	JMP	(A0)
  END);
END arctan;

(*$ EntryExitCode:=FALSE *)
PROCEDURE exp	(x: REAL): REAL;
BEGIN
  ASSEMBLE(
	FETOX.S	x-4(A7),FP0 (* Parallelität ausnutzen! *)
	MOVE.L	(A7)+,A0
	ADDQ.L	#4,A7
	FMOVE.S	FP0,D0
	JMP	(A0)
  END);
END exp;

(*$ EntryExitCode:=FALSE *)
PROCEDURE ln	(x: REAL): REAL;
BEGIN
  ASSEMBLE(
	FLOGN.S	x-4(A7),FP0 (* Parallelität ausnutzen! *)
	MOVE.L	(A7)+,A0
	ADDQ.L	#4,A7
	FMOVE.S	FP0,D0
	JMP	(A0)
  END);
END ln;

(*$ EntryExitCode:=FALSE *)
PROCEDURE sqrt	(x: REAL): REAL;
BEGIN
  ASSEMBLE(
	FSQRT.S	x-4(A7),FP0 (* Parallelität ausnutzen! *)
	MOVE.L	(A7)+,A0
	ADDQ.L	#4,A7
	FMOVE.S	FP0,D0
	JMP	(A0)
  END);
END sqrt;


(*$ ELSE *) (* normaler Prozessor *)
IMPORT o:OptMathIEEESingTransL;

TYPE IEEE=RECORD CASE :INTEGER OF
          | 0:r:REAL;
          | 1:hi,lo:BITSET;
          END END;

CONST
 oneDIVsqrt2= 0.7071067814;
 piDIV2     = 1.570796327;
 ln2H       = 0.693359;
 ln2L       = -2.118194E-4;
 bigArg     =
      "Argument zu gross";
(*    "Argument too big"; *)
 mathErr    =
     "Fehler in MathLib0";
(*   "Error in MathLib0"; *)
 lnErr="ln(x): x<=0.0";
 sqrtErr="sqrt(x): x<0.0";

PROCEDURE round(x:REAL):INTEGER;
 BEGIN
  IF x<0.0 THEN RETURN TRUNC(x-0.5)
  ELSE RETURN TRUNC(x+0.5)
  END
 END round;

PROCEDURE expIn(VAR x:REAL; exp:INTEGER);
 VAR t:IEEE;
 BEGIN
  IF exp>255 THEN Arts.Error(ADR(mathErr),ADR(bigArg)); END;
  WITH t DO
   r:=x;
   hi:=hi*{0..6,15}+CAST(BITSET,exp*128);
   x:=r
  END
 END expIn;

PROCEDURE expOut(x:REAL):INTEGER;
 VAR t:IEEE;
 BEGIN
  WITH t DO
   r:=x;
   RETURN CAST(INTEGER,hi*{7..14}) DIV 128;
  END
 END expOut;

PROCEDURE sico(x:REAL; sgn:INTEGER):REAL;
 (*
               9    7    5    3
     sin x = Ax + Bx + Cx + Dx + x  , x in [-pi/2..pi/2]
 *)
 CONST
  oneDIVpi=0.31830989;
  A=2.601903E-6; B=-1.980742E-4; C=8.333025E-3; D=-0.166667;

 VAR
  xQuad: REAL;
  f: INTEGER;

 BEGIN
  IF x>1.0E6 THEN Arts.Error(ADR(mathErr),ADR(bigArg)); END;
  f:=round(x*oneDIVpi);
  IF ODD(f) THEN sgn:=-sgn END;
  x:=x-(FLOAT(f)*3.140625)-(FLOAT(f)*9.676537E-4);
  IF ABS(x)>=2.5E-4 THEN
   xQuad:=x*x;
   x:=x*(((A*xQuad + B)*xQuad + C)*xQuad + D)*xQuad + x
  END;
  RETURN x*FLOAT(sgn)
 END sico;

PROCEDURE sin(x:REAL):REAL;
BEGIN
  IF o.mathieeesingtransBase#NIL THEN
    RETURN o.Sin(x)
  ELSE
    IF x<0.0 THEN RETURN sico(-x,-1)
    ELSE RETURN sico(x,1)
    END
  END;
END sin;

PROCEDURE cos(x:REAL):REAL;
BEGIN
  IF o.mathieeesingtransBase#NIL THEN
    RETURN o.Cos(x)
  ELSE
    RETURN sico(ABS(x) + piDIV2,1)
  END;
END cos;

PROCEDURE exp(x:REAL):REAL;
 (*
      x    y+f*ln2/2    y  f+1
     e  = e          = e *2     , y in [-ln2/2..ln2/2]
                3
      y       Ax + Cx
     e  = ------------------ + E
             3    2
          -Ax + Bx - Cx + D
 *)

CONST
   oneDIVln2=1.442695041;
   A=4.160289E-3; B=4.998718E-2; C=0.25; D=0.5; E=0.5;

VAR
  x1, xQuad: REAL;
  f: INTEGER;

BEGIN
  IF x=0.0 THEN RETURN 1.0
  ELSIF o.mathieeesingtransBase#NIL THEN
    RETURN o.Exp(x)
  ELSE
    f:=round(x*oneDIVln2);
    x:=x-FLOAT(f)*ln2H-FLOAT(f)*ln2L;
    xQuad:=x*x;
    x1:=x*(A*xQuad + C); x1:=E + x1/(B*xQuad + D - x1);
    expIn(x1,expOut(x1)+1+f);
    RETURN x1
  END
END exp;

PROCEDURE arctan(x:REAL):REAL;
 (*
                      3      5      7
      arctan x = x - x /3 + x /5 - x /7 +-...
                  mit Hilfe einer Doppelbruchentwicklung auf 8 Glieder.
 *)

VAR
  x1, fac: REAL;
  i: INTEGER;
BEGIN
  IF o.mathieeesingtransBase#NIL THEN
    RETURN o.Atan(x)
  ELSE
    x1:=x*x;
    IF ABS(x)>1.0 THEN x1:=1.0/x1 END;
    fac:=17.0;
    FOR i:=8 TO 1 BY -1 DO
      fac:=FLOAT(2*i-1) + FLOAT(i*i)*x1/fac
    END;
    IF ABS(x)>1.0 THEN
      x1:=piDIV2 - 1.0/(ABS(x)*fac);
      IF x<0.0 THEN RETURN -x1 END;
      RETURN x1
    END;
    RETURN x/fac
  END;
END arctan;

PROCEDURE ln(x:REAL):REAL;
 (*
                  e
      ln(x) = ln(2 *m) = e*ln2 + ln(m) , m in [1/2..1]

                                   3
                 1 + t          A*t
      ln(m) = ln(-----) = t + --------   , t in [-0.34..0.34]
                 1 - t              2
                               B - t
 *)

 CONST
  A=0.552707; B=6.632718;
 VAR
  temp, quad, rexp: REAL;
  exp: INTEGER;

BEGIN
  IF x<=0.0 THEN Arts.Error(ADR(mathErr),ADR(lnErr)); END;
  IF o.mathieeesingtransBase#NIL THEN
    RETURN o.Log(x)
  ELSE
    exp:=expOut(x)-126; expIn(x,126);
    IF x>oneDIVsqrt2 THEN temp:=(x - 1.0)/(x*0.5 + 0.5)
    ELSE
      DEC(exp);
      temp:=x-0.5; temp:=temp/(temp*0.5 + 0.5)
    END;
    quad:=temp*temp; rexp:=REAL(exp);
    RETURN A*quad/(B-quad)*temp+temp+ln2H*rexp+ln2L*rexp;
  END;
END ln;

PROCEDURE sqrt(x:REAL):REAL;
 (*
              e     e/2
      sqrt(m*2 ) = 2   * sqrt(m) , m in [1/2..1]

      sqrt(m) berechnet durch lineare Naeherung und 3 Newtonschritte
 *)

VAR
  root: REAL;
  exp, i: INTEGER;
BEGIN
  IF x<0.0 THEN Arts.Error(ADR(mathErr),ADR(sqrtErr)); END;
  IF x=0.0 THEN RETURN 0.0
  ELSIF o.mathieeesingtransBase#NIL THEN
    RETURN o.Sqrt(x)
  ELSE
    exp:=expOut(x)-126; expIn(x,126);
    root:=0.41731 + 0.59016*x;                      (* Naeherung +    *)
    FOR i:=0 TO 2 DO root:=0.5*(root + x/root) END; (* Newtonschritte *)
    IF ODD(exp) THEN INC(exp); root:=root*oneDIVsqrt2 END;
    expIn(root,expOut(root) + exp DIV 2);
    RETURN root
  END;
END sqrt;
(*$ ENDIF *)

END MathLib0.