Page 1 of 1

### Some modifications to the MATH module

Posted: Tue Mar 30, 2021 5:57 pm
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.