/* LibMImplB.c */
/* NFS */
/* cabs(), hypot(), cbrt(), floor(), ceil(), rint() */
asm (" export Libm");
mesa double LibmSupport←copysign(), LibmSupport←scalb();
mesa double LibmSupport←drem(), LibmSupport←logb(), Libm←sqrt();
mesa int LibmSupport←finite();
mesa double Libm←modf();
#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))
#define sqrt(x) (Libm←sqrt(x))
/* CABS(Z)
* RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER Z = X + iY
* DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
* CODED IN C BY K.C. NG, 11/28/84.
* REVISED BY K.C. NG, 7/12/85.
*
* Required kernel function :
* hypot(x,y)
*
* Method :
* cabs(z) = hypot(x,y) .
*/
asm (" exportproc ←cabs, Libm");
double cabs(z)
struct { double x, y;} z;
{
double hypot();
return(hypot(z.x,z.y));
}
/* HYPOT(X,Y)
* RETURN THE SQUARE ROOT OF X↑2 + Y↑2 WHERE Z=X+iY
* DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
* CODED IN C BY K.C. NG, 11/28/84;
* REVISED BY K.C. NG, 7/12/85.
*
* Required system supported functions :
* copysign(x,y)
* finite(x)
* scalb(x,N)
* sqrt(x)
*
* Method :
* 1. replace x by |x| and y by |y|, and swap x and
* y if y > x (hence x is never smaller than y).
* 2. Hypot(x,y) is computed by:
* Case I, x/y > 2
*
* y
* hypot = x + -----------------------------
* 2
* sqrt ( 1 + [x/y] ) + x/y
*
* Case II, x/y <= 2
* y
* hypot = x + --------------------------------------------------
* 2
* [x/y] - 2
* (sqrt(2)+1) + (x-y)/y + -----------------------------
* 2
* sqrt ( 1 + [x/y] ) + sqrt(2)
*
*
*
* Special cases:
* hypot(x,y) is INF if x or y is +INF or -INF; else
* hypot(x,y) is NAN if x or y is NAN.
*
* Accuracy:
* hypot(x,y) returns the sqrt(x↑2+y↑2) with error less than 1 ulps (units
* in the last place). See Kahan's "Interval Arithmetic Options in the
* Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
* 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
* code follows in comments.) In a test run with 500,000 random arguments
* on a VAX, the maximum observed error was .959 ulps.
*
* 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
r2p1hi = 2.4142135623730949234E0 , /*Hex 2↑1 * 1.3504F333F9DE6 */
r2p1lo = 1.2537167179050217666E-16 , /*Hex 2↑-53 * 1.21165F626CDD5 */
sqrt2 = 1.4142135623730951455E0 ; /*Hex 2↑ 0 * 1.6A09E667F3BCD */
asm (" exportproc ←hypot, Libm");
double hypot(x,y)
double x, y;
{
static double zero=0, one=1,
small=1.0E-18; /* fl(1+small)==1 */
static ibig=30; /* fl(1+2**(2*ibig))==1 */
double LibmSupport←copysign(),
LibmSupport←scalb(),
LibmSupport←logb(),
Libm←sqrt(),t,r;
int LibmSupport←finite(), exp;
if(finite(x))
if(finite(y))
{
x=copysign(x,one);
y=copysign(y,one);
if(y > x)
{ t=x; x=y; y=t; }
if(x == zero) return(zero);
if(y == zero) return(x);
exp= logb(x);
if(exp-(int)logb(y) > ibig )
/* raise inexact flag and return |x| */
{ one+small; return(x); }
/* start computing sqrt(x↑2 + y↑2) */
r=x-y;
if(r>y) { /* x/y > 2 */
r=x/y;
r=r+sqrt(one+r*r); }
else { /* 1 <= x/y <= 2 */
r/=y; t=r*(r+2.0);
r+=t/(sqrt2+sqrt(2.0+t));
r+=r2p1lo; r+=r2p1hi; }
r=y/r;
return(x+r);
}
else if(y==y) /* y is +-INF */
return(copysign(y,one));
else
return(y); /* y is NaN and x is finite */
else if(x==x) /* x is +-INF */
return (copysign(x,one));
else if(finite(y))
return(x); /* x is NaN, y is finite */
else if(y!=y) return(y); /* x and y is NaN */
else return(copysign(y,one)); /* y is INF */
}
/* A faster but less accurate version of cabs(x,y) */
#if 0
double hypot(x,y)
double x, y;
{
static double zero=0, one=1;
small=1.0E-18; /* fl(1+small)==1 */
static ibig=30; /* fl(1+2**(2*ibig))==1 */
double LibmSupport←copysign(),
LibmSupport←scalb(),
LibmSupport←logb(),
Libm←sqrt(),temp;
int LibmSupport←finite(), exp;
if(finite(x))
if(finite(y))
{
x=copysign(x,one);
y=copysign(y,one);
if(y > x)
{ temp=x; x=y; y=temp; }
if(x == zero) return(zero);
if(y == zero) return(x);
exp= logb(x);
x=scalb(x,-exp);
if(exp-(int)logb(y) > ibig )
/* raise inexact flag and return |x| */
{ one+small; return(scalb(x,exp)); }
else y=scalb(y,-exp);
return(scalb(sqrt(x*x+y*y),exp));
}
else if(y==y) /* y is +-INF */
return(copysign(y,one));
else
return(y); /* y is NaN and x is finite */
else if(x==x) /* x is +-INF */
return (copysign(x,one));
else if(finite(y))
return(x); /* x is NaN, y is finite */
else if(y!=y) return(y); /* x and y is NaN */
else return(copysign(y,one)); /* y is INF */
}
#endif
/* CBRT */
/* kahan's cube root (53 bits IEEE double precision)
* for IEEE machines only
* coded in C by K.C. Ng, 4/30/85
*
* Accuracy:
* better than 0.667 ulps according to an error analysis. Maximum
* error observed was 0.666 ulps in an 1,000,000 random arguments test.
*
* Warning: this code is semi machine dependent; the ordering of words in
* a floating point number must be known in advance. I assume that the
* long interger at the address of a floating point number will be the
* leading 32 bits of that floating point number (i.e., sign, exponent,
* and the 20 most significant bits).
* On a National machine, it has different ordering; therefore, this code
* must be compiled with flag -DNATIONAL.
*/
static unsigned long B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
static double
C= 19./35.,
D= -864./1225.,
E= 99./70.,
F= 45./28.,
G= 5./14.;
asm (" exportproc ←cbrt, Libm");
double cbrt(x)
double x;
{
double r,s,t=0.0,w;
unsigned long *px = (unsigned long *) &x,
*pt = (unsigned long *) &t,
mexp,sign;
int n0=0,n1=1;
mexp=px[n0]&0x7ff00000;
if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */
if(x==0.0) return(x); /* cbrt(0) is itself */
sign=px[n0]&0x80000000; /* sign= sign(x) */
px[n0] ↑= sign; /* x=|x| */
/* rough cbrt to 5 bits */
if(mexp==0) /* subnormal number */
{pt[n0]=0x43500000; /* set t= 2**54 */
t*=x; pt[n0]=pt[n0]/3+B2;
}
else
pt[n0]=px[n0]/3+B1;
/* new cbrt to 23 bits, may be implemented in single precision */
r=t*t/x;
s=C+r*t;
t*=G+F/(s+E+D/s);
/* chopped to 20 bits and make it larger than cbrt(x) */
pt[n1]=0; pt[n0]+=0x00000001;
/* one step newton iteration to 53 bits with error less than 0.667 ulps */
s=t*t; /* t*t is exact */
r=x/s;
w=t+t;
r=(r-t)/(w+r); /* r-s is exact */
t=t+t*r;
/* retore the sign bit */
pt[n0] |= sign;
return(t);
}
/* @(#)floor.c 4.2 9/11/85 */
/*
* floor and ceil-- greatest integer <= arg
* (resp least >=)
*/
double Libm←modf();
asm (" exportproc ←floor, Libm");
double
floor(d)
double d;
{
double fract;
if (d<0.0) {
d = -d;
fract = Libm←modf(d, &d);
if (fract != 0.0)
d += 1;
d = -d;
} else
Libm←modf(d, &d);
return(d);
}
asm (" exportproc ←ceil, Libm");
double
ceil(d)
double d;
{
return(-floor(-d));
}
/*
* algorithm for rint(x) in pseudo-pascal form ...
*
* real rint(x): real x;
* ... delivers integer nearest x in direction of prevailing rounding
* ... mode
* const L = (last consecutive integer)/2
* = 2**55; for VAX D
* = 2**52; for IEEE 754 Double
* real s,t;
* begin
* if x != x then return x; ... NaN
* if |x| >= L then return x; ... already an integer
* s := copysign(L,x);
* t := x + s; ... = (x+s) rounded to integer
* return t - s
* end;
*
* Note: Inexact will be signaled if x is not an integer, as is
* customary for IEEE 754. No other signal can be emitted.
*/
asm (" exportproc ←rint, Libm");
static double L = 4503599627370496.0E0; /* 2**52 */
double
rint(x)
double x;
{
double s,t,one = 1.0,LibmSupport←copysign();
if (x != x) /* NaN */
return (x);
if (copysign(x,one) >= L) /* already an integer */
return (x);
s = copysign(L,x);
t = x + s; /* x+s rounded to integer */
return (t - s);
}