diff --git a/makefile b/makefile index b862560d..9ee8fb5b 100644 --- a/makefile +++ b/makefile @@ -133,6 +133,7 @@ stage6: $(VOCSTATIC) -sP ooc2IntStr.Mod $(VOCSTATIC) -sP ooc2Real0.Mod $(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod + $(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.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 1a1b9d70..71d012b9 100644 --- a/makefile.gnuc.armv6j +++ b/makefile.gnuc.armv6j @@ -133,6 +133,7 @@ stage6: $(VOCSTATIC) -sP ooc2IntStr.Mod $(VOCSTATIC) -sP ooc2Real0.Mod $(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod + $(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.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 7e9ce7dc..f843a2c3 100644 --- a/makefile.gnuc.armv6j_hardfp +++ b/makefile.gnuc.armv6j_hardfp @@ -133,6 +133,7 @@ stage6: $(VOCSTATIC) -sP ooc2IntStr.Mod $(VOCSTATIC) -sP ooc2Real0.Mod $(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod + $(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.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 2a7b60bd..12182320 100644 --- a/makefile.gnuc.armv7a_hardfp +++ b/makefile.gnuc.armv7a_hardfp @@ -133,6 +133,7 @@ stage6: $(VOCSTATIC) -sP ooc2IntStr.Mod $(VOCSTATIC) -sP ooc2Real0.Mod $(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod + $(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.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 c6b7ebf8..5d1e812e 100644 --- a/makefile.gnuc.x86 +++ b/makefile.gnuc.x86 @@ -133,6 +133,7 @@ stage6: $(VOCSTATIC) -sP ooc2IntStr.Mod $(VOCSTATIC) -sP ooc2Real0.Mod $(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod + $(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.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 b862560d..9ee8fb5b 100644 --- a/makefile.gnuc.x86_64 +++ b/makefile.gnuc.x86_64 @@ -133,6 +133,7 @@ stage6: $(VOCSTATIC) -sP ooc2IntStr.Mod $(VOCSTATIC) -sP ooc2Real0.Mod $(VOCSTATIC) -sP oocLowReal.Mod oocLowLReal.Mod + $(VOCSTATIC) -sP oocRealMath.Mod oocOakMath.Mod $(VOCSTATIC) -sP oocwrapperlibc.Mod $(VOCSTATIC) -sP ulmSYSTEM.Mod $(VOCSTATIC) -sP ulmASCII.Mod ulmSets.Mod diff --git a/ocat b/ocat index 443db97c..85a685c8 100755 Binary files a/ocat and b/ocat differ diff --git a/showdef.REMOVED.git-id b/showdef.REMOVED.git-id index 8def490d..ce9c3c1d 100644 --- a/showdef.REMOVED.git-id +++ b/showdef.REMOVED.git-id @@ -1 +1 @@ -fdd63e86437d5144b706952e0eacff1129b2fa9f \ No newline at end of file +fbc7697598433a1043bcaa070da8870e6f7f07c2 \ No newline at end of file diff --git a/src/lib/ooc/oocOakMath.Mod b/src/lib/ooc/oocOakMath.Mod new file mode 100644 index 00000000..fdfb2818 --- /dev/null +++ b/src/lib/ooc/oocOakMath.Mod @@ -0,0 +1,137 @@ +(* $Id: OakMath.Mod,v 1.1 1997/02/07 07:45:32 oberon1 Exp $ *) +MODULE oocOakMath; + +IMPORT RealMath := oocRealMath; + + +CONST + pi* = RealMath.pi; + e* = RealMath.exp1; + +PROCEDURE sqrt* (x: REAL): REAL; +(* sqrt(x) returns the square root of x, where x must be positive. *) + BEGIN + RETURN RealMath.sqrt (x) + END sqrt; + +PROCEDURE power* (x, base: REAL): REAL; +(* power(x, base) returns the x to the power base. *) + BEGIN + RETURN RealMath.power (x, base) + END power; + +PROCEDURE exp* (x: REAL): REAL; +(* exp(x) is the exponential of x base e. x must not be so small that this + exponential underflows nor so large that it overflows. *) + BEGIN + RETURN RealMath.exp (x) + END exp; + +PROCEDURE ln* (x: REAL): REAL; +(* ln(x) returns the natural logarithm (base e) of x. *) + BEGIN + RETURN RealMath.ln (x) + END ln; + +PROCEDURE log* (x, base: REAL): REAL; +(* log(x,base) is the logarithm of x base b. All positive arguments are + allowed. The base b must be positive. *) + BEGIN + RETURN RealMath.log (x, base) + END log; + +PROCEDURE round* (x: REAL): REAL; +(* round(x) if fraction part of x is in range 0.0 to 0.5 then the result is + the largest integer not greater than x, otherwise the result is x rounded + up to the next highest whole number. Note that integer values cannot always + be exactly represented in REAL or REAL format. *) + BEGIN + RETURN RealMath.round (x) + END round; + +PROCEDURE sin* (x: REAL): REAL; + BEGIN + RETURN RealMath.sin (x) + END sin; + +PROCEDURE cos* (x: REAL): REAL; + BEGIN + RETURN RealMath.cos (x) + END cos; + +PROCEDURE tan* (x: REAL): REAL; +(* sin, cos, tan(x) returns the sine, cosine or tangent value of x, where x is + in radians. *) + BEGIN + RETURN RealMath.tan (x) + END tan; + +PROCEDURE arcsin* (x: REAL): REAL; + BEGIN + RETURN RealMath.arcsin (x) + END arcsin; + +PROCEDURE arccos* (x: REAL): REAL; + BEGIN + RETURN RealMath.arccos (x) + END arccos; + +PROCEDURE arctan* (x: REAL): REAL; +(* arcsin, arcos, arctan(x) returns the arcsine, arcos, arctan value in radians + of x, where x is in the sine, cosine or tangent value. *) + BEGIN + RETURN RealMath.arctan (x) + END arctan; + +PROCEDURE arctan2* (xn, xd: REAL): REAL; +(* 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. *) + BEGIN + RETURN RealMath.arctan2 (xn, xd) + END arctan2; + + +PROCEDURE sinh* (x: REAL): REAL; +(* sinh(x) is the hyperbolic sine of x. The argument x must not be so large + that exp(|x|) overflows. *) + BEGIN + RETURN RealMath.sinh (x) + END sinh; + +PROCEDURE cosh* (x: REAL): REAL; +(* cosh(x) is the hyperbolic cosine of x. The argument x must not be so large + that exp(|x|) overflows. *) + BEGIN + RETURN RealMath.cosh (x) + END cosh; + +PROCEDURE tanh* (x: REAL): REAL; +(* tanh(x) is the hyperbolic tangent of x. All arguments are legal. *) + BEGIN + RETURN RealMath.tanh (x) + END tanh; + +PROCEDURE arcsinh* (x: REAL): REAL; +(* arcsinh(x) is the arc hyperbolic sine of x. All arguments are legal. *) + BEGIN + RETURN RealMath.arcsinh (x) + END arcsinh; + +PROCEDURE arccosh* (x: REAL): REAL; +(* arccosh(x) is the arc hyperbolic cosine of x. All arguments greater than + or equal to 1 are legal. *) + BEGIN + RETURN RealMath.arccosh (x) + END arccosh; + +PROCEDURE arctanh* (x: REAL): REAL; +(* 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. *) + BEGIN + RETURN RealMath.arctanh (x) + END arctanh; + + +END oocOakMath. diff --git a/src/lib/ooc/oocRealMath.Mod b/src/lib/ooc/oocRealMath.Mod new file mode 100644 index 00000000..611f1f96 --- /dev/null +++ b/src/lib/ooc/oocRealMath.Mod @@ -0,0 +1,609 @@ +(* $Id: RealMath.Mod,v 1.6 1999/09/02 13:19:17 acken Exp $ *) +MODULE oocRealMath; + +(* + RealMath - Target independent mathematical functions for REAL + (IEEE single-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) 1995 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 := oocLowReal, S := SYSTEM; + +CONST + pi* = 3.1415926535897932384626433832795028841972; + exp1* = 2.7182818284590452353602874713526624977572; + + ZERO=0.0; ONE=1.0; HALF=0.5; TWO=2.0; (* local constants *) + + (* internally-used constants *) + huge=l.large; (* largest number this package accepts *) + miny=ONE/huge; (* smallest number this package accepts *) + sqrtHalf=0.70710678118654752440; + Limit=2.4414062E-4; (* 2**(-MantBits/2) *) + eps=2.9802322E-8; (* 2**(-MantBits-1) *) + piInv=0.31830988618379067154; (* 1/pi *) + piByTwo=1.57079632679489661923132; + piByFour=0.78539816339744830962; + lnv=0.6931610107421875; (* should be exact *) + vbytwo=0.13830277879601902638E-4; (* used in sinh/cosh *) + ln2Inv=1.44269504088896340735992468100189213; + + (* error/exception codes *) + NoError*=0; IllegalRoot*=1; IllegalLog*=2; Overflow*=3; IllegalPower*=4; IllegalLogBase*=5; + IllegalTrig*=6; IllegalInvTrig*=7; HypInvTrigClipped*=8; IllegalHypInvTrig*=9; + LossOfAccuracy*=10; Underflow*=11; + +VAR + a1: ARRAY 18 OF REAL; (* lookup table for power function *) + a2: ARRAY 9 OF REAL; (* lookup table for power function *) + em: REAL; (* largest number such that 1+epsilon > 1.0 *) + LnInfinity: REAL; (* natural log of infinity *) + LnSmall: REAL; (* natural log of very small number *) + SqrtInfinity: REAL; (* square root of infinity *) + TanhMax: REAL; (* maximum Tanh value *) + t: REAL; (* internal variables *) + +(* internally used support routines *) + +PROCEDURE SinCos (x, y, sign: REAL): REAL; + CONST + ymax=9099; (* ENTIER(pi*2**(MantBits/2)) *) + r1=-0.1666665668E+0; + r2= 0.8333025139E-2; + r3=-0.1980741872E-3; + r4= 0.2601903036E-5; + VAR + n: LONGINT; xn, f, g: REAL; +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 *) + f:=SHORT(ABS(LONG(x)) - LONG(xn)*pi); + + (* Pre: |f| <= pi/2 *) + IF ABS(f)= 0 *) + CONST + P0=0.41731; P1=0.59016; + VAR + xMant, yEst, z: REAL; xExp: INTEGER; +BEGIN + (* optimize zeros and check for illegal negative roots *) + IF x=ZERO THEN RETURN ZERO END; + IF x=LnInfinity THEN l.ErrorHandler(Overflow); RETURN huge + ELSIF x 0 *) + CONST + c1=355.0/512.0; c2=-2.121944400546905827679E-4; + A0=-0.5527074855E+0; B0=-0.6632718214E+1; + VAR f, zn, zd, r, z, w, xn: REAL; 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; r:=z+z*(w*A0/(w+B0)); + + (* 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: REAL): REAL; + (* Returns the sine of x for all x *) +BEGIN + IF xymax THEN l.ErrorHandler(LossOfAccuracy); RETURN ZERO END; + + (* determine n and the fraction f *) + n:=round(x*twoByPi); xn:=n; + f:=SHORT(LONG(x)-LONG(xn)*piByTwo); + + (* check for underflow *) + IF ABS(f)HALF THEN + i:=1-flag; + IF y>ONE THEN l.ErrorHandler(IllegalInvTrig); res:=huge; RETURN END; + + (* reduce the input argument *) + g:=(ONE-y)*HALF; r:=-sqrt(g); y:=r+r; + + (* compute approximation *) + r:=((P2*g+P1)*g)/((g+Q1)*g+Q0); + res:=y+(y*r) + ELSE + i:=flag; + IF yONE THEN f:=ONE/f; n:=2 + ELSE n:=0 + END; + + (* check if f should be scaled *) + IF f>rt32 THEN f:=(((a*f-HALF)-HALF)+f)/(rt3+f); INC(n) END; + + (* check for underflow *) + IF ABS(f)1 THEN res:=-res END; + CASE n OF + | 1: res:=res+piBySix + | 2: res:=res+piByTwo + | 3: res:=res+piByThree + | ELSE (* do nothing *) + END; + RETURN res +END atan; + +PROCEDURE arctan*(x: REAL): REAL; + (* Returns the arctangent of x, in the range [-pi/2, pi/2] for all x *) +BEGIN + IF x<0 THEN RETURN -atan(-x) + ELSE RETURN atan(x) + END +END arctan; + +PROCEDURE power*(base, exponent: REAL): REAL; + (* Returns the value of the number base raised to the power exponent + for base > 0 *) + CONST P1=0.83357541E-1; K=0.4426950409; + Q1=0.69314675; Q2=0.24018510; Q3=0.54360383E-1; + OneOver16=0.0625; XMAX=16*(l.expoMax+1)-1; (*XMIN=16*l.expoMin;*) XMIN=-2016; (* to make it easier for voc; -- noch *) + VAR z, g, R, v, u2, u1, w1, w2: REAL; w: LONGREAL; + m, p, i: INTEGER; mp, pp, iw1: LONGINT; +BEGIN + (* handle all possible error conditions *) + IF base<=ZERO THEN + IF base#ZERO THEN l.ErrorHandler(IllegalPower); base:=-base + ELSIF exponent>ZERO THEN RETURN ZERO + ELSE l.ErrorHandler(IllegalPower); 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:=P1*v*z; R:=R+K*R; u2:=(R+z*K)+z; + u1:=(m*16-p)*OneOver16; w:=LONG(exponent)*(LONG(u1)+LONG(u2)); (* need extra precision *) + + (* calculations below were modified to work properly -- incorrect in cited reference? *) + iw1:=ENTIER(16*w); w1:=iw1*OneOver16; w2:=SHORT(w-w1); + + (* 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:=((Q3*w2+Q2)*w2+Q1)*w2; z:=a1[pp+1]+a1[pp+1]*z; + RETURN l.scale(z, SHORT(mp)) +END power; + +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: REAL): REAL; +(* log(x,base) is the logarithm of x base 'base'. All positive arguments are + allowed but base > 0 and base # 1 *) +BEGIN + (* log(x, base) = ln(x) / ln(base) *) + IF base<=ZERO THEN l.ErrorHandler(IllegalLogBase); RETURN -huge + ELSE RETURN ln(x)/ln(base) + END +END log; + +PROCEDURE ipower* (x: REAL; base: INTEGER): REAL; +(* ipower(x, base) returns the x to the integer power base where Log2(x) < expoMax *) + VAR Exp: INTEGER; y: REAL; neg: BOOLEAN; + + PROCEDURE Adjust(xadj: REAL): REAL; + 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: REAL; VAR Sin, Cos: REAL); +(* 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: REAL): REAL; +(* 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. *) +VAR + res: REAL; xpdiff: LONGINT; +BEGIN + (* check for error conditions *) + IF xd=ZERO THEN + IF xn=ZERO THEN l.ErrorHandler(IllegalTrig); RETURN ZERO + ELSIF xn<0 THEN RETURN -piByTwo + ELSE RETURN piByTwo + END; + ELSE + xpdiff:=l.exponent(xn)-l.exponent(xd); + IF ABS(xpdiff)>=l.expoMax-3 THEN + (* overflow detected *) + IF xn<0 THEN RETURN -piByTwo + ELSE RETURN piByTwo + END + ELSE + res:=ABS(xn/xd); + IF res#ZERO THEN res:=atan(res) END; + 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: REAL): REAL; +(* cosh(x) is the hyperbolic cosine of x. The argument x must not be so large + that exp(|x|) overflows. *) + VAR y, f: REAL; +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: REAL): REAL; +(* tanh(x) is the hyperbolic tangent of x. All arguments are legal. *) + CONST P0=-0.8237728127; P1=-0.3831010665E-2; Q0=2.471319654; ln3over2=0.5493061443; + BIG=9.010913347; (* (ln(2)+(t+1)*ln(B))/2 where t=mantissa bits, B=base *) + VAR f, t: REAL; +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: REAL): REAL; +(* 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.999984991; (* Tanh(5.9) *) + VAR t: REAL; +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; + +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) *) + a1[1] :=ONE; + a1[2] :=S.VAL(REAL, 3F75257DH); + a1[3] :=S.VAL(REAL, 3F6AC0C7H); + a1[4] :=S.VAL(REAL, 3F60CCDFH); + a1[5] :=S.VAL(REAL, 3F5744FDH); + a1[6] :=S.VAL(REAL, 3F4E248CH); + a1[7] :=S.VAL(REAL, 3F45672AH); + a1[8] :=S.VAL(REAL, 3F3D08A4H); + a1[9] :=S.VAL(REAL, 3F3504F3H); + a1[10]:=S.VAL(REAL, 3F2D583FH); + a1[11]:=S.VAL(REAL, 3F25FED7H); + a1[12]:=S.VAL(REAL, 3F1EF532H); + a1[13]:=S.VAL(REAL, 3F1837F0H); + a1[14]:=S.VAL(REAL, 3F11C3D3H); + a1[15]:=S.VAL(REAL, 3F0B95C2H); + a1[16]:=S.VAL(REAL, 3F05AAC3H); + a1[17]:=HALF; + + (* a2[i]=2**[(1-2i)/16] - a1[2i]; delta resolution *) + a2[1]:=S.VAL(REAL, 31A92436H); + a2[2]:=S.VAL(REAL, 336C2A95H); + a2[3]:=S.VAL(REAL, 31A8FC24H); + a2[4]:=S.VAL(REAL, 331F580CH); + a2[5]:=S.VAL(REAL, 336A42A1H); + a2[6]:=S.VAL(REAL, 32C12342H); + a2[7]:=S.VAL(REAL, 32E75624H); + a2[8]:=S.VAL(REAL, 32CF9890H) +END oocRealMath. diff --git a/voc.REMOVED.git-id b/voc.REMOVED.git-id index da3d00ff..a494cc2b 100644 --- a/voc.REMOVED.git-id +++ b/voc.REMOVED.git-id @@ -1 +1 @@ -fa4f569630633bf18b885d2dba015044a9039f7a \ No newline at end of file +eb7855c3da28b25bac66fa2b1632c6523dc6b801 \ No newline at end of file diff --git a/vocstatic.REMOVED.git-id b/vocstatic.REMOVED.git-id index da3d00ff..a494cc2b 100644 --- a/vocstatic.REMOVED.git-id +++ b/vocstatic.REMOVED.git-id @@ -1 +1 @@ -fa4f569630633bf18b885d2dba015044a9039f7a \ No newline at end of file +eb7855c3da28b25bac66fa2b1632c6523dc6b801 \ No newline at end of file