diff options
Diffstat (limited to 'lib/libm/ieee')
| -rw-r--r-- | lib/libm/ieee/cabs.c | 233 | ||||
| -rw-r--r-- | lib/libm/ieee/cbrt.c | 121 | ||||
| -rw-r--r-- | lib/libm/ieee/support.c | 527 |
3 files changed, 0 insertions, 881 deletions
diff --git a/lib/libm/ieee/cabs.c b/lib/libm/ieee/cabs.c deleted file mode 100644 index 5b84530f02f5..000000000000 --- a/lib/libm/ieee/cabs.c +++ /dev/null @@ -1,233 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)cabs.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* 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. - */ -#include "mathimpl.h" - -vc(r2p1hi, 2.4142135623730950345E0 ,8279,411a,ef32,99fc, 2, .9A827999FCEF32) -vc(r2p1lo, 1.4349369327986523769E-17 ,597d,2484,754b,89b3, -55, .84597D89B3754B) -vc(sqrt2, 1.4142135623730950622E0 ,04f3,40b5,de65,33f9, 1, .B504F333F9DE65) - -ic(r2p1hi, 2.4142135623730949234E0 , 1, 1.3504F333F9DE6) -ic(r2p1lo, 1.2537167179050217666E-16 , -53, 1.21165F626CDD5) -ic(sqrt2, 1.4142135623730951455E0 , 0, 1.6A09E667F3BCD) - -#ifdef vccast -#define r2p1hi vccast(r2p1hi) -#define r2p1lo vccast(r2p1lo) -#define sqrt2 vccast(sqrt2) -#endif - -double -hypot(x,y) -double x, y; -{ - static const double zero=0, one=1, - small=1.0E-18; /* fl(1+small)==1 */ - static const ibig=30; /* fl(1+2**(2*ibig))==1 */ - double t,r; - int 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 */ -#if !defined(vax)&&!defined(tahoe) - else if(y!=y) return(y); /* x and y is NaN */ -#endif /* !defined(vax)&&!defined(tahoe) */ - else return(copysign(y,one)); /* y is INF */ -} - -/* 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) . - */ - -struct complex { double x, y; }; - -double -cabs(z) -struct complex z; -{ - return hypot(z.x,z.y); -} - -double -z_abs(z) -struct complex *z; -{ - return hypot(z->x,z->y); -} - -/* A faster but less accurate version of cabs(x,y) */ -#if 0 -double hypot(x,y) -double x, y; -{ - static const double zero=0, one=1; - small=1.0E-18; /* fl(1+small)==1 */ - static const ibig=30; /* fl(1+2**(2*ibig))==1 */ - double temp; - int 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 diff --git a/lib/libm/ieee/cbrt.c b/lib/libm/ieee/cbrt.c deleted file mode 100644 index 425a7af4a271..000000000000 --- a/lib/libm/ieee/cbrt.c +++ /dev/null @@ -1,121 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)cbrt.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* 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. - */ -#if !defined(vax)&&!defined(tahoe) - -static const unsigned long - B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ - B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ -static const double - C= 19./35., - D= -864./1225., - E= 99./70., - F= 45./28., - G= 5./14.; - -double cbrt(x) -double x; -{ - double r,s,t=0.0,w; - unsigned long *px = (unsigned long *) &x, - *pt = (unsigned long *) &t, - mexp,sign; - -#ifdef national /* ordering of words in a floating points number */ - const int n0=1,n1=0; -#else /* national */ - const int n0=0,n1=1; -#endif /* national */ - - 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-t is exact */ - t=t+t*r; - - - /* retore the sign bit */ - pt[n0] |= sign; - return(t); -} -#endif diff --git a/lib/libm/ieee/support.c b/lib/libm/ieee/support.c deleted file mode 100644 index f7bda5a3bd20..000000000000 --- a/lib/libm/ieee/support.c +++ /dev/null @@ -1,527 +0,0 @@ -/* - * Copyright (c) 1985, 1993 - * The Regents of the University of California. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#ifndef lint -static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93"; -#endif /* not lint */ - -/* - * Some IEEE standard 754 recommended functions and remainder and sqrt for - * supporting the C elementary functions. - ****************************************************************************** - * WARNING: - * These codes are developed (in double) to support the C elementary - * functions temporarily. They are not universal, and some of them are very - * slow (in particular, drem and sqrt is extremely inefficient). Each - * computer system should have its implementation of these functions using - * its own assembler. - ****************************************************************************** - * - * IEEE 754 required operations: - * drem(x,p) - * returns x REM y = x - [x/y]*y , where [x/y] is the integer - * nearest x/y; in half way case, choose the even one. - * sqrt(x) - * returns the square root of x correctly rounded according to - * the rounding mod. - * - * IEEE 754 recommended functions: - * (a) copysign(x,y) - * returns x with the sign of y. - * (b) scalb(x,N) - * returns x * (2**N), for integer values N. - * (c) logb(x) - * returns the unbiased exponent of x, a signed integer in - * double precision, except that logb(0) is -INF, logb(INF) - * is +INF, and logb(NAN) is that NAN. - * (d) finite(x) - * returns the value TRUE if -INF < x < +INF and returns - * FALSE otherwise. - * - * - * CODED IN C BY K.C. NG, 11/25/84; - * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85. - */ - -#include "mathimpl.h" - -#if defined(vax)||defined(tahoe) /* VAX D format */ -#include <errno.h> - static const unsigned short msign=0x7fff , mexp =0x7f80 ; - static const short prep1=57, gap=7, bias=129 ; - static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ; -#else /* defined(vax)||defined(tahoe) */ - static const unsigned short msign=0x7fff, mexp =0x7ff0 ; - static const short prep1=54, gap=4, bias=1023 ; - static const double novf=1.7E308, nunf=3.0E-308,zero=0.0; -#endif /* defined(vax)||defined(tahoe) */ - -double scalb(x,N) -double x; int N; -{ - int k; - -#ifdef national - unsigned short *px=(unsigned short *) &x + 3; -#else /* national */ - unsigned short *px=(unsigned short *) &x; -#endif /* national */ - - if( x == zero ) return(x); - -#if defined(vax)||defined(tahoe) - if( (k= *px & mexp ) != ~msign ) { - if (N < -260) - return(nunf*nunf); - else if (N > 260) { - return(copysign(infnan(ERANGE),x)); - } -#else /* defined(vax)||defined(tahoe) */ - if( (k= *px & mexp ) != mexp ) { - if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf); - if( k == 0 ) { - x *= scalb(1.0,(int)prep1); N -= prep1; return(scalb(x,N));} -#endif /* defined(vax)||defined(tahoe) */ - - if((k = (k>>gap)+ N) > 0 ) - if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap); - else x=novf+novf; /* overflow */ - else - if( k > -prep1 ) - /* gradual underflow */ - {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);} - else - return(nunf*nunf); - } - return(x); -} - - -double copysign(x,y) -double x,y; -{ -#ifdef national - unsigned short *px=(unsigned short *) &x+3, - *py=(unsigned short *) &y+3; -#else /* national */ - unsigned short *px=(unsigned short *) &x, - *py=(unsigned short *) &y; -#endif /* national */ - -#if defined(vax)||defined(tahoe) - if ( (*px & mexp) == 0 ) return(x); -#endif /* defined(vax)||defined(tahoe) */ - - *px = ( *px & msign ) | ( *py & ~msign ); - return(x); -} - -double logb(x) -double x; -{ - -#ifdef national - short *px=(short *) &x+3, k; -#else /* national */ - short *px=(short *) &x, k; -#endif /* national */ - -#if defined(vax)||defined(tahoe) - return (int)(((*px&mexp)>>gap)-bias); -#else /* defined(vax)||defined(tahoe) */ - if( (k= *px & mexp ) != mexp ) - if ( k != 0 ) - return ( (k>>gap) - bias ); - else if( x != zero) - return ( -1022.0 ); - else - return(-(1.0/zero)); - else if(x != x) - return(x); - else - {*px &= msign; return(x);} -#endif /* defined(vax)||defined(tahoe) */ -} - -finite(x) -double x; -{ -#if defined(vax)||defined(tahoe) - return(1); -#else /* defined(vax)||defined(tahoe) */ -#ifdef national - return( (*((short *) &x+3 ) & mexp ) != mexp ); -#else /* national */ - return( (*((short *) &x ) & mexp ) != mexp ); -#endif /* national */ -#endif /* defined(vax)||defined(tahoe) */ -} - -double drem(x,p) -double x,p; -{ - short sign; - double hp,dp,tmp; - unsigned short k; -#ifdef national - unsigned short - *px=(unsigned short *) &x +3, - *pp=(unsigned short *) &p +3, - *pd=(unsigned short *) &dp +3, - *pt=(unsigned short *) &tmp+3; -#else /* national */ - unsigned short - *px=(unsigned short *) &x , - *pp=(unsigned short *) &p , - *pd=(unsigned short *) &dp , - *pt=(unsigned short *) &tmp; -#endif /* national */ - - *pp &= msign ; - -#if defined(vax)||defined(tahoe) - if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */ -#else /* defined(vax)||defined(tahoe) */ - if( ( *px & mexp ) == mexp ) -#endif /* defined(vax)||defined(tahoe) */ - return (x-p)-(x-p); /* create nan if x is inf */ - if (p == zero) { -#if defined(vax)||defined(tahoe) - return(infnan(EDOM)); -#else /* defined(vax)||defined(tahoe) */ - return zero/zero; -#endif /* defined(vax)||defined(tahoe) */ - } - -#if defined(vax)||defined(tahoe) - if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */ -#else /* defined(vax)||defined(tahoe) */ - if( ( *pp & mexp ) == mexp ) -#endif /* defined(vax)||defined(tahoe) */ - { if (p != p) return p; else return x;} - - else if ( ((*pp & mexp)>>gap) <= 1 ) - /* subnormal p, or almost subnormal p */ - { double b; b=scalb(1.0,(int)prep1); - p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} - else if ( p >= novf/2) - { p /= 2 ; x /= 2; return(drem(x,p)*2);} - else - { - dp=p+p; hp=p/2; - sign= *px & ~msign ; - *px &= msign ; - while ( x > dp ) - { - k=(*px & mexp) - (*pd & mexp) ; - tmp = dp ; - *pt += k ; - -#if defined(vax)||defined(tahoe) - if( x < tmp ) *pt -= 128 ; -#else /* defined(vax)||defined(tahoe) */ - if( x < tmp ) *pt -= 16 ; -#endif /* defined(vax)||defined(tahoe) */ - - x -= tmp ; - } - if ( x > hp ) - { x -= p ; if ( x >= hp ) x -= p ; } - -#if defined(vax)||defined(tahoe) - if (x) -#endif /* defined(vax)||defined(tahoe) */ - *px ^= sign; - return( x); - - } -} - - -double sqrt(x) -double x; -{ - double q,s,b,r; - double t; - double const zero=0.0; - int m,n,i; -#if defined(vax)||defined(tahoe) - int k=54; -#else /* defined(vax)||defined(tahoe) */ - int k=51; -#endif /* defined(vax)||defined(tahoe) */ - - /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ - if(x!=x||x==zero) return(x); - - /* sqrt(negative) is invalid */ - if(x<zero) { -#if defined(vax)||defined(tahoe) - return (infnan(EDOM)); /* NaN */ -#else /* defined(vax)||defined(tahoe) */ - return(zero/zero); -#endif /* defined(vax)||defined(tahoe) */ - } - - /* sqrt(INF) is INF */ - if(!finite(x)) return(x); - - /* scale x to [1,4) */ - n=logb(x); - x=scalb(x,-n); - if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */ - m += n; - n = m/2; - if((n+n)!=m) {x *= 2; m -=1; n=m/2;} - - /* generate sqrt(x) bit by bit (accumulating in q) */ - q=1.0; s=4.0; x -= 1.0; r=1; - for(i=1;i<=k;i++) { - t=s+1; x *= 4; r /= 2; - if(t<=x) { - s=t+t+2, x -= t; q += r;} - else - s *= 2; - } - - /* generate the last bit and determine the final rounding */ - r/=2; x *= 4; - if(x==zero) goto end; 100+r; /* trigger inexact flag */ - if(s<x) { - q+=r; x -=s; s += 2; s *= 2; x *= 4; - t = (x-s)-5; - b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */ - b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */ - if(t>=0) q+=r; } /* else: Round-to-nearest */ - else { - s *= 2; x *= 4; - t = (x-s)-1; - b=1.0+3*r/4; if(b==1.0) goto end; - b=1.0+r/4; if(b>1.0) t=1; - if(t>=0) q+=r; } - -end: return(scalb(q,n)); -} - -#if 0 -/* DREM(X,Y) - * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * INTENDED FOR ASSEMBLY LANGUAGE - * CODED IN C BY K.C. NG, 3/23/85, 4/8/85. - * - * Warning: this code should not get compiled in unless ALL of - * the following machine-dependent routines are supplied. - * - * Required machine dependent functions (not on a VAX): - * swapINX(i): save inexact flag and reset it to "i" - * swapENI(e): save inexact enable and reset it to "e" - */ - -double drem(x,y) -double x,y; -{ - -#ifdef national /* order of words in floating point number */ - static const n0=3,n1=2,n2=1,n3=0; -#else /* VAX, SUN, ZILOG, TAHOE */ - static const n0=0,n1=1,n2=2,n3=3; -#endif - - static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390; - static const double zero=0.0; - double hy,y1,t,t1; - short k; - long n; - int i,e; - unsigned short xexp,yexp, *px =(unsigned short *) &x , - nx,nf, *py =(unsigned short *) &y , - sign, *pt =(unsigned short *) &t , - *pt1 =(unsigned short *) &t1 ; - - xexp = px[n0] & mexp ; /* exponent of x */ - yexp = py[n0] & mexp ; /* exponent of y */ - sign = px[n0] &0x8000; /* sign of x */ - -/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */ - if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */ - if( xexp == mexp ) return(zero/zero); /* x is INF */ - if(y==zero) return(y/y); - -/* save the inexact flag and inexact enable in i and e respectively - * and reset them to zero - */ - i=swapINX(0); e=swapENI(0); - -/* subnormal number */ - nx=0; - if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;} - -/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */ - if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} - - nf=nx; - py[n0] &= 0x7fff; - px[n0] &= 0x7fff; - -/* mask off the least significant 27 bits of y */ - t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t; - -/* LOOP: argument reduction on x whenever x > y */ -loop: - while ( x > y ) - { - t=y; - t1=y1; - xexp=px[n0]&mexp; /* exponent of x */ - k=xexp-yexp-m25; - if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ - {pt[n0]+=k;pt1[n0]+=k;} - n=x/t; x=(x-n*t1)-n*(t-t1); - } - /* end while (x > y) */ - - if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} - -/* final adjustment */ - - hy=y/2.0; - if(x>hy||((x==hy)&&n%2==1)) x-=y; - px[n0] ^= sign; - if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} - -/* restore inexact flag and inexact enable */ - swapINX(i); swapENI(e); - - return(x); -} -#endif - -#if 0 -/* SQRT - * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT - * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE - * CODED IN C BY K.C. NG, 3/22/85. - * - * Warning: this code should not get compiled in unless ALL of - * the following machine-dependent routines are supplied. - * - * Required machine dependent functions: - * swapINX(i) ...return the status of INEXACT flag and reset it to "i" - * swapRM(r) ...return the current Rounding Mode and reset it to "r" - * swapENI(e) ...return the status of inexact enable and reset it to "e" - * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer - * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer - */ - -static const unsigned long table[] = { -0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, -58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, -21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; - -double newsqrt(x) -double x; -{ - double y,z,t,addc(),subc() - double const b54=134217728.*134217728.; /* b54=2**54 */ - long mx,scalx; - long const mexp=0x7ff00000; - int i,j,r,e,swapINX(),swapRM(),swapENI(); - unsigned long *py=(unsigned long *) &y , - *pt=(unsigned long *) &t , - *px=(unsigned long *) &x ; -#ifdef national /* ordering of word in a floating point number */ - const int n0=1, n1=0; -#else - const int n0=0, n1=1; -#endif -/* Rounding Mode: RN ...round-to-nearest - * RZ ...round-towards 0 - * RP ...round-towards +INF - * RM ...round-towards -INF - */ - const int RN=0,RZ=1,RP=2,RM=3; - /* machine dependent: work on a Zilog Z8070 - * and a National 32081 & 16081 - */ - -/* exceptions */ - if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ - if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */ - if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */ - -/* save, reset, initialize */ - e=swapENI(0); /* ...save and reset the inexact enable */ - i=swapINX(0); /* ...save INEXACT flag */ - r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */ - scalx=0; - -/* subnormal number, scale up x to x*2**54 */ - if(mx==0) {x *= b54 ; scalx-=0x01b00000;} - -/* scale x to avoid intermediate over/underflow: - * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */ - if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;} - if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;} - -/* magic initial approximation to almost 8 sig. bits */ - py[n0]=(px[n0]>>1)+0x1ff80000; - py[n0]=py[n0]-table[(py[n0]>>15)&31]; - -/* Heron's rule once with correction to improve y to almost 18 sig. bits */ - t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; - -/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ - t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; - t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; - -/* twiddle last bit to force y correctly rounded */ - swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ - swapINX(0); /* ...clear INEXACT flag */ - swapENI(e); /* ...restore inexact enable status */ - t=x/y; /* ...chopped quotient, possibly inexact */ - j=swapINX(i); /* ...read and restore inexact flag */ - if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */ - b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */ - if(r==RN) t=addc(t); /* ...t=t+ulp */ - else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */ - y=y+t; /* ...chopped sum */ - py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */ -end: py[n0]=py[n0]+scalx; /* ...scale back y */ - swapRM(r); /* ...restore Rounding Mode */ - return(y); -} -#endif |