lostecho/oldlinux-files
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add directory gnu

Browse Source
This commit is contained in:
gohigh
2024-02-19 00:24:47 -05:00
parent 32616db5a4
commit a40f4cadb0
5086 changed files with 1860970 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

214
gnu/glibc/glibc-1.03/math/bsd/ieee/cabs.c Normal file
View File

@@ -0,0 +1,214 @@
/*
* Copyright (c) 1985 Regents of the University of California.
* All rights reserved.
*
* Redistribution and use in source and binary forms are permitted provided
* that: (1) source distributions retain this entire copyright notice and
* comment, and (2) distributions including binaries display the following
* acknowledgement: ``This product includes software developed by the
* University of California, Berkeley and its contributors'' in the
* documentation or other materials provided with the distribution and in
* all advertising materials mentioning features or use of this software.
* 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 ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
* WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*/
#ifndef lint
static char sccsid[] = "@(#)cabs.c 5.6 (Berkeley) 10/9/90";
#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) .
*/
double
cabs(z)
struct/* { double x, y;}*/ __complex z;
{
return hypot(z.__x,z.__y);
}
double
z_abs(z)
struct/* { double x,y;}*/ __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

106
gnu/glibc/glibc-1.03/math/bsd/ieee/cbrt.c Normal file
View File

@@ -0,0 +1,106 @@
/*
* Copyright (c) 1985 Regents of the University of California.
* All rights reserved.
*
* Redistribution and use in source and binary forms are permitted provided
* that: (1) source distributions retain this entire copyright notice and
* comment, and (2) distributions including binaries display the following
* acknowledgement: ``This product includes software developed by the
* University of California, Berkeley and its contributors'' in the
* documentation or other materials provided with the distribution and in
* all advertising materials mentioning features or use of this software.
* 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 ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
* WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*/
#ifndef lint
static char sccsid[] = "@(#)cbrt.c 5.8 (Berkeley) 10/9/90";
#endif /* not lint */
#include <sys/stdc.h>
/* 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

510
gnu/glibc/glibc-1.03/math/bsd/ieee/support.c Normal file
View File

@@ -0,0 +1,510 @@
/*
* Copyright (c) 1985 Regents of the University of California.
* All rights reserved.
*
* Redistribution and use in source and binary forms are permitted provided
* that: (1) source distributions retain this entire copyright notice and
* comment, and (2) distributions including binaries display the following
* acknowledgement: ``This product includes software developed by the
* University of California, Berkeley and its contributors'' in the
* documentation or other materials provided with the distribution and in
* all advertising materials mentioning features or use of this software.
* 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 ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
* WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*/
#ifndef lint
static char sccsid[] = "@(#)support.c 5.6 (Berkeley) 10/9/90";
#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