Astrobe Forum

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.

Astrobe Forum

Some modifications to the MATH module

Some modifications to the MATH module