/* LibMImplA.c */
/* NFS 25-Nov-85 11:43:58 */
/* log(), log1p(), exp(), expm1(), frexp() */
asm (" export Libm");
asm (" export LibmSupport");
mesa double LibmSupport←copysign(), LibmSupport←scalb();
mesa double LibmSupport←drem(), LibmSupport←logb();
mesa int LibmSupport←finite();
#define copysign(x,y) (LibmSupport←copysign(x,y))
#define scalb(x,n) (LibmSupport←scalb(x,n))
#define logb(x) (LibmSupport←logb(x))
#define finite(x) (LibmSupport←finite(x))
#define drem(x,p) (LibmSupport←drem(x,p))
/* LOG(X)
* RETURN THE LOGARITHM OF x
* DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
* CODED IN C BY K.C. NG, 1/19/85;
* REVISED BY K.C. NG on 2/7/85, 3/7/85, 3/24/85, 4/16/85.
*
* Required system supported functions:
* scalb(x,n)
* copysign(x,y)
* logb(x)
* finite(x)
*
* Required kernel function:
* log←←L(z)
*
* Method :
* 1. Argument Reduction: find k and f such that
* x = 2↑k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* log(1+f) is computed by
*
* log(1+f) = 2s + s*log←←L(s*s)
* where
* log←←L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
*
* See log←←L() for the values of the coefficients.
*
* 3. Finally, log(x) = k*ln2 + log(1+f). (Here n*ln2 will be stored
* in two floating point number: n*ln2hi + n*ln2lo, n*ln2hi is exact
* since the last 20 bits of ln2hi is 0.)
*
* Special cases:
* log(x) is NaN with signal if x < 0 (including -INF) ;
* log(+INF) is +INF; log(0) is -INF with signal;
* log(NaN) is that NaN with no signal.
*
* Accuracy:
* log(x) returns the exact log(x) nearly rounded. In a test run with
* 1,536,000 random arguments on a VAX, the maximum observed error was
* .826 ulps (units in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
double static
ln2hi = 6.9314718036912381649E-1 , /*Hex 2↑ -1 * 1.62E42FEE00000 */
ln2lo = 1.9082149292705877000E-10 , /*Hex 2↑-33 * 1.A39EF35793C76 */
sqrt2 = 1.4142135623730951455E0 ; /*Hex 2↑ 0 * 1.6A09E667F3BCD */
asm (" exportproc ←log, Libm");
double log(x)
double x;
{
static double zero=0.0, negone= -1.0, half=1.0/2.0;
double LibmSupport←logb(),
LibmSupport←scalb(),
LibmSupport←copysign(),
log←←L(),s,z,t;
int k,n,finite();
if(x!=x) return(x); /* x is NaN */
if(finite(x)) {
if( x > zero ) {
/* argument reduction */
k=logb(x); x=scalb(x,-k);
if(k == -1022) /* subnormal no. */
{n=logb(x); x=scalb(x,-n); k+=n;}
if(x >= sqrt2 ) {k += 1; x *= half;}
x += negone ;
/* compute log(1+x) */
s=x/(2+x); t=x*x*half;
z=k*ln2lo+s*(t+log←←L(s*s));
x += (z - t) ;
return(k*ln2hi+x);
}
/* end of if (x > zero) */
else {
/* zero argument, return -INF with signal */
if ( x == zero )
return( negone/zero );
/* negative argument, return NaN with signal */
else
return ( zero / zero );
}
}
/* end of if (finite(x)) */
/* NOT REACHED ifdef VAX */
/* log(-INF) is NaN with signal */
else if (x<0)
return(zero/zero);
/* log(+INF) is +INF */
else return(x);
}
/* LOG1P(x)
* RETURN THE LOGARITHM OF 1+x
* DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS)
* CODED IN C BY K.C. NG, 1/19/85;
* REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85.
*
* Required system supported functions:
* scalb(x,n)
* copysign(x,y)
* logb(x)
* finite(x)
*
* Required kernel function:
* log←←L(z)
*
* Method :
* 1. Argument Reduction: find k and f such that
* 1+x = 2↑k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* log(1+f) is computed by
*
* log(1+f) = 2s + s*log←←L(s*s)
* where
* log←←L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
*
* See log←←L() for the values of the coefficients.
*
* 3. Finally, log(1+x) = k*ln2 + log(1+f).
*
* Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers
* n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last
* 20 bits (for VAX D format), or the last 21 bits ( for IEEE
* double) is 0. This ensures n*ln2hi is exactly representable.
* 2. In step 1, f may not be representable. A correction term c
* for f is computed. It follows that the correction term for
* f - t (the leading term of log(1+f) in step 2) is c-c*x. We
* add this correction term to n*ln2lo to attenuate the error.
*
*
* Special cases:
* log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal;
* log1p(INF) is +INF; log1p(-1) is -INF with signal;
* only log1p(0)=0 is exact for finite argument.
*
* Accuracy:
* log1p(x) returns the exact log(1+x) nearly rounded. In a test run
* with 1,536,000 random arguments on a VAX, the maximum observed
* error was .846 ulps (units in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
asm (" exportproc ←log1p, Libm");
double log1p(x)
double x;
{
static double zero=0.0, negone= -1.0, one=1.0,
half=1.0/2.0, small=1.0E-20; /* 1+small == 1 */
double LibmSupport←logb(),
LibmSupport←copysign(),
LibmSupport←scalb(), log←←L(),z,s,t,c;
int k,finite();
if(x!=x) return(x); /* x is NaN */
if(finite(x)) {
if( x > negone ) {
/* argument reduction */
if(copysign(x,one)<small) return(x);
k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k);
if(z+t >= sqrt2 )
{ k += 1 ; z *= half; t *= half; }
t += negone; x = z + t;
c = (t-x)+z ; /* correction term for x */
/* compute log(1+x) */
s = x/(2+x); t = x*x*half;
c += (k*ln2lo-c*x);
z = c+s*(t+log←←L(s*s));
x += (z - t) ;
return(k*ln2hi+x);
}
/* end of if (x > negone) */
else {
/* x = -1, return -INF with signal */
if ( x == negone ) return( negone/zero );
/* negative argument for log, return NaN with signal */
else return ( zero / zero );
}
}
/* end of if (finite(x)) */
/* log(-INF) is NaN */
else if(x<0)
return(zero/zero);
/* log(+INF) is INF */
else return(x);
}
/* LOG10(X)
* RETURN THE BASE 10 LOGARITHM OF x
* DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
* CODED IN C BY K.C. NG, 1/20/85;
* REVISED BY K.C. NG on 1/23/85, 3/7/85, 4/16/85.
*
* Required kernel function:
* log(x)
*
* Method :
* log(x)
* log10(x) = --------- or [1/log(10)]*log(x)
* log(10)
*
* Note:
* [log(10)] rounded to 56 bits has error .0895 ulps,
* [1/log(10)] rounded to 53 bits has error .198 ulps;
* therefore, for better accuracy, in VAX D format, we divide
* log(x) by log(10), but in IEEE Double format, we multiply
* log(x) by [1/log(10)].
*
* Special cases:
* log10(x) is NaN with signal if x < 0;
* log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
* log10(NaN) is that NaN with no signal.
*
* Accuracy:
* log10(X) returns the exact log10(x) nearly rounded. In a test run
* with 1,536,000 random arguments on a VAX, the maximum observed
* error was 1.74 ulps (units in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
static double
ivln10 = 4.3429448190325181667E-1 ; /*Hex 2↑ -2 * 1.BCB7B1526E50E */
asm (" exportproc ←log10, Libm");
double log10(x)
double x;
{
double log();
return(ivln10*log(x));
}
/* log←←L(Z)
* LOG(1+X) - 2S X
* RETURN --------------- WHERE Z = S*S, S = ------- , 0 <= Z <= .0294...
* S 2 + X
*
* DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
* KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS
* CODED IN C BY K.C. NG, 1/19/85;
* REVISED BY K.C. Ng, 2/3/85, 4/16/85.
*
* Method :
* 1. Polynomial approximation: let s = x/(2+x).
* Based on log(1+x) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
*
* (log(1+x) - 2s)/s is computed by
*
* z*(L1 + z*(L2 + z*(... (L7 + z*L8)...)))
*
* where z=s*s. (See the listing below for Lk's values.) The
* coefficients are obtained by a special Remez algorithm.
*
* Accuracy:
* Assuming no rounding error, the maximum magnitude of the approximation
* error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63)
* for VAX D format.
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
static double
L1 = 6.6666666666667340202E-1 , /*Hex 2↑ -1 * 1.5555555555592 */
L2 = 3.9999999999416702146E-1 , /*Hex 2↑ -2 * 1.999999997FF24 */
L3 = 2.8571428742008753154E-1 , /*Hex 2↑ -2 * 1.24924941E07B4 */
L4 = 2.2222198607186277597E-1 , /*Hex 2↑ -3 * 1.C71C52150BEA6 */
L5 = 1.8183562745289935658E-1 , /*Hex 2↑ -3 * 1.74663CC94342F */
L6 = 1.5314087275331442206E-1 , /*Hex 2↑ -3 * 1.39A1EC014045B */
L7 = 1.4795612545334174692E-1 ; /*Hex 2↑ -3 * 1.2F039F0085122 */
double log←←L(z)
double z;
{
return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7)))))));
}
/* EXP(X)
* RETURN THE EXPONENTIAL OF X
* DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
* CODED IN C BY K.C. NG, 1/19/85;
* REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85.
*
* Required system supported functions:
* scalb(x,n)
* copysign(x,y)
* finite(x)
*
* Kernel function:
* exp←←E(x,c)
*
* Method:
* 1. Argument Reduction: given the input x, find r and integer k such
* that
* x = k*ln2 + r, |r| <= 0.5*ln2 .
* r will be represented as r := z+c for better accuracy.
*
* 2. Compute expm1(r)=exp(r)-1 by
*
* expm1(r=z+c) := z + exp←←E(z,r)
*
* 3. exp(x) = 2↑k * ( expm1(r) + 1 ).
*
* Special cases:
* exp(INF) is INF, exp(NaN) is NaN;
* exp(-INF)= 0;
* for finite argument, only exp(0)=1 is exact.
*
* Accuracy:
* exp(x) returns the exponential of x nearly rounded. In a test run
* with 1,156,000 random arguments on a VAX, the maximum observed
* error was .768 ulps (units in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
double static
lnhuge = 7.1602103751842355450E2 , /*Hex 2↑ 9 * 1.6602B15B7ECF2 */
lntiny = -7.5137154372698068983E2 , /*Hex 2↑ 9 * -1.77AF8EBEAE354 */
invln2 = 1.4426950408889633870E0 ; /*Hex 2↑ 0 * 1.71547652B82FE */
asm (" exportproc ←exp, Libm");
double exp(x)
double x;
{
double LibmSupport←scalb(),
LibmSupport←copysign(),
exp←←E(), z,hi,lo,c;
int k,LibmSupport←finite();
if(x!=x) return(x); /* x is NaN */
if( x <= lnhuge ) {
if( x >= lntiny ) {
/* argument reduction : x --> x - k*ln2 */
k=invln2*x+copysign(0.5,x); /* k=NINT(x/ln2) */
/* express x-k*ln2 as z+c */
hi=x-k*ln2hi;
z=hi-(lo=k*ln2lo);
c=(hi-z)-lo;
/* return 2↑k*[expm1(x) + 1] */
z += exp←←E(z,c);
return (scalb(z+1.0,k));
}
/* end of x > lntiny */
else
/* exp(-big#) underflows to zero */
if(finite(x)) return(scalb(1.0,-5000));
/* exp(-INF) is zero */
else return(0.0);
}
/* end of x < lnhuge */
else
/* exp(INF) is INF, exp(+big#) overflows to INF */
return( finite(x) ? scalb(1.0,5000) : x);
}
/* EXPM1(X)
* RETURN THE EXPONENTIAL OF X MINUS ONE
* DOUBLE PRECISION (IEEE 53 BITS, VAX D FORMAT 56 BITS)
* CODED IN C BY K.C. NG, 1/19/85;
* REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/21/85, 4/16/85.
*
* Required system supported functions:
* scalb(x,n)
* copysign(x,y)
* finite(x)
*
* Kernel function:
* exp←←E(x,c)
*
* Method:
* 1. Argument Reduction: given the input x, find r and integer k such
* that
* x = k*ln2 + r, |r| <= 0.5*ln2 .
* r will be represented as r := z+c for better accuracy.
*
* 2. Compute EXPM1(r)=exp(r)-1 by
*
* EXPM1(r=z+c) := z + exp←←E(z,c)
*
* 3. EXPM1(x) = 2↑k * ( EXPM1(r) + 1-2↑-k ).
*
* Remarks:
* 1. When k=1 and z < -0.25, we use the following formula for
* better accuracy:
* EXPM1(x) = 2 * ( (z+0.5) + exp←←E(z,c) )
* 2. To avoid rounding error in 1-2↑-k where k is large, we use
* EXPM1(x) = 2↑k * { [z+(exp←←E(z,c)-2↑-k )] + 1 }
* when k>56.
*
* Special cases:
* EXPM1(INF) is INF, EXPM1(NaN) is NaN;
* EXPM1(-INF)= -1;
* for finite argument, only EXPM1(0)=0 is exact.
*
* Accuracy:
* EXPM1(x) returns the exact (exp(x)-1) nearly rounded. In a test run with
* 1,166,000 random arguments on a VAX, the maximum observed error was
* .872 ulps (units of the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
asm (" exportproc ←expm1, Libm");
double expm1(x)
double x;
{
double static one=1.0, half=1.0/2.0;
double LibmSupport←scalb(),
LibmSupport←copysign(), exp←←E(), z,hi,lo,c;
int k,LibmSupport←finite();
static prec=53;
if(x!=x) return(x); /* x is NaN */
if( x <= lnhuge ) {
if( x >= -40.0 ) {
/* argument reduction : x - k*ln2 */
k= invln2 *x+copysign(0.5,x); /* k=NINT(x/ln2) */
hi=x-k*ln2hi ;
z=hi-(lo=k*ln2lo);
c=(hi-z)-lo;
if(k==0) return(z+exp←←E(z,c));
if(k==1)
if(z< -0.25)
{x=z+half;x +=exp←←E(z,c); return(x+x);}
else
{z+=exp←←E(z,c); x=half+z; return(x+x);}
/* end of k=1 */
else {
if(k<=prec)
{ x=one-scalb(one,-k); z += exp←←E(z,c);}
else if(k<100)
{ x = exp←←E(z,c)-scalb(one,-k); x+=z; z=one;}
else
{ x = exp←←E(z,c)+z; z=one;}
return (scalb(x+z,k));
}
}
/* end of x > lnunfl */
else
/* expm1(-big#) rounded to -1 (inexact) */
if(finite(x))
{ ln2hi+ln2lo; return(-one);}
/* expm1(-INF) is -1 */
else return(-one);
}
/* end of x < lnhuge */
else
/* expm1(INF) is INF, expm1(+big#) overflows to INF */
return( finite(x) ? scalb(one,5000) : x);
}
/* POW(X,Y)
* RETURN X**Y
* DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
* CODED IN C BY K.C. NG, 1/8/85;
* REVISED BY K.C. NG on 7/10/85.
*
* Required system supported functions:
* scalb(x,n)
* logb(x)
* copysign(x,y)
* finite(x)
* drem(x,y)
*
* Required kernel functions:
* exp←←E(a,c) ...return exp(a+c) - 1 - a*a/2
* log←←L(x) ...return (log(1+x) - 2s)/s, s=x/(2+x)
* pow←p(x,y) ...return +(anything)**(finite non zero)
*
* Method
* 1. Compute and return log(x) in three pieces:
* log(x) = n*ln2 + hi + lo,
* where n is an integer.
* 2. Perform y*log(x) by simulating muti-precision arithmetic and
* return the answer in three pieces:
* y*log(x) = m*ln2 + hi + lo,
* where m is an integer.
* 3. Return x**y = exp(y*log(x))
* = 2↑m * ( exp(hi+lo) ).
*
* Special cases:
* (anything) ** 0 is 1 ;
* (anything) ** 1 is itself;
* (anything) ** NaN is NaN;
* NaN ** (anything except 0) is NaN;
* +-(anything > 1) ** +INF is +INF;
* +-(anything > 1) ** -INF is +0;
* +-(anything < 1) ** +INF is +0;
* +-(anything < 1) ** -INF is +INF;
* +-1 ** +-INF is NaN and signal INVALID;
* +0 ** +(anything except 0, NaN) is +0;
* -0 ** +(anything except 0, NaN, odd integer) is +0;
* +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO;
* -0 ** -(anything except 0, NaN, odd integer) is +INF with signal;
* -0 ** (odd integer) = -( +0 ** (odd integer) );
* +INF ** +(anything except 0,NaN) is +INF;
* +INF ** -(anything except 0,NaN) is +0;
* -INF ** (odd integer) = -( +INF ** (odd integer) );
* -INF ** (even integer) = ( +INF ** (even integer) );
* -INF ** -(anything except integer,NaN) is NaN with signal;
* -(x=anything) ** (k=integer) is (-1)**k * (x ** k);
* -(anything except 0) ** (non-integer) is NaN with signal;
*
* Accuracy:
* pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
* and a Zilog Z8000,
* pow(integer,integer)
* always returns the correct integer provided it is representable.
* In a test run with 100,000 random arguments with 0 < x, y < 20.0
* on a VAX, the maximum observed error was 1.79 ulps (units in the
* last place).
*
* Constants :
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
double static zero=0.0, half=1.0/2.0, one=1.0, two=2.0, negone= -1.0;
asm (" exportproc ←pow, Libm");
double pow(x,y)
double x,y;
{
double LibmSupport←drem(),pow←p(),LibmSupport←copysign(),t;
int LibmSupport←finite();
if (y==zero) return(one);
else if(y==one
||x!=x
) return( x ); /* if x is NaN or y=1 */
else if(y!=y) return( y ); /* if y is NaN */
else if(!finite(y)) /* if y is INF */
if((t=copysign(x,one))==one) return(zero/zero);
else if(t>one) return((y>zero)?y:zero);
else return((y<zero)?-y:zero);
else if(y==two) return(x*x);
else if(y==negone) return(one/x);
/* sign(x) = 1 */
else if(copysign(one,x)==one) return(pow←p(x,y));
/* sign(x)= -1 */
/* if y is an even integer */
else if ( (t=drem(y,two)) == zero) return( pow←p(-x,y) );
/* if y is an odd integer */
else if (copysign(t,one) == one) return( -pow←p(-x,y) );
/* Henceforth y is not an integer */
else if(x==zero) /* x is -0 */
return((y>zero)?-x:one/(-x));
else { /* return NaN */
return(zero/zero);
}
}
/* pow←p(x,y) return x**y for x with sign=1 and finite y */
static double pow←p(x,y)
double x,y;
{
double logb(),
LibmSupport←scalb(),
LibmSupport←copysign(),log←←L(),exp←←E();
double c,s,t,z,tx,ty;
float sx,sy;
long k=0;
int n,m;
if(x==zero||!finite(x)) { /* if x is +INF or +0 */
return((y>zero)?x:one/x);
}
if(x==1.0) return(x); /* if x=1.0, return 1 since y is finite */
/* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */
z=scalb(x,-(n=logb(x)));
if(n <= -1022) {n += (m=logb(z)); z=scalb(z,-m);}
if(z >= sqrt2 ) {n += 1; z *= half;} z -= one ;
/* log(x) = nlog2+log(1+z) ~ nlog2 + t + tx */
s=z/(two+z); c=z*z*half; tx=s*(c+log←←L(s*s));
t= z-(c-tx); tx += (z-t)-c;
/* if y*log(x) is neither too big nor too small */
if((s=logb(y)+logb(n+t)) < 12.0)
if(s>-60.0) {
/* compute y*log(x) ~ mlog2 + t + c */
s=y*(n+invln2*t);
m=s+copysign(half,s); /* m := nint(y*log(x)) */
k=y;
if((double)k==y) { /* if y is an integer */
k = m-k*n;
sx=t; tx+=(t-sx); }
else { /* if y is not an integer */
k =m;
tx+=n*ln2lo;
sx=(c=n*ln2hi)+t; tx+=(c-sx)+t; }
/* end of checking whether k==y */
sy=y; ty=y-sy; /* y ~ sy + ty */
s=(double)sx*sy-k*ln2hi; /* (sy+ty)*(sx+tx)-kln2 */
z=(tx*ty-k*ln2lo);
tx=tx*sy; ty=sx*ty;
t=ty+z; t+=tx; t+=s;
c= -((((t-s)-tx)-ty)-z);
/* return exp(y*log(x)) */
t += exp←←E(t,c); return(scalb(one+t,m));
}
/* end of if log(y*log(x)) > -60.0 */
else
/* exp(+- tiny) = 1 with inexact flag */
{ln2hi+ln2lo; return(one);}
else if(copysign(one,y)*(n+invln2*t) <zero)
/* exp(-(big#)) underflows to zero */
return(scalb(one,-5000));
else
/* exp(+(big#)) overflows to INF */
return(scalb(one, 5000));
}
/* exp←←E(x,c)
* ASSUMPTION: c << x SO THAT fl(x+c)=x.
* (c is the correction term for x)
* exp←←E RETURNS
*
* / exp(x+c) - 1 - x , 1E-19 < |x| < .3465736
* exp←←E(x,c) = |
* \ 0 , |x| < 1E-19.
*
* DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
* KERNEL FUNCTION OF EXP, EXPM1, POW FUNCTIONS
* CODED IN C BY K.C. NG, 1/31/85;
* REVISED BY K.C. NG on 3/16/85, 4/16/85.
*
* Required system supported function:
* copysign(x,y)
*
* Method:
* 1. Rational approximation. Let r=x+c.
* Based on
* 2 * sinh(r/2)
* exp(r) - 1 = ---------------------- ,
* cosh(r/2) - sinh(r/2)
* exp←←E(r) is computed using
* x*x (x/2)*W - ( Q - ( 2*P + x*P ) )
* --- + (c + x*[---------------------------------- + c ])
* 2 1 - W
* where P := p1*x↑2 + p2*x↑4,
* Q := q1*x↑2 + q2*x↑4 (for 56 bits precision, add q3*x↑6)
* W := x/2-(Q-x*P),
*
* (See the listing below for the values of p1,p2,q1,q2,q3. The poly-
* nomials P and Q may be regarded as the approximations to sinh
* and cosh :
* sinh(r/2) = r/2 + r * P , cosh(r/2) = 1 + Q . )
*
* The coefficients were obtained by a special Remez algorithm.
*
* Approximation error:
*
* | exp(x) - 1 | 2**(-57), (IEEE double)
* | ------------ - (exp←←E(x,0)+x)/x | <=
* | x | 2**(-69). (VAX D)
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
static double
p1 = 1.3887401997267371720E-2 , /*Hex 2↑ -7 * 1.C70FF8B3CC2CF */
p2 = 3.3044019718331897649E-5 , /*Hex 2↑-15 * 1.15317DF4526C4 */
q1 = 1.1110813732786649355E-1 , /*Hex 2↑ -4 * 1.C719538248597 */
q2 = 9.9176615021572857300E-4 ; /*Hex 2↑-10 * 1.03FC4CB8C98E8 */
double exp←←E(x,c)
double x,c;
{
double static zero=0.0, one=1.0, half=1.0/2.0, small=1.0E-19;
double LibmSupport←copysign(),z,p,q,xp,xh,w;
if(copysign(x,one)>small) {
z = x*x ;
p = z*( p1 +z* p2 );
q = z*( q1 +z* q2 );
xp= x*p ;
xh= x*half ;
w = xh-(q-xp) ;
p = p+p;
c += x*((xh*w-(q-(p+xp)))/(one-w)+c);
return(z*half+c);
}
/* end of |x| > small */
else {
if(x!=zero) one+small; /* raise the inexact flag */
return(copysign(zero,x));
}
}
/* expE just calls exp←←E, but is needed for exporting to mesa interface */
asm (" exportproc ←expE, LibmSupport");
double expE(x,c)
double x,c;
{
return (exp←←E(x,c));
}
/* FREXP
* the call
* x = frexp(arg,&exp);
* must return a double fp quantity x which is <1.0
* and the corresponding binary exponent "exp".
* such that
* arg = x*2↑exp
* if the argument is 0.0, return 0.0 mantissa and 0 exponent.
*/
asm (" exportproc ←frexp, Libm");
double
frexp(x,i)
double x;
int *i;
{
int neg;
int j;
j = 0;
neg = 0;
if(x<0){
x = -x;
neg = 1;
}
if(x>=1.0)
while(x>=1.0){
j = j+1;
x = x/2;
}
else if(x<0.5 && x != 0.0)
while(x<0.5){
j = j-1;
x = 2*x;
}
*i = j;
if(neg) x = -x;
return(x);
}