vishapoberon/compiler
1
0
Fork 0
mirror of https://github.com/vishapoberon/compiler.git synced 2026-04-05 22:12:24 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

ported four more modules from ooc lib

Browse source
This commit is contained in:
Norayr Chilingarian 2013-10-18 18:43:15 +04:00
parent 24b30cdae7
commit c6a299a1fe
10 changed files with 1539 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

2
makefile
View file

@ -134,6 +134,8 @@ stage6:
$(VOCSTATIC) -sP ooc2Real0.Mod
$(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod
$(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.Mod
$(VOCSTATIC) -sP oocLRealMath.Mod
$(VOCSTATIC) -sP oocLongInts.Mod oocLRealConv.Mod oocLRealStr.Mod
$(VOCSTATIC) -sP oocwrapperlibc.Mod
$(VOCSTATIC) -sP ulmSYSTEM.Mod
$(VOCSTATIC) -sP ulmASCII.Mod ulmSets.Mod

2
makefile.gnuc.armv6j
View file

@ -134,6 +134,8 @@ stage6:
$(VOCSTATIC) -sP ooc2Real0.Mod
$(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod
$(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.Mod
$(VOCSTATIC) -sP oocLRealMath.Mod
$(VOCSTATIC) -sP oocLongInts.Mod oocLRealConv.Mod oocLRealStr.Mod
$(VOCSTATIC) -sP oocwrapperlibc.Mod
$(VOCSTATIC) -sP ulmSYSTEM.Mod
$(VOCSTATIC) -sP ulmASCII.Mod ulmSets.Mod

2
makefile.gnuc.armv6j_hardfp
View file

@ -134,6 +134,8 @@ stage6:
$(VOCSTATIC) -sP ooc2Real0.Mod
$(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod
$(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.Mod
$(VOCSTATIC) -sP oocLRealMath.Mod
$(VOCSTATIC) -sP oocLongInts.Mod oocLRealConv.Mod oocLRealStr.Mod
$(VOCSTATIC) -sP oocwrapperlibc.Mod
$(VOCSTATIC) -sP ulmSYSTEM.Mod
$(VOCSTATIC) -sP ulmASCII.Mod ulmSets.Mod

2
makefile.gnuc.armv7a_hardfp
View file

@ -134,6 +134,8 @@ stage6:
$(VOCSTATIC) -sP ooc2Real0.Mod
$(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod
$(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.Mod
$(VOCSTATIC) -sP oocLRealMath.Mod
$(VOCSTATIC) -sP oocLongInts.Mod oocLRealConv.Mod oocLRealStr.Mod
$(VOCSTATIC) -sP oocwrapperlibc.Mod
$(VOCSTATIC) -sP ulmSYSTEM.Mod
$(VOCSTATIC) -sP ulmASCII.Mod ulmSets.Mod

2
makefile.gnuc.x86
View file

@ -134,6 +134,8 @@ stage6:
$(VOCSTATIC) -sP ooc2Real0.Mod
$(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod
$(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.Mod
$(VOCSTATIC) -sP oocLRealMath.Mod
$(VOCSTATIC) -sP oocLongInts.Mod oocLRealConv.Mod oocLRealStr.Mod
$(VOCSTATIC) -sP oocwrapperlibc.Mod
$(VOCSTATIC) -sP ulmSYSTEM.Mod
$(VOCSTATIC) -sP ulmASCII.Mod ulmSets.Mod

2
makefile.gnuc.x86_64
View file

@ -134,6 +134,8 @@ stage6:
$(VOCSTATIC) -sP ooc2Real0.Mod
$(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod
$(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.Mod
$(VOCSTATIC) -sP oocLRealMath.Mod
$(VOCSTATIC) -sP oocLongInts.Mod oocLRealConv.Mod oocLRealStr.Mod
$(VOCSTATIC) -sP oocwrapperlibc.Mod
$(VOCSTATIC) -sP ulmSYSTEM.Mod
$(VOCSTATIC) -sP ulmASCII.Mod ulmSets.Mod

414
src/lib/ooc/oocLRealConv.Mod Normal file
View file

@ -0,0 +1,414 @@
(* $Id: LRealConv.Mod,v 1.7 1999/10/12 07:17:54 ooc-devel Exp $ *)
MODULE oocLRealConv;
(*
LRealConv - Low-level LONGREAL/string conversions.
Copyright (C) 1996 Michael Griebling
This module is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.
This module is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*)
IMPORT
Char := oocCharClass, Low := oocLowLReal, Str := oocStrings, Conv := oocConvTypes,
LInt := oocLongInts, SYSTEM;
CONST
ZERO=0.0D0;
SigFigs*=15; (* accuracy of LONGREALs *)
DEBUG = FALSE;
TYPE
ConvResults*= Conv.ConvResults; (* strAllRight, strOutOfRange, strWrongFormat, strEmpty *)
LongInt=LInt.LongInt;
CONST
strAllRight*=Conv.strAllRight; (* the string format is correct for the corresponding conversion *)
strOutOfRange*=Conv.strOutOfRange; (* the string is well-formed but the value cannot be represented *)
strWrongFormat*=Conv.strWrongFormat; (* the string is in the wrong format for the conversion *)
strEmpty*=Conv.strEmpty; (* the given string is empty *)
VAR
RS, P, F, E, SE, WE, SR: Conv.ScanState;
PROCEDURE IsExponent (ch: CHAR) : BOOLEAN;
BEGIN
RETURN (ch="E") OR (ch="D")
END IsExponent;
PROCEDURE IsSign (ch: CHAR): BOOLEAN;
(* Return TRUE for '+' or '-' *)
BEGIN
RETURN (ch='+')OR(ch='-')
END IsSign;
(* internal state machine procedures *)
PROCEDURE RSState(inputCh: CHAR; VAR chClass: Conv.ScanClass; VAR nextState: Conv.ScanState);
BEGIN
IF Char.IsNumeric(inputCh) THEN chClass:=Conv.valid; nextState:=P
ELSE chClass:=Conv.invalid; nextState:=RS
END
END RSState;
PROCEDURE PState(inputCh: CHAR; VAR chClass: Conv.ScanClass; VAR nextState: Conv.ScanState);
BEGIN
IF Char.IsNumeric(inputCh) THEN chClass:=Conv.valid; nextState:=P
ELSIF inputCh="." THEN chClass:=Conv.valid; nextState:=F
ELSIF IsExponent(inputCh) THEN chClass:=Conv.valid; nextState:=E
ELSE chClass:=Conv.terminator; nextState:=NIL
END
END PState;
PROCEDURE FState(inputCh: CHAR; VAR chClass: Conv.ScanClass; VAR nextState: Conv.ScanState);
BEGIN
IF Char.IsNumeric(inputCh) THEN chClass:=Conv.valid; nextState:=F
ELSIF IsExponent(inputCh) THEN chClass:=Conv.valid; nextState:=E
ELSE chClass:=Conv.terminator; nextState:=NIL
END
END FState;
PROCEDURE EState(inputCh: CHAR; VAR chClass: Conv.ScanClass; VAR nextState: Conv.ScanState);
BEGIN
IF IsSign(inputCh) THEN chClass:=Conv.valid; nextState:=SE
ELSIF Char.IsNumeric(inputCh) THEN chClass:=Conv.valid; nextState:=WE
ELSE chClass:=Conv.invalid; nextState:=E
END
END EState;
PROCEDURE SEState(inputCh: CHAR; VAR chClass: Conv.ScanClass; VAR nextState: Conv.ScanState);
BEGIN
IF Char.IsNumeric(inputCh) THEN chClass:=Conv.valid; nextState:=WE
ELSE chClass:=Conv.invalid; nextState:=SE
END
END SEState;
PROCEDURE WEState(inputCh: CHAR; VAR chClass: Conv.ScanClass; VAR nextState: Conv.ScanState);
BEGIN
IF Char.IsNumeric(inputCh) THEN chClass:=Conv.valid; nextState:=WE
ELSE chClass:=Conv.terminator; nextState:=NIL
END
END WEState;
PROCEDURE Real (VAR x: LongInt; exp, digits: LONGINT; VAR outOfRange: BOOLEAN): LONGREAL;
CONST BR=LInt.B+ZERO; InvLOGB=0.221461873; (* real version *)
VAR cnt, len, scale, Bscale, start, bexp, max: LONGINT; r: LONGREAL;
BEGIN
(* scale by the exponent *)
scale:=exp+digits;
IF scale>=ABS(digits) THEN
Bscale:=0; LInt.TenPower(x, SHORT(scale))
ELSE
Bscale:=ENTIER(-scale*InvLOGB)+6;
LInt.BPower(x, SHORT(Bscale)); (* x*B^Bscale *)
LInt.TenPower(x, SHORT(scale)); (* x*B^BScale*10^scale *)
END;
(* prescale to left-justify the number *)
start:=LInt.MinDigit(x); bexp:=0; (* find starting digit *)
IF (start=LEN(x)-1)&(x[start]=0) THEN (* exit here for zero *)
outOfRange:=FALSE; RETURN ZERO
END;
WHILE x[start]<LInt.B DIV 2 DO
LInt.MultDigit(x, 2, 0); INC(bexp) (* normalize *)
END;
(* convert to a LONGREAL *)
r:=ZERO; len:=LEN(x)-1; max:=start+3;
IF max>len THEN max:=len END;
FOR cnt:=start TO max DO r:=r*BR+x[cnt] END;
(* post scaling *)
INC(bexp, (Bscale-len+max)*15);
(* quick check for overflow *)
max:=Low.exponent(r)-SHORT(bexp);
IF (max>Low.expoMax) OR (max<Low.expoMin) THEN
outOfRange:=TRUE;
RETURN ZERO
ELSE
outOfRange:=FALSE;
RETURN Low.scale(r, -SHORT(bexp))
END
END Real;
PROCEDURE ScanReal*(inputCh: CHAR; VAR chClass: Conv.ScanClass; VAR nextState: Conv.ScanState);
(*
Represents the start state of a finite state scanner for real numbers - assigns
class of inputCh to chClass and a procedure representing the next state to
nextState.
The call of ScanReal(inputCh,chClass,nextState) shall assign values to
`chClass' and `nextState' depending upon the value of `inputCh' as
shown in the following table.
Procedure inputCh chClass nextState (a procedure
with behaviour of)
--------- --------- -------- ---------
ScanReal space padding ScanReal
sign valid RSState
decimal digit valid PState
other invalid ScanReal
RSState decimal digit valid PState
other invalid RSState
PState decimal digit valid PState
"." valid FState
"E", "D" valid EState
other terminator --
FState decimal digit valid FState
"E", "D" valid EState
other terminator --
EState sign valid SEState
decimal digit valid WEState
other invalid EState
SEState decimal digit valid WEState
other invalid SEState
WEState decimal digit valid WEState
other terminator --
For examples of how to use ScanReal, refer to FormatReal and
ValueReal below.
*)
BEGIN
IF Char.IsWhiteSpace(inputCh) THEN chClass:=Conv.padding; nextState:=SR
ELSIF IsSign(inputCh) THEN chClass:=Conv.valid; nextState:=RS
ELSIF Char.IsNumeric(inputCh) THEN chClass:=Conv.valid; nextState:=P
ELSE chClass:=Conv.invalid; nextState:=SR
END
END ScanReal;
PROCEDURE FormatReal*(str: ARRAY OF CHAR): ConvResults;
(* Returns the format of the string value for conversion to LONGREAL. *)
VAR
ch: CHAR;
rn: LONGREAL;
len, index, digit, nexp, exp: INTEGER;
state: Conv.ScanState;
inExp, posExp, decExp, outOfRange: BOOLEAN;
prev, class: Conv.ScanClass;
int: LongInt;
BEGIN
state:=SR; rn:=0.0; exp:=0; nexp:= 0;
class:=Conv.padding; prev:=class;
inExp:=FALSE; posExp:=TRUE; decExp:=FALSE;
(*FOR len:=0 TO SHORT(LEN(int))-1 DO int[len]:=0 END;*)
FOR len:=0 TO (LEN(int))-1 DO int[len]:=0 END; (* I don't understand why to SHORT it. LEN(int) will return 170 (defined in LongInts as ARRAY 170 OF INTEGER) both with voc and oo2c the same way; -- noch *)
len:=Str.Length(str); index:=0;
LOOP
IF index=len THEN EXIT END;
ch:=str[index];
state.p(ch, class, state);
CASE class OF
| Conv.padding: (* nothing to do *)
| Conv.valid:
IF inExp THEN
IF IsSign(ch) THEN posExp:=ch="+"
ELSE (* must be digits *)
digit:=ORD(ch)-ORD("0");
IF posExp THEN exp:=exp*10+digit
ELSE exp:=exp*10-digit
END
END
ELSIF IsExponent(ch) THEN inExp:=TRUE
ELSIF ch="." THEN decExp:=TRUE
ELSE (* must be a digit *)
LInt.MultDigit(int, 10, ORD(ch)-ORD("0"));
IF decExp THEN DEC(nexp) END;
END
| Conv.invalid, Conv.terminator: EXIT
END;
prev:=class; INC(index)
END;
IF class IN {Conv.invalid, Conv.terminator} THEN
RETURN strWrongFormat
ELSIF prev=Conv.padding THEN
RETURN strEmpty
ELSE
rn:=Real(int, exp, nexp, outOfRange);
IF outOfRange THEN RETURN strOutOfRange
ELSE RETURN strAllRight
END
END
END FormatReal;
PROCEDURE ValueReal*(str: ARRAY OF CHAR): LONGREAL;
VAR
ch: CHAR;
rn: LONGREAL;
len, index, digit, nexp, exp: INTEGER;
state: Conv.ScanState;
inExp, positive, posExp, decExp, outOfRange: BOOLEAN;
prev, class: Conv.ScanClass;
int: LongInt;
BEGIN
state:=SR; rn:=0.0; exp:=0; nexp:= 0;
class:=Conv.padding; prev:=class;
positive:=TRUE; inExp:=FALSE; posExp:=TRUE; decExp:=FALSE;
(*FOR len:=0 TO SHORT(LEN(int))-1 DO int[len]:=0 END;*)
FOR len:=0 TO (LEN(int))-1 DO int[len]:=0 END; (* I don't understand why to SHORT it; -- noch *)
len:=Str.Length(str); index:=0;
LOOP
IF index=len THEN EXIT END;
ch:=str[index];
state.p(ch, class, state);
CASE class OF
| Conv.padding: (* nothing to do *)
| Conv.valid:
IF inExp THEN
IF IsSign(ch) THEN posExp:=ch="+"
ELSE (* must be digits *)
digit:=ORD(ch)-ORD("0");
IF posExp THEN exp:=exp*10+digit
ELSE exp:=exp*10-digit
END
END
ELSIF IsExponent(ch) THEN inExp:=TRUE
ELSIF IsSign(ch) THEN positive:=ch="+"
ELSIF ch="." THEN decExp:=TRUE
ELSE (* must be a digit *)
LInt.MultDigit(int, 10, ORD(ch)-ORD("0"));
IF decExp THEN DEC(nexp) END;
END
| Conv.invalid, Conv.terminator: EXIT
END;
prev:=class; INC(index)
END;
IF class IN {Conv.invalid, Conv.terminator} THEN
RETURN ZERO
ELSIF prev=Conv.padding THEN
RETURN ZERO
ELSE
rn:=Real(int, exp, nexp, outOfRange);
IF outOfRange THEN RETURN Low.large END
END;
IF ~positive THEN rn:=-rn END;
RETURN rn
END ValueReal;
PROCEDURE LengthFloatReal*(real: LONGREAL; sigFigs: INTEGER): INTEGER;
(*
Returns the number of characters in the floating-point string
representation of real with sigFigs significant figures.
This value corresponds to the capacity of an array `str' which
is of the minimum capacity needed to avoid truncation of the
result in the call LongStr.RealToFloat(real,sigFigs,str).
*)
VAR
len, exp: INTEGER;
BEGIN
IF Low.IsNaN(real) THEN RETURN 3
ELSIF Low.IsInfinity(real) THEN
IF real<ZERO THEN RETURN 9 ELSE RETURN 8 END
END;
IF sigFigs=0 THEN sigFigs:=SigFigs END; len:=sigFigs; (* default digits -- if none given *)
IF real<ZERO THEN INC(len); real:=-real END; (* account for the sign *)
exp:=Low.exponent10(real);
IF sigFigs>1 THEN INC(len) END; (* account for the decimal point *)
IF exp>10 THEN INC(len, 4) (* account for the exponent *)
ELSIF exp#0 THEN INC(len, 3)
END;
RETURN len
END LengthFloatReal;
PROCEDURE LengthEngReal*(real: LONGREAL; sigFigs: INTEGER): INTEGER;
(*
Returns the number of characters in the floating-point engineering
string representation of real with sigFigs significant figures.
This value corresponds to the capacity of an array `str' which is
of the minimum capacity needed to avoid truncation of the result in
the call LongStr.RealToEng(real,sigFigs,str).
*)
VAR
len, exp, off: INTEGER;
BEGIN
IF Low.IsNaN(real) THEN RETURN 3
ELSIF Low.IsInfinity(real) THEN
IF real<ZERO THEN RETURN 9 ELSE RETURN 8 END
END;
IF sigFigs=0 THEN sigFigs:=SigFigs END; len:=sigFigs; (* default digits -- if none given *)
IF real<ZERO THEN INC(len); real:=-real END; (* account for the sign *)
exp:=Low.exponent10(real); off:=exp MOD 3; (* account for the exponent *)
IF exp-off>10 THEN INC(len, 4)
ELSIF exp-off#0 THEN INC(len, 3)
END;
IF sigFigs>off+1 THEN INC(len) END; (* account for the decimal point *)
IF off+1-sigFigs>0 THEN INC(len, off+1-sigFigs) END; (* account for extra padding digits *)
RETURN len
END LengthEngReal;
PROCEDURE LengthFixedReal*(real: LONGREAL; place: INTEGER): INTEGER;
(* Returns the number of characters in the fixed-point string
representation of real rounded to the given place relative
to the decimal point.
This value corresponds to the capacity of an array `str' which
is of the minimum capacity needed to avoid truncation of the
result in the call LongStr.RealToFixed(real,sigFigs,str).
*)
VAR
len, exp: INTEGER; addDecPt: BOOLEAN;
BEGIN
IF Low.IsNaN(real) THEN RETURN 3
ELSIF Low.IsInfinity(real) THEN
IF real<ZERO THEN RETURN 9 ELSE RETURN 8 END
END;
exp:=Low.exponent10(real); addDecPt:=place>=0;
IF place<0 THEN INC(place, 2) ELSE INC(place) END;
IF exp<0 THEN (* account for digits *)
IF place<=0 THEN len:=1 ELSE len:=place END
ELSE len:=exp+place;
IF 1-place>0 THEN INC(len, 1-place) END
END;
IF real<ZERO THEN INC(len) END; (* account for the sign *)
IF addDecPt THEN INC(len) END; (* account for decimal point *)
RETURN len
END LengthFixedReal;
PROCEDURE IsRConvException*(): BOOLEAN;
(* Returns TRUE if the current coroutine is in the exceptional execution state because
of the raising of the RealConv exception; otherwise returns FALSE.
*)
BEGIN
RETURN FALSE
END IsRConvException;
PROCEDURE Test;
VAR res: INTEGER; f: LONGREAL;
BEGIN
f:=ValueReal("-1.8770465240919248E+246");
f:=ValueReal("5.1059259362558051E-111");
f:=ValueReal("2.4312432637500083E-88");
res:=LengthFixedReal(100, 0);
res:=LengthEngReal(100, 0);
res:=LengthFloatReal(100, 0);
res:=LengthFixedReal(-100.123, 0);
res:=LengthEngReal(-100.123, 0);
res:=LengthFloatReal(-100.123, 0);
res:=LengthFixedReal(-1.0D20, 0);
res:=LengthEngReal(-1.0D20, 0);
res:=LengthFloatReal(-1.0D20, 0)
END Test;
BEGIN
NEW(RS); NEW(P); NEW(F); NEW(E); NEW(SE); NEW(WE); NEW(SR);
RS.p:=RSState; P.p:=PState; F.p:=FState; E.p:=EState;
SE.p:=SEState; WE.p:=WEState; SR.p:=ScanReal;
IF DEBUG THEN Test END
END oocLRealConv.

561
src/lib/ooc/oocLRealMath.Mod Normal file
View file

@ -0,0 +1,561 @@
MODULE oocLRealMath;
(*
LRealMath - Target independent mathematical functions for LONGREAL
(IEEE double-precision) numbers.
Numerical approximations are taken from "Software Manual for the
Elementary Functions" by Cody & Waite and "Computer Approximations"
by Hart et al.
Copyright (C) 1996-1998 Michael Griebling
This module is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.
This module is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*)
IMPORT l := oocLowLReal, m := oocRealMath;
CONST
pi* = 3.1415926535897932384626433832795028841972D0;
exp1* = 2.7182818284590452353602874713526624977572D0;
ZERO=0.0D0; ONE=1.0D0; HALF=0.5D0; TWO=2.0D0; (* local constants *)
(* internally-used constants *)
huge=l.large; (* largest number this package accepts *)
miny=l.small; (* smallest number this package accepts *)
sqrtHalf=0.70710678118654752440D0;
Limit=1.0536712D-8; (* 2**(-MantBits/2) *)
eps=5.5511151D-17; (* 2**(-MantBits-1) *)
piInv=0.31830988618379067154D0; (* 1/pi *)
piByTwo=1.57079632679489661923D0;
lnv=0.6931610107421875D0; (* should be exact *)
vbytwo=0.13830277879601902638D-4; (* used in sinh/cosh *)
ln2Inv=1.44269504088896340735992468100189213D0;
(* error/exception codes *)
NoError*=m.NoError; IllegalRoot*=m.IllegalRoot; IllegalLog*=m.IllegalLog; Overflow*=m.Overflow;
IllegalPower*=m.IllegalPower; IllegalLogBase*=m.IllegalLogBase; IllegalTrig*=m.IllegalTrig;
IllegalInvTrig*=m.IllegalInvTrig; HypInvTrigClipped*=m.HypInvTrigClipped;
IllegalHypInvTrig*=m.IllegalHypInvTrig; LossOfAccuracy*=m.LossOfAccuracy;
VAR
a1: ARRAY 18 OF LONGREAL; (* lookup table for power function *)
a2: ARRAY 9 OF LONGREAL; (* lookup table for power function *)
em: LONGREAL; (* largest number such that 1+epsilon > 1.0 *)
LnInfinity: LONGREAL; (* natural log of infinity *)
LnSmall: LONGREAL; (* natural log of very small number *)
SqrtInfinity: LONGREAL; (* square root of infinity *)
TanhMax: LONGREAL; (* maximum Tanh value *)
t: LONGREAL; (* internal variables *)
(* internally used support routines *)
PROCEDURE SinCos (x, y, sign: LONGREAL): LONGREAL;
CONST
ymax=210828714; (* ENTIER(pi*2**(MantBits/2)) *)
c1=3.1416015625D0;
c2=-8.908910206761537356617D-6;
r1=-0.16666666666666665052D+0;
r2= 0.83333333333331650314D-2;
r3=-0.19841269841201840457D-3;
r4= 0.27557319210152756119D-5;
r5=-0.25052106798274584544D-7;
r6= 0.16058936490371589114D-9;
r7=-0.76429178068910467734D-12;
r8= 0.27204790957888846175D-14;
VAR
n: LONGINT; xn, f, x1, g: LONGREAL;
BEGIN
IF y>=ymax THEN l.ErrorHandler(LossOfAccuracy); RETURN ZERO END;
(* determine the reduced number *)
n:=ENTIER(y*piInv+HALF); xn:=n;
IF ODD(n) THEN sign:=-sign END;
x:=ABS(x);
IF x#y THEN xn:=xn-HALF END;
(* fractional part of reduced number *)
x1:=ENTIER(x);
f:=((x1-xn*c1)+(x-x1))-xn*c2;
(* Pre: |f| <= pi/2 *)
IF ABS(f)<Limit THEN RETURN sign*f END;
(* evaluate polynomial approximation of sin *)
g:=f*f; g:=(((((((r8*g+r7)*g+r6)*g+r5)*g+r4)*g+r3)*g+r2)*g+r1)*g;
g:=f+f*g; (* don't use less accurate f(1+g) *)
RETURN sign*g
END SinCos;
PROCEDURE div (x, y : LONGINT) : LONGINT;
(* corrected MOD function *)
BEGIN
IF x < 0 THEN RETURN -ABS(x) DIV y ELSE RETURN x DIV y END
END div;
(* forward declarations *)
PROCEDURE^ arctan2* (xn, xd: LONGREAL): LONGREAL;
PROCEDURE^ sincos* (x: LONGREAL; VAR Sin, Cos: LONGREAL);
PROCEDURE sqrt*(x: LONGREAL): LONGREAL;
(* Returns the positive square root of x where x >= 0 *)
CONST
P0=0.41731; P1=0.59016;
VAR
xMant, yEst, z: LONGREAL; xExp: INTEGER;
BEGIN
(* optimize zeros and check for illegal negative roots *)
IF x=ZERO THEN RETURN ZERO END;
IF x<ZERO THEN l.ErrorHandler(IllegalRoot); x:=-x END;
(* reduce the input number to the range 0.5 <= x <= 1.0 *)
xMant:=l.fraction(x)*HALF; xExp:=l.exponent(x)+1;
(* initial estimate of the square root *)
yEst:=P0+P1*xMant;
(* perform three newtonian iterations *)
z:=(yEst+xMant/yEst); yEst:=0.25*z+xMant/z;
yEst:=HALF*(yEst+xMant/yEst);
(* adjust for odd exponents *)
IF ODD(xExp) THEN yEst:=yEst*sqrtHalf; INC(xExp) END;
(* single Newtonian iteration to produce real number accuracy *)
RETURN l.scale(yEst, xExp DIV 2)
END sqrt;
PROCEDURE exp*(x: LONGREAL): LONGREAL;
(* Returns the exponential of x for x < Ln(MAX(REAL) *)
CONST
c1=0.693359375D0; c2=-2.1219444005469058277D-4;
P0=0.249999999999999993D+0; P1=0.694360001511792852D-2; P2=0.165203300268279130D-4;
Q1=0.555538666969001188D-1; Q2=0.495862884905441294D-3;
VAR xn, g, p, q, z: LONGREAL; n: INTEGER;
BEGIN
(* Ensure we detect overflows and return 0 for underflows *)
IF x>LnInfinity THEN l.ErrorHandler(Overflow); RETURN huge
ELSIF x<LnSmall THEN RETURN ZERO
ELSIF ABS(x)<eps THEN RETURN ONE
END;
(* Decompose and scale the number *)
IF x>=ZERO THEN n:=SHORT(ENTIER(ln2Inv*x+HALF))
ELSE n:=SHORT(ENTIER(ln2Inv*x-HALF))
END;
xn:=n; g:=(x-xn*c1)-xn*c2;
(* Calculate exp(g)/2 from "Software Manual for the Elementary Functions" *)
z:=g*g; p:=((P2*z+P1)*z+P0)*g; q:=(Q2*z+Q1)*z+HALF;
RETURN l.scale(HALF+p/(q-p), n+1)
END exp;
PROCEDURE ln*(x: LONGREAL): LONGREAL;
(* Returns the natural logarithm of x for x > 0 *)
CONST
c1=355.0D0/512.0D0; c2=-2.121944400546905827679D-4;
P0=-0.64124943423745581147D+2; P1=0.16383943563021534222D+2; P2=-0.78956112887491257267D+0;
Q0=-0.76949932108494879777D+3; Q1=0.31203222091924532844D+3; Q2=-0.35667977739034646171D+2;
VAR f, zn, zd, r, z, w, p, q, xn: LONGREAL; n: INTEGER;
BEGIN
(* ensure illegal inputs are trapped and handled *)
IF x<=ZERO THEN l.ErrorHandler(IllegalLog); RETURN -huge END;
(* reduce the range of the input *)
f:=l.fraction(x)*HALF; n:=l.exponent(x)+1;
IF f>sqrtHalf THEN zn:=(f-HALF)-HALF; zd:=f*HALF+HALF
ELSE zn:=f-HALF; zd:=zn*HALF+HALF; DEC(n)
END;
(* evaluate rational approximation from "Software Manual for the Elementary Functions" *)
z:=zn/zd; w:=z*z; q:=((w+Q2)*w+Q1)*w+Q0; p:=w*((P2*w+P1)*w+P0); r:=z+z*(p/q);
(* scale the output *)
xn:=n;
RETURN (xn*c2+r)+xn*c1
END ln;
(* The angle in all trigonometric functions is measured in radians *)
PROCEDURE sin* (x: LONGREAL): LONGREAL;
BEGIN
IF x<ZERO THEN RETURN SinCos(x, -x, -ONE)
ELSE RETURN SinCos(x, x, ONE)
END
END sin;
PROCEDURE cos* (x: LONGREAL): LONGREAL;
BEGIN
RETURN SinCos(x, ABS(x)+piByTwo, ONE)
END cos;
PROCEDURE tan*(x: LONGREAL): LONGREAL;
(* Returns the tangent of x where x cannot be an odd multiple of pi/2 *)
VAR Sin, Cos: LONGREAL;
BEGIN
sincos(x, Sin, Cos);
IF ABS(Cos)<miny THEN l.ErrorHandler(IllegalTrig); RETURN huge
ELSE RETURN Sin/Cos
END
END tan;
PROCEDURE arcsin*(x: LONGREAL): LONGREAL;
(* Returns the arcsine of x, in the range [-pi/2, pi/2] where -1 <= x <= 1 *)
BEGIN
IF ABS(x)>ONE THEN l.ErrorHandler(IllegalInvTrig); RETURN huge
ELSE RETURN arctan2(x, sqrt(ONE-x*x))
END
END arcsin;
PROCEDURE arccos*(x: LONGREAL): LONGREAL;
(* Returns the arccosine of x, in the range [0, pi] where -1 <= x <= 1 *)
BEGIN
IF ABS(x)>ONE THEN l.ErrorHandler(IllegalInvTrig); RETURN huge
ELSE RETURN arctan2(sqrt(ONE-x*x), x)
END
END arccos;
PROCEDURE arctan*(x: LONGREAL): LONGREAL;
(* Returns the arctangent of x, in the range [-pi/2, pi/2] for all x *)
BEGIN
RETURN arctan2(x, ONE)
END arctan;
PROCEDURE power*(base, exponent: LONGREAL): LONGREAL;
(* Returns the value of the number base raised to the power exponent
for base > 0 *)
CONST
P1=0.83333333333333211405D-1; P2=0.12500000000503799174D-1;
P3=0.22321421285924258967D-2; P4=0.43445775672163119635D-3;
K=0.44269504088896340736D0;
Q1=0.69314718055994529629D+0; Q2=0.24022650695909537056D+0;
Q3=0.55504108664085595326D-1; Q4=0.96181290595172416964D-2;
Q5=0.13333541313585784703D-2; Q6=0.15400290440989764601D-3;
Q7=0.14928852680595608186D-4;
OneOver16=0.0625D0; XMAX=16*l.expoMax-1; (*XMIN=16*l.expoMin+1;*) XMIN=-16351; (* noch *)
VAR z, g, R, v, u2, u1, w1, w2, y1, y2, w: LONGREAL; m, p, i: INTEGER; mp, pp, iw1: LONGINT;
BEGIN
(* handle all possible error conditions *)
IF ABS(exponent)<miny THEN RETURN ONE (* base**0 = 1 *)
ELSIF base<ZERO THEN l.ErrorHandler(IllegalPower); RETURN -huge
ELSIF ABS(base)<miny THEN
IF exponent>ZERO THEN RETURN ZERO ELSE l.ErrorHandler(Overflow); RETURN -huge END
END;
(* extract the exponent of base to m and clear exponent of base in g *)
g:=l.fraction(base)*HALF; m:=l.exponent(base)+1;
(* determine p table offset with an unrolled binary search *)
p:=1;
IF g<=a1[9] THEN p:=9 END;
IF g<=a1[p+4] THEN INC(p, 4) END;
IF g<=a1[p+2] THEN INC(p, 2) END;
(* compute scaled z so that |z| <= 0.044 *)
z:=((g-a1[p+1])-a2[(p+1) DIV 2])/(g+a1[p+1]); z:=z+z;
(* approximation for log2(z) from "Software Manual for the Elementary Functions" *)
v:=z*z; R:=(((P4*v+P3)*v+P2)*v+P1)*v*z; R:=R+K*R; u2:=(R+z*K)+z; u1:=(m*16-p)*OneOver16;
(* generate w with extra precision calculations *)
y1:=ENTIER(16*exponent)*OneOver16; y2:=exponent-y1; w:=u2*exponent+u1*y2;
w1:=ENTIER(16*w)*OneOver16; w2:=w-w1; w:=w1+u1*y1;
w1:=ENTIER(16*w)*OneOver16; w2:=w2+(w-w1); w:=ENTIER(16*w2)*OneOver16;
iw1:=ENTIER(16*(w+w1)); w2:=w2-w;
(* check for overflow/underflow *)
IF iw1>XMAX THEN l.ErrorHandler(Overflow); RETURN huge
ELSIF iw1<XMIN THEN RETURN ZERO (* underflow *)
END;
(* final approximation 2**w2-1 where -0.0625 <= w2 <= 0 *)
IF w2>ZERO THEN INC(iw1); w2:=w2-OneOver16 END; IF iw1<0 THEN i:=0 ELSE i:=1 END;
mp:=div(iw1, 16)+i; pp:=16*mp-iw1;
z:=((((((Q7*w2+Q6)*w2+Q5)*w2+Q4)*w2+Q3)*w2+Q2)*w2+Q1)*w2; z:=a1[pp+1]+a1[pp+1]*z;
RETURN l.scale(z, SHORT(mp))
END power;
PROCEDURE round*(x: LONGREAL): LONGINT;
(* Returns the value of x rounded to the nearest integer *)
BEGIN
IF x<ZERO THEN RETURN -ENTIER(HALF-x)
ELSE RETURN ENTIER(x+HALF)
END
END round;
PROCEDURE IsRMathException*(): BOOLEAN;
(* Returns TRUE if the current coroutine is in the exceptional execution state
because of the raising of the RealMath exception; otherwise returns FALSE.
*)
BEGIN
RETURN FALSE
END IsRMathException;
(*
Following routines are provided as extensions to the ISO standard.
They are either used as the basis of other functions or provide
useful functions which are not part of the ISO standard.
*)
PROCEDURE log* (x, base: LONGREAL): LONGREAL;
(* log(x,base) is the logarithm of x base b. All positive arguments are
allowed but base > 0 and base # 1. *)
BEGIN
(* log(x, base) = log2(x) / log2(base) *)
IF base<=ZERO THEN l.ErrorHandler(IllegalLogBase); RETURN -huge
ELSE RETURN ln(x)/ln(base)
END
END log;
PROCEDURE ipower* (x: LONGREAL; base: INTEGER): LONGREAL;
(* ipower(x, base) returns the x to the integer power base where base*Log2(x) < Log2(Max) *)
VAR y: LONGREAL; neg: BOOLEAN; Exp: LONGINT;
PROCEDURE Adjust(xadj: LONGREAL): LONGREAL;
BEGIN
IF (x<ZERO)&ODD(base) THEN RETURN -xadj ELSE RETURN xadj END
END Adjust;
BEGIN
(* handle all possible error conditions *)
IF base=0 THEN RETURN ONE (* x**0 = 1 *)
ELSIF ABS(x)<miny THEN
IF base>0 THEN RETURN ZERO ELSE l.ErrorHandler(Overflow); RETURN Adjust(huge) END
END;
(* trap potential overflows and underflows *)
Exp:=(l.exponent(x)+1)*base; y:=LnInfinity*ln2Inv;
IF Exp>y THEN l.ErrorHandler(Overflow); RETURN Adjust(huge)
ELSIF Exp<-y THEN RETURN ZERO
END;
(* compute x**base using an optimised algorithm from Knuth, slightly
altered : p442, The Art Of Computer Programming, Vol 2 *)
y:=ONE; IF base<0 THEN neg:=TRUE; base := -base ELSE neg:= FALSE END;
LOOP
IF ODD(base) THEN y:=y*x END;
base:=base DIV 2; IF base=0 THEN EXIT END;
x:=x*x;
END;
IF neg THEN RETURN ONE/y ELSE RETURN y END
END ipower;
PROCEDURE sincos* (x: LONGREAL; VAR Sin, Cos: LONGREAL);
(* More efficient sin/cos implementation if both values are needed. *)
BEGIN
Sin:=sin(x); Cos:=sqrt(ONE-Sin*Sin)
END sincos;
PROCEDURE arctan2* (xn, xd: LONGREAL): LONGREAL;
(* arctan2(xn,xd) is the quadrant-correct arc tangent atan(xn/xd). If the
denominator xd is zero, then the numerator xn must not be zero. All
arguments are legal except xn = xd = 0. *)
CONST
P0=0.216062307897242551884D+3; P1=0.3226620700132512059245D+3;
P2=0.13270239816397674701D+3; P3=0.1288838303415727934D+2;
Q0=0.2160623078972426128957D+3; Q1=0.3946828393122829592162D+3;
Q2=0.221050883028417680623D+3; Q3=0.3850148650835119501D+2;
PiOver2=pi/2; Sqrt3=1.7320508075688772935D0;
VAR atan, z, z2, p, q: LONGREAL; xnExp, xdExp: INTEGER; Quadrant: SHORTINT;
BEGIN
IF ABS(xd)<miny THEN
IF ABS(xn)<miny THEN l.ErrorHandler(IllegalInvTrig); atan:=ZERO
ELSE l.ErrorHandler(Overflow); atan:=PiOver2
END
ELSE xnExp:=l.exponent(xn); xdExp:=l.exponent(xd);
IF xnExp-xdExp>=l.expoMax-3 THEN l.ErrorHandler(Overflow); atan:=PiOver2
ELSIF xnExp-xdExp<l.expoMin+3 THEN atan:=ZERO
ELSE
(* ensure division of xn/xd always produces a number < 1 & resolve quadrant *)
IF ABS(xn)>ABS(xd) THEN z:=ABS(xd/xn); Quadrant:=2
ELSE z:=ABS(xn/xd); Quadrant:=0
END;
(* further reduce range to within 0 to 2-sqrt(3) *)
IF z>TWO-Sqrt3 THEN z:=(z*Sqrt3-ONE)/(Sqrt3+z); INC(Quadrant) END;
(* approximation from "Computer Approximations" table ARCTN 5075 *)
IF ABS(z)<Limit THEN atan:=z (* for small values of z2, return this value *)
ELSE z2:=z*z; p:=(((P3*z2+P2)*z2+P1)*z2+P0)*z; q:=(((z2+Q3)*z2+Q2)*z2+Q1)*z2+Q0; atan:=p/q;
END;
(* adjust for z's quadrant *)
IF Quadrant>1 THEN atan:=-atan END;
CASE Quadrant OF
1: atan:=atan+pi/6
| 2: atan:=atan+PiOver2
| 3: atan:=atan+pi/3
| ELSE (* angle is correct *)
END
END;
(* map negative xds into the correct quadrant *)
IF xd<ZERO THEN atan:=pi-atan END
END;
(* map negative xns into the correct quadrant *)
IF xn<ZERO THEN atan:=-atan END;
RETURN atan
END arctan2;
PROCEDURE sinh* (x: LONGREAL): LONGREAL;
(* sinh(x) is the hyperbolic sine of x. The argument x must not be so large
that exp(|x|) overflows. *)
CONST
P0=-0.35181283430177117881D+6; P1=-0.11563521196851768270D+5;
P2=-0.16375798202630751372D+3; P3=-0.78966127417357099479D+0;
Q0=-0.21108770058106271242D+7; Q1= 0.36162723109421836460D+5;
Q2=-0.27773523119650701667D+3;
VAR y, f, p, q: LONGREAL;
BEGIN y:=ABS(x);
IF y<=ONE THEN (* handle small arguments *)
IF y<Limit THEN RETURN x END;
(* use approximation from "Software Manual for the Elementary Functions" *)
f:=y*y; p:=((P3*f+P2)*f+P1)*f+P0; q:=((f+Q2)*f+Q1)*f+Q0; y:=f*(p/q); RETURN x+x*y
ELSIF y>LnInfinity THEN (* handle exp overflows *)
y:=y-lnv;
IF y>LnInfinity-lnv+0.69 THEN l.ErrorHandler(Overflow);
IF x>ZERO THEN RETURN huge ELSE RETURN -huge END
ELSE f:=exp(y); f:=f+f*vbytwo (* don't change to f(1+vbytwo) *)
END
ELSE f:=exp(y); f:=(f-ONE/f)*HALF
END;
(* reach here when 1 < ABS(x) < LnInfinity-lnv+0.69 *)
IF x>ZERO THEN RETURN f ELSE RETURN -f END
END sinh;
PROCEDURE cosh* (x: LONGREAL): LONGREAL;
(* cosh(x) is the hyperbolic cosine of x. The argument x must not be so large
that exp(|x|) overflows. *)
VAR y, f: LONGREAL;
BEGIN y:=ABS(x);
IF y>LnInfinity THEN (* handle exp overflows *)
y:=y-lnv;
IF y>LnInfinity-lnv+0.69 THEN l.ErrorHandler(Overflow);
IF x>ZERO THEN RETURN huge ELSE RETURN -huge END
ELSE f:=exp(y); RETURN f+f*vbytwo (* don't change to f(1+vbytwo) *)
END
ELSE f:=exp(y); RETURN (f+ONE/f)*HALF
END
END cosh;
PROCEDURE tanh* (x: LONGREAL): LONGREAL;
(* tanh(x) is the hyperbolic tangent of x. All arguments are legal. *)
CONST
P0=-0.16134119023996228053D+4; P1=-0.99225929672236083313D+2; P2=-0.96437492777225469787D+0;
Q0= 0.48402357071988688686D+4; Q1= 0.22337720718962312926D+4; Q2= 0.11274474380534949335D+3;
ln3over2=0.54930614433405484570D0;
BIG=19.06154747D0; (* (ln(2)+(t+1)*ln(B))/2 where t=mantissa bits, B=base *)
VAR f, t: LONGREAL;
BEGIN f:=ABS(x);
IF f>BIG THEN t:=ONE
ELSIF f>ln3over2 THEN t:=ONE-TWO/(exp(TWO*f)+ONE)
ELSIF f<Limit THEN t:=f
ELSE (* approximation from "Software Manual for the Elementary Functions" *)
t:=f*f; t:=t*(((P2*t+P1)*t+P0)/(((t+Q2)*t+Q1)*t+Q0)); t:=f+f*t
END;
IF x<ZERO THEN RETURN -t ELSE RETURN t END
END tanh;
PROCEDURE arcsinh* (x: LONGREAL): LONGREAL;
(* arcsinh(x) is the arc hyperbolic sine of x. All arguments are legal. *)
BEGIN
IF ABS(x)>SqrtInfinity*HALF THEN l.ErrorHandler(HypInvTrigClipped);
IF x>ZERO THEN RETURN ln(SqrtInfinity) ELSE RETURN -ln(SqrtInfinity) END;
ELSIF x<ZERO THEN RETURN -ln(-x+sqrt(x*x+ONE))
ELSE RETURN ln(x+sqrt(x*x+ONE))
END
END arcsinh;
PROCEDURE arccosh* (x: LONGREAL): LONGREAL;
(* arccosh(x) is the arc hyperbolic cosine of x. All arguments greater than
or equal to 1 are legal. *)
BEGIN
IF x<ONE THEN l.ErrorHandler(IllegalHypInvTrig); RETURN ZERO
ELSIF x>SqrtInfinity*HALF THEN l.ErrorHandler(HypInvTrigClipped); RETURN ln(SqrtInfinity)
ELSE RETURN ln(x+sqrt(x*x-ONE))
END
END arccosh;
PROCEDURE arctanh* (x: LONGREAL): LONGREAL;
(* arctanh(x) is the arc hyperbolic tangent of x. |x| < 1 - sqrt(em), where
em is machine epsilon. Note that |x| must not be so close to 1 that the
result is less accurate than half precision. *)
CONST TanhLimit=0.999984991D0; (* Tanh(5.9) *)
VAR t: LONGREAL;
BEGIN t:=ABS(x);
IF (t>=ONE) OR (t>(ONE-TWO*em)) THEN l.ErrorHandler(IllegalHypInvTrig);
IF x<ZERO THEN RETURN -TanhMax ELSE RETURN TanhMax END
ELSIF t>TanhLimit THEN l.ErrorHandler(LossOfAccuracy)
END;
RETURN arcsinh(x/sqrt(ONE-x*x))
END arctanh;
PROCEDURE ToLONGREAL (hi, lo: LONGINT): LONGREAL;
VAR ra: ARRAY 2 OF LONGINT;
BEGIN ra[0]:=hi; ra[1]:=lo;
RETURN l.Real(ra)
END ToLONGREAL;
BEGIN
(* determine some fundamental constants used by hyperbolic trig functions *)
em:=l.ulp(ONE);
LnInfinity:=ln(huge);
LnSmall:=ln(miny);
SqrtInfinity:=sqrt(huge);
t:=l.pred(ONE)/sqrt(em); TanhMax:=ln(t+sqrt(t*t+ONE));
(* initialize some tables for the power() function a1[i]=2**((1-i)/16) *)
(* disable compiler warnings about 32-bit negative integers *)
(*<* PUSH; Warnings := FALSE *>*)
a1[ 1]:=ONE;
a1[ 2]:=ToLONGREAL(3FEEA4AFH, 0A2A490DAH);
a1[ 3]:=ToLONGREAL(3FED5818H, 0DCFBA487H);
a1[ 4]:=ToLONGREAL(3FEC199BH, 0DD85529CH);
a1[ 5]:=ToLONGREAL(3FEAE89FH, 0995AD3ADH);
a1[ 6]:=ToLONGREAL(3FE9C491H, 082A3F090H);
a1[ 7]:=ToLONGREAL(3FE8ACE5H, 0422AA0DBH);
a1[ 8]:=ToLONGREAL(3FE7A114H, 073EB0186H);
a1[ 9]:=ToLONGREAL(3FE6A09EH, 0667F3BCCH);
a1[10]:=ToLONGREAL(3FE5AB07H, 0DD485429H);
a1[11]:=ToLONGREAL(3FE4BFDAH, 0D5362A27H);
a1[12]:=ToLONGREAL(3FE3DEA6H, 04C123422H);
a1[13]:=ToLONGREAL(3FE306FEH, 00A31B715H);
a1[14]:=ToLONGREAL(3FE2387AH, 06E756238H);
a1[15]:=ToLONGREAL(3FE172B8H, 03C7D517AH);
a1[16]:=ToLONGREAL(3FE0B558H, 06CF9890FH);
a1[17]:=HALF;
(* a2[i]=2**[(1-2i)/16] - a1[2i]; delta resolution *)
a2[1]:=ToLONGREAL(3C90B1EEH, 074320000H);
a2[2]:=ToLONGREAL(3C711065H, 089500000H);
a2[3]:=ToLONGREAL(3C6C7C46H, 0B0700000H);
a2[4]:=ToLONGREAL(3C9AFAA2H, 0047F0000H);
a2[5]:=ToLONGREAL(3C86324CH, 005460000H);
a2[6]:=ToLONGREAL(3C7ADA09H, 011F00000H);
a2[7]:=ToLONGREAL(3C89B07EH, 0B6C80000H);
a2[8]:=ToLONGREAL(3C88A62EH, 04ADC0000H);
(* reenable compiler warnings *)
(*<* POP *>*)
END oocLRealMath.

451
src/lib/ooc/oocLRealStr.Mod Normal file
View file

@ -0,0 +1,451 @@
(* $Id: LRealStr.Mod,v 1.8 2001/07/15 14:59:29 ooc-devel Exp $ *)
MODULE oocLRealStr;
(*
LRealStr - LONGREAL/string conversions.
Copyright (C) 1996, 2001 Michael Griebling
This module is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.
This module is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*)
IMPORT
Low := oocLowLReal, Conv := oocConvTypes, RC := oocLRealConv, Str := oocStrings,
LInt := oocLongInts;
CONST
ZERO=0.0D0; B=8000H;
TYPE
ConvResults*= Conv.ConvResults; (* strAllRight, strOutOfRange, strWrongFormat, strEmpty *)
CONST
strAllRight*=Conv.strAllRight; (* the string format is correct for the corresponding conversion *)
strOutOfRange*=Conv.strOutOfRange; (* the string is well-formed but the value cannot be represented *)
strWrongFormat*=Conv.strWrongFormat; (* the string is in the wrong format for the conversion *)
strEmpty*=Conv.strEmpty; (* the given string is empty *)
(* the string form of a signed fixed-point real number is
["+" | "-"], decimal digit, {decimal digit}, [".", {decimal digit}]
*)
(* the string form of a signed floating-point real number is
signed fixed-point real number, "E"|"e", ["+" | "-"], decimal digit, {decimal digit}
*)
PROCEDURE StrToReal*(str: ARRAY OF CHAR; VAR real: LONGREAL; VAR res: ConvResults);
(*
Ignores any leading spaces in str. If the subsequent characters in str
are in the format of a signed real number, and shall assign values to
`res' and `real' as follows:
strAllRight
if the remainder of `str' represents a complete signed real number
in the range of the type of `real' -- the value of this number shall
be assigned to `real';
strOutOfRange
if the remainder of `str' represents a complete signed real number
but its value is out of the range of the type of `real' -- the
maximum or minimum value of the type of `real' shall be assigned to
`real' according to the sign of the number;
strWrongFormat
if there are remaining characters in `str' but these are not in the
form of a complete signed real number -- the value of `real' is not
defined;
strEmpty
if there are no remaining characters in `str' -- the value of `real'
is not defined.
*)
BEGIN
res:=RC.FormatReal(str);
IF res IN {strAllRight, strOutOfRange} THEN real:=RC.ValueReal(str) END
END StrToReal;
PROCEDURE AppendChar(ch: CHAR; VAR str: ARRAY OF CHAR);
VAR ds: ARRAY 2 OF CHAR;
BEGIN
ds[0]:=ch; ds[1]:=0X; Str.Append(ds, str)
END AppendChar;
PROCEDURE AppendDigit(dig: LONGINT; VAR str: ARRAY OF CHAR);
BEGIN
AppendChar(CHR(dig+ORD("0")), str)
END AppendDigit;
PROCEDURE AppendExponent(exp: INTEGER; VAR str: ARRAY OF CHAR);
BEGIN
Str.Append("E", str);
IF exp<0 THEN exp:=-exp; Str.Append("-", str)
ELSE Str.Append("+", str)
END;
IF exp>=100 THEN AppendDigit(exp DIV 100, str) END;
IF exp>=10 THEN AppendDigit((exp DIV 10) MOD 10, str) END;
AppendDigit(exp MOD 10, str)
END AppendExponent;
PROCEDURE AppendFraction(VAR n: LInt.LongInt; sigFigs, place: INTEGER; VAR str: ARRAY OF CHAR);
VAR digs, end: INTEGER; d: LONGINT; lstr: ARRAY 64 OF CHAR;
BEGIN
(* write significant digits *)
lstr:="";
FOR digs:=1 TO sigFigs DO
LInt.DivDigit(n, 10, d); AppendDigit(d, lstr);
END;
(* reverse the real digits and append to str *)
end:=sigFigs-1;
FOR digs:=0 TO sigFigs-1 DO
IF digs=place THEN Str.Append(".", str) END;
AppendChar(lstr[end], str); DEC(end)
END;
(* pad out digits to the decimal position *)
FOR digs:=sigFigs TO place-1 DO Str.Append("0", str) END
END AppendFraction;
PROCEDURE RemoveLeadingZeros(VAR str: ARRAY OF CHAR);
VAR len: LONGINT;
BEGIN
len:=Str.Length(str);
WHILE (len>1)&(str[0]="0")&(str[1]#".") DO Str.Delete(str, 0, 1); DEC(len) END
END RemoveLeadingZeros;
PROCEDURE MaxDigit (VAR n: LInt.LongInt) : LONGINT;
VAR
i, max : LONGINT;
BEGIN
(* return the maximum digit in the specified LongInt number *)
FOR i:=0 TO LEN(n)-1 DO
IF n[i] # 0 THEN
max := n[i];
WHILE max>=10 DO max:=max DIV 10 END;
RETURN max;
END;
END;
RETURN 0;
END MaxDigit;
PROCEDURE Scale (x: LONGREAL; VAR n: LInt.LongInt; sigFigs: INTEGER; exp: INTEGER; VAR overflow : BOOLEAN);
CONST
MaxDigits=4; LOG2B=15;
VAR
i, m, ln, d: LONGINT; e1, e2: INTEGER;
max: LONGINT;
BEGIN
(* extract fraction & exponent *)
m:=0; overflow := FALSE;
WHILE Low.exponent(x)=Low.expoMin DO (* scale up subnormal numbers *)
x:=x*2.0D0; DEC(m)
END;
m:=m+Low.exponent(x); x:=Low.fraction(x);
x:=Low.scale(x, SHORT(m MOD LOG2B)); (* scale up the number *)
m:=m DIV LOG2B; (* base B exponent *)
(* convert to an extended integer MOD B *)
ln:=LEN(n)-1;
FOR i:=ln-MaxDigits TO ln DO
n[i]:=SHORT(ENTIER(x)); (* convert/store the number *)
x:=(x-n[i])*B
END;
FOR i:=0 TO ln-MaxDigits-1 DO n[i]:=0 END; (* zero the other digits *)
(* scale to get the number of significant digits *)
e1:=SHORT(m)-MaxDigits; e2:= sigFigs-exp-1;
IF e1>=0 THEN
LInt.BPower(n, e1+1); LInt.TenPower(n, e2);
max := MaxDigit(n); (* remember the original digit so we can check for round-up *)
LInt.AddDigit(n, B DIV 2); LInt.DivDigit(n, B, d) (* round *)
ELSIF e2>0 THEN
LInt.TenPower(n, e2);
IF e1>0 THEN LInt.BPower(n, e1-1) ELSE LInt.BPower(n, e1+1) END;
max := MaxDigit(n); (* remember the original digit so we can check for round-up *)
LInt.AddDigit(n, B DIV 2); LInt.DivDigit(n, B, d) (* round *)
ELSE (* e1<=0, e2<=0 *)
LInt.TenPower(n, e2); LInt.BPower(n, e1+1);
max := MaxDigit(n); (* remember the original digit so we can check for round-up *)
LInt.AddDigit(n, B DIV 2); LInt.DivDigit(n, B, d) (* round *)
END;
(* check if the upper digit was changed by rounding up *)
IF (max = 9) & (max # MaxDigit(n)) THEN
overflow := TRUE;
END
END Scale;
PROCEDURE RealToFloat*(real: LONGREAL; sigFigs: INTEGER; VAR str: ARRAY OF CHAR);
(*
The call RealToFloat(real,sigFigs,str) shall assign to `str' the
possibly truncated string corresponding to the value of `real' in
floating-point form. A sign shall be included only for negative
values. One significant digit shall be included in the whole number
part. The signed exponent part shall be included only if the exponent
value is not 0. If the value of `sigFigs' is greater than 0, that
number of significant digits shall be included, otherwise an
implementation-defined number of significant digits shall be
included. The decimal point shall not be included if there are no
significant digits in the fractional part.
For example:
value: 3923009 39.23009 0.0003923009
sigFigs
1 4E+6 4E+1 4E-4
2 3.9E+6 3.9E+1 3.9E-4
5 3.9230E+6 3.9230E+1 3.9230E-4
*)
VAR
exp: INTEGER; in: LInt.LongInt;
lstr: ARRAY 64 OF CHAR;
overflow: BOOLEAN;
d: LONGINT;
BEGIN
(* set significant digits, extract sign & exponent *)
lstr:="";
IF sigFigs<=0 THEN sigFigs:=RC.SigFigs END;
(* check for illegal numbers *)
IF Low.IsNaN(real) THEN COPY("NaN", str); RETURN END;
IF real<ZERO THEN Str.Append("-", lstr); real:=-real END;
IF Low.IsInfinity(real) THEN Str.Append("Infinity", lstr); COPY(lstr, str); RETURN END;
exp:=Low.exponent10(real);
(* round the number and extract exponent again *)
Scale(real, in, sigFigs, exp, overflow);
IF overflow THEN
IF exp>=0 THEN INC(exp) ELSE DEC(exp) END;
LInt.DivDigit(in, 10, d)
END;
(* output number like x[.{x}][E+n[n]] *)
AppendFraction(in, sigFigs, 1, lstr);
IF exp#0 THEN AppendExponent(exp, lstr) END;
(* possibly truncate the result *)
COPY(lstr, str)
END RealToFloat;
PROCEDURE RealToEng*(real: LONGREAL; sigFigs: INTEGER; VAR str: ARRAY OF CHAR);
(*
Converts the value of real to floating-point string form, with
sigFigs significant figures, and copies the possibly truncated
result to str. The number is scaled with one to three digits in
the whole number part and with an exponent that is a multiple of
three.
For example:
value: 3923009 39.23009 0.0003923009
sigFigs
1 4E+6 40 400E-6
2 3.9E+6 39 390E-6
5 3.9230E+6 39.230 392.30E-6
*)
VAR
in: LInt.LongInt; exp, offset: INTEGER;
lstr: ARRAY 64 OF CHAR;
d: LONGINT;
overflow: BOOLEAN;
BEGIN
(* set significant digits, extract sign & exponent *)
lstr:="";
IF sigFigs<=0 THEN sigFigs:=RC.SigFigs END;
(* check for illegal numbers *)
IF Low.IsNaN(real) THEN COPY("NaN", str); RETURN END;
IF real<ZERO THEN Str.Append("-", lstr); real:=-real END;
IF Low.IsInfinity(real) THEN Str.Append("Infinity", lstr); COPY(lstr, str); RETURN END;
exp:=Low.exponent10(real);
(* round the number and extract exponent again (ie. 9.9 => 10.0) *)
Scale(real, in, sigFigs, exp, overflow);
IF overflow THEN
IF exp>=0 THEN INC(exp) ELSE DEC(exp) END;
LInt.DivDigit(in, 10, d)
END;
(* find the offset to make the exponent a multiple of three *)
offset:=exp MOD 3;
(* output number like x[x][x][.{x}][E+n[n]] *)
AppendFraction(in, sigFigs, offset+1, lstr);
exp:=exp-offset;
IF exp#0 THEN AppendExponent(exp, lstr) END;
(* possibly truncate the result *)
COPY(lstr, str)
END RealToEng;
PROCEDURE RealToFixed*(real: LONGREAL; place: INTEGER; VAR str: ARRAY OF CHAR);
(*
The call RealToFixed(real,place,str) shall assign to `str' the
possibly truncated string corresponding to the value of `real' in
fixed-point form. A sign shall be included only for negative values.
At least one digit shall be included in the whole number part. The
value shall be rounded to the given value of `place' relative to the
decimal point. The decimal point shall be suppressed if `place' is
less than 0.
For example:
value: 3923009 3.923009 0.0003923009
sigFigs
-5 3920000 0 0
-2 3923010 0 0
-1 3923009 4 0
0 3923009. 4. 0.
1 3923009.0 3.9 0.0
4 3923009.0000 3.9230 0.0004
*)
VAR
in: LInt.LongInt; exp, digs: INTEGER;
overflow, addDecPt: BOOLEAN;
lstr: ARRAY 256 OF CHAR;
BEGIN
(* set significant digits, extract sign & exponent *)
lstr:=""; addDecPt:=place=0;
(* check for illegal numbers *)
IF Low.IsNaN(real) THEN COPY("NaN", str); RETURN END;
IF real<ZERO THEN Str.Append("-", lstr); real:=-real END;
IF Low.IsInfinity(real) THEN Str.Append("Infinity", lstr); COPY(lstr, str); RETURN END;
exp:=Low.exponent10(real);
IF place<0 THEN digs:=place+exp+2 ELSE digs:=place+exp+1 END;
(* round the number and extract exponent again (ie. 9.9 => 10.0) *)
Scale(real, in, digs, exp, overflow);
IF overflow THEN
INC(digs); INC(exp);
addDecPt := place=0;
END;
(* output number like x[{x}][.{x}] *)
IF exp<0 THEN
IF place<0 THEN AppendFraction(in, 1, 1, lstr)
ELSE AppendFraction(in, place+1, 1, lstr)
END
ELSE AppendFraction(in, digs, exp+1, lstr);
RemoveLeadingZeros(lstr)
END;
(* special formatting *)
IF addDecPt THEN Str.Append(".", lstr) END;
(* possibly truncate the result *)
COPY(lstr, str)
END RealToFixed;
PROCEDURE RealToStr*(real: LONGREAL; VAR str: ARRAY OF CHAR);
(*
If the sign and magnitude of `real' can be shown within the capacity
of `str', the call RealToStr(real,str) shall behave as the call
RealToFixed(real,place,str), with a value of `place' chosen to fill
exactly the remainder of `str'. Otherwise, the call shall behave as
the call RealToFloat(real,sigFigs,str), with a value of `sigFigs' of
at least one, but otherwise limited to the number of significant
digits that can be included together with the sign and exponent part
in `str'.
*)
VAR
cap, exp, fp, len, pos: INTEGER;
found: BOOLEAN;
BEGIN
cap:=SHORT(LEN(str))-1; (* determine the capacity of the string with space for trailing 0X *)
(* check for illegal numbers *)
IF Low.IsNaN(real) THEN COPY("NaN", str); RETURN END;
IF real<ZERO THEN COPY("-", str); fp:=-1 ELSE COPY("", str); fp:=0 END;
IF Low.IsInfinity(ABS(real)) THEN Str.Append("Infinity", str); RETURN END;
(* extract exponent *)
exp:=Low.exponent10(real);
(* format number *)
INC(fp, RC.SigFigs-exp-2);
len:=RC.LengthFixedReal(real, fp);
IF cap>=len THEN
RealToFixed(real, fp, str);
(* pad with remaining zeros *)
IF fp<0 THEN Str.Append(".", str); INC(len) END; (* add decimal point *)
WHILE len<cap DO Str.Append("0", str); INC(len) END
ELSE
fp:=RC.LengthFloatReal(real, RC.SigFigs); (* check actual length *)
IF fp<=cap THEN
RealToFloat(real, RC.SigFigs, str);
(* pad with remaining zeros *)
Str.FindNext("E", str, 2, found, pos);
WHILE fp<cap DO Str.Insert("0", pos, str); INC(fp) END
ELSE fp:=RC.SigFigs-fp+cap;
IF fp<1 THEN fp:=1 END;
RealToFloat(real, fp, str)
END
END
END RealToStr;
END oocLRealStr.

101
src/lib/ooc/oocLongInts.Mod Normal file
View file

@ -0,0 +1,101 @@
(* $Id: LongInts.Mod,v 1.3 1999/09/02 13:14:52 acken Exp $ *)
MODULE oocLongInts;
(*
LongInts - Simple extended integer implementation.
Copyright (C) 1996 Michael Griebling
This module is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.
This module is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*)
CONST
B*=8000H;
TYPE
LongInt*=ARRAY 170 OF INTEGER;
PROCEDURE MinDigit * (VAR w: LongInt) : LONGINT;
VAR min, l: LONGINT;
BEGIN
min:=1; l:=LEN(w)-1;
WHILE (min<l) & (w[min]=0) DO INC(min) END;
RETURN min
END MinDigit;
PROCEDURE MultDigit * (VAR w: LongInt; digit, k: LONGINT);
VAR i, t, min: LONGINT;
BEGIN
i:=LEN(w)-1; min:=MinDigit(w)-2;
REPEAT
t:=w[i]*digit+k; (* multiply *)
w[i]:=SHORT(t MOD B); k:=t DIV B; (* generate result & carry *)
DEC(i)
UNTIL i=min
END MultDigit;
PROCEDURE AddDigit * (VAR w: LongInt; k: LONGINT);
VAR i, t, min: LONGINT;
BEGIN
i:=LEN(w)-1; min:=MinDigit(w)-2;
REPEAT
t:=w[i]+k; (* add *)
w[i]:=SHORT(t MOD B); k:=t DIV B; (* generate result & carry *)
DEC(i)
UNTIL i=min
END AddDigit;
PROCEDURE DivDigit * (VAR w: LongInt; digit: LONGINT; VAR r: LONGINT);
VAR j, t, m: LONGINT;
BEGIN
j:=MinDigit(w)-1; r:=0; m:=LEN(w)-1;
REPEAT
t:=r*B+w[j];
w[j]:=SHORT(t DIV digit); r:=t MOD digit; (* generate result & remainder *)
INC(j)
UNTIL j>m
END DivDigit;
PROCEDURE TenPower * (VAR x: LongInt; power: INTEGER);
VAR exp, i: INTEGER; d: LONGINT;
BEGIN
IF power>0 THEN
exp:=power DIV 4; power:=power MOD 4;
FOR i:=1 TO exp DO MultDigit(x, 10000, 0) END;
FOR i:=1 TO power DO MultDigit(x, 10, 0) END
ELSIF power<0 THEN
power:=-power;
exp:=power DIV 4; power:=power MOD 4;
FOR i:=1 TO exp DO DivDigit(x, 10000, d) END;
FOR i:=1 TO power DO DivDigit(x, 10, d) END
END
END TenPower;
PROCEDURE BPower * (VAR x: LongInt; power: INTEGER);
VAR i, lx: LONGINT;
BEGIN
lx:=LEN(x);
IF power>0 THEN
FOR i:=1 TO lx-1-power DO x[i]:=x[i+power] END;
FOR i:=lx-power TO lx-1 DO x[i]:=0 END
ELSIF power<0 THEN
power:=-power;
FOR i:=lx-1-power TO 1 BY -1 DO x[i+power]:=x[i] END;
FOR i:=1 TO power DO x[i]:=0 END
END
END BPower;
END oocLongInts.