Some modifications to the MATH module

[bscholte]

While extensively using the MATH module, we made a few modifications that made the code a bit faster and more robust. Maybe these are useful.....

1.) Avoiding unnecessary iterations during Exp
The following uses some simple big steps to get a factor close to 1 and then falls back to the original ExpSeries:

Code: Select all

PROCEDURE Exp*(x: REAL): REAL;
CONST exp10 = 22026.46579;
      exp5  = 148.4131591; 
      exp1  = 2.718282;
VAR Result: REAL;
BEGIN 
  IF x=0.0 THEN Result:=1.0
  ELSIF x<0.0 THEN IF    x< -30.0 THEN Result:=1.0E-13
                   ELSIF x< -10.0 THEN Result:=Exp(x+10.0) / exp10
                   ELSIF x<  -5.0 THEN Result:=Exp(x+ 5.0) / exp5
                   ELSIF x<  -1.0 THEN Result:=Exp(x+ 1.0) / exp1
                   ELSE Result:=ExpSeries(x)
                   END
  ELSE (*x>0.0*)   IF    x>  30.0 THEN Result:=1.0E+13 
                   ELSIF x>  10.0 THEN Result:=Exp(x-10.0) * exp10
                   ELSIF x>   5.0 THEN Result:=Exp(x- 5.0) * exp5
                   ELSIF x>   1.0 THEN Result:=Exp(x- 1.0) * exp1
                   ELSE Result:=ExpSeries(x)
                   END;
  END;
  RETURN Result
END Exp;

2.) We found a rather large deviation of the Ln when going to small numbers. We added a different approximation, next to the original LnSeries, for values less than 0.5:

Code: Select all

CONST ln2*       = 0.69314718;
      MaxReal*      = 3.40282E+38;
      SmallestReal* = 1.17549E-38;
      MaxInt*       = 7FFFFH;

(*
  Newton's method for better log-approximation for x < 0.5
*)
PROCEDURE LogNewton(x: REAL): REAL;
VAR
  yn, Result: REAL;
BEGIN
  yn:=x;
  Result:=x-1.0;
  WHILE ABS(yn-Result)>0.0001 DO 
    yn:=Result;
    Result:=yn+2.0*(x-Exp(yn))/(x+Exp(x));
  END;
  RETURN Result
END LogNewton;


PROCEDURE Ln*(x: REAL): REAL;
VAR
  e: INTEGER;
  Result:REAL;
BEGIN
  (*ASSERT(x > 0.0, 22);*)
  (*Out.String("ln("); Out.Real(x); Out.Ln;*)
  IF     x=1.0  THEN Result:=0.0 (*definition*)
  ELSIF (x>0.0) THEN
	IF x<0.5 THEN
		Result:=LogNewton(x);
	ELSE 
		UNPK(x, e);
		PACK(x, 0);
        IF x > 1.5 THEN x := x * 0.5; INC(e, 1) END;
        Result:=(FLT(e) * ln2) + LnSeries(x)
	END;
  ELSE Out.String("ln("); Out.Real(x); Out.Ln; Result := -1.0 * MaxReal (*zero or smaller returns minimum real value to avoid runtime error*)
  END;
  RETURN Result
END Ln;

3.) to avoid overrun and underrun errors while dividing and multiplying Reals, we are using two seperate functions. In case of overruns the procedures return the largest possible values (+MaxReal or -MaxReal). In case ove underrruns, the smallest possible values (+ or - SmallestReal). This can probably be done more efficient by not using EXP10 exponents, but this seems to work fine:

Code: Select all

CONST MaxReal*      = 3.40282E+38;
      SmallestReal* = 1.17549E-38;

PROCEDURE rDIV*(A,B:REAL):REAL;
(*Real division A/B with overflow protection*)
VAR Result:REAL;
BEGIN
  IF    ABS(B) >= 1.0                 THEN Result:= A / B      (*Since |A|<=MaxReal, any |B|>=1.0 will always yield a result*)                              
  ELSIF (ABS(A) < (ABS(B) * MaxReal)) THEN Result:= A / B      (*Checking |A/B| < MaxReal with |A|<|B|*MaxReal: Since SmallestReal>|B|>1.0, B * MaxReal will always yield result*)
  ELSIF (Reals.Exp10(A)-Reals.Exp10(B)< -38) THEN IF BFX(A, 31) = BFX(B, 31) THEN Result:= SmallestReal ELSE Result:= -SmallestReal END;(*|A/B| would be <SmallestReal*) 
  ELSE IF BFX(A, 31) = BFX(B, 31) THEN Result:= MaxReal; psSysRecord.MathDivError ELSE Result:= -MaxReal END; (*Division by zero or overflow: just return MaxReal or MinReal*)
  END;
  RETURN Result
END rDIV;

PROCEDURE rMUL*(A,B:REAL):REAL;
(*Real Multiplication with overflow protection*)
VAR Result:REAL;
    ExpSum:INTEGER;
BEGIN
  ExpSum:=Reals.Exp10(A) + Reals.Exp10(B);
  IF    ABS(ExpSum)<=  37 THEN Result := A * B (*Result is ok. Cuts off results at 9.99E+/-37*)
  ELSIF     ExpSum >   37 THEN IF BFX(A, 31) = BFX(B, 31) THEN Result := MaxReal ELSE Result := -MaxReal END; (*|A*B|>MaxReal*)
  ELSIF     ExpSum <  -37 THEN IF BFX(A, 31) = BFX(B, 31) THEN Result := SmallestReal ELSE Result := - SmallestReal END; (*|A*B|<SmallestReal*) 
  END;     
  RETURN Result
END rMUL;

In our application, we actually also log such over and underruns in an error record to catch them so we can try to avoid it.