diff --git a/makefile b/makefile index 9ee8fb5b..475d19ab 100644 --- a/makefile +++ b/makefile @@ -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 diff --git a/makefile.gnuc.armv6j b/makefile.gnuc.armv6j index 71d012b9..d88f9695 100644 --- a/makefile.gnuc.armv6j +++ b/makefile.gnuc.armv6j @@ -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 diff --git a/makefile.gnuc.armv6j_hardfp b/makefile.gnuc.armv6j_hardfp index f843a2c3..c98d3dc0 100644 --- a/makefile.gnuc.armv6j_hardfp +++ b/makefile.gnuc.armv6j_hardfp @@ -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 diff --git a/makefile.gnuc.armv7a_hardfp b/makefile.gnuc.armv7a_hardfp index 12182320..af52b661 100644 --- a/makefile.gnuc.armv7a_hardfp +++ b/makefile.gnuc.armv7a_hardfp @@ -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 diff --git a/makefile.gnuc.x86 b/makefile.gnuc.x86 index 5d1e812e..6203ec3a 100644 --- a/makefile.gnuc.x86 +++ b/makefile.gnuc.x86 @@ -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 diff --git a/makefile.gnuc.x86_64 b/makefile.gnuc.x86_64 index 9ee8fb5b..475d19ab 100644 --- a/makefile.gnuc.x86_64 +++ b/makefile.gnuc.x86_64 @@ -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 diff --git a/src/lib/ooc/oocLRealConv.Mod b/src/lib/ooc/oocLRealConv.Mod new file mode 100644 index 00000000..a596e6de --- /dev/null +++ b/src/lib/ooc/oocLRealConv.Mod @@ -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]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 (max1 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 real10 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=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 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)= 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 xLnInfinity THEN l.ErrorHandler(Overflow); RETURN huge + ELSIF 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 xONE 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)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 iw1ZERO 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 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 (x0 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)=l.expoMax-3 THEN l.ErrorHandler(Overflow); atan:=PiOver2 + ELSIF xnExp-xdExpABS(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)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 xdLnInfinity 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 fSqrtInfinity*HALF THEN l.ErrorHandler(HypInvTrigClipped); + IF x>ZERO THEN RETURN ln(SqrtInfinity) ELSE RETURN -ln(SqrtInfinity) END; + ELSIF xSqrtInfinity*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 xTanhLimit 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. + diff --git a/src/lib/ooc/oocLRealStr.Mod b/src/lib/ooc/oocLRealStr.Mod new file mode 100644 index 00000000..42cdc26e --- /dev/null +++ b/src/lib/ooc/oocLRealStr.Mod @@ -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=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 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 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=len THEN + RealToFixed(real, fp, str); + + (* pad with remaining zeros *) + IF fp<0 THEN Str.Append(".", str); INC(len) END; (* add decimal point *) + WHILE lenm +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.