xref: /csrg-svn/lib/libm/ieee/support.c (revision 34925)
134123Sbostic /*
224581Szliu  * Copyright (c) 1985 Regents of the University of California.
334123Sbostic  * All rights reserved.
434123Sbostic  *
534123Sbostic  * Redistribution and use in source and binary forms are permitted
6*34925Sbostic  * provided that the above copyright notice and this paragraph are
7*34925Sbostic  * duplicated in all such forms and that any documentation,
8*34925Sbostic  * advertising materials, and other materials related to such
9*34925Sbostic  * distribution and use acknowledge that the software was developed
10*34925Sbostic  * by the University of California, Berkeley.  The name of the
11*34925Sbostic  * University may not be used to endorse or promote products derived
12*34925Sbostic  * from this software without specific prior written permission.
13*34925Sbostic  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
14*34925Sbostic  * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
15*34925Sbostic  * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
1634123Sbostic  *
1734123Sbostic  * All recipients should regard themselves as participants in an ongoing
1834123Sbostic  * research project and hence should feel obligated to report their
1934123Sbostic  * experiences (good or bad) with these elementary function codes, using
2034123Sbostic  * the sendbug(8) program, to the authors.
2124581Szliu  */
2224581Szliu 
2324581Szliu #ifndef lint
24*34925Sbostic static char sccsid[] = "@(#)support.c	5.3 (Berkeley) 06/30/88";
2534123Sbostic #endif /* not lint */
2624581Szliu 
2724581Szliu /*
2831828Szliu  * Some IEEE standard 754 recommended functions and remainder and sqrt for
2924581Szliu  * supporting the C elementary functions.
3024581Szliu  ******************************************************************************
3124581Szliu  * WARNING:
3224581Szliu  *      These codes are developed (in double) to support the C elementary
3324581Szliu  * functions temporarily. They are not universal, and some of them are very
3424581Szliu  * slow (in particular, drem and sqrt is extremely inefficient). Each
3524581Szliu  * computer system should have its implementation of these functions using
3624581Szliu  * its own assembler.
3724581Szliu  ******************************************************************************
3824581Szliu  *
3931828Szliu  * IEEE 754 required operations:
4024581Szliu  *     drem(x,p)
4124581Szliu  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
4224581Szliu  *              nearest x/y; in half way case, choose the even one.
4324581Szliu  *     sqrt(x)
4424581Szliu  *              returns the square root of x correctly rounded according to
4524581Szliu  *		the rounding mod.
4624581Szliu  *
4731828Szliu  * IEEE 754 recommended functions:
4824581Szliu  * (a) copysign(x,y)
4924581Szliu  *              returns x with the sign of y.
5024581Szliu  * (b) scalb(x,N)
5124581Szliu  *              returns  x * (2**N), for integer values N.
5224581Szliu  * (c) logb(x)
5324581Szliu  *              returns the unbiased exponent of x, a signed integer in
5424581Szliu  *              double precision, except that logb(0) is -INF, logb(INF)
5524581Szliu  *              is +INF, and logb(NAN) is that NAN.
5624581Szliu  * (d) finite(x)
5724581Szliu  *              returns the value TRUE if -INF < x < +INF and returns
5824581Szliu  *              FALSE otherwise.
5924581Szliu  *
6024581Szliu  *
6124581Szliu  * CODED IN C BY K.C. NG, 11/25/84;
6224581Szliu  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
6324581Szliu  */
6424581Szliu 
6524581Szliu 
6631856Szliu #if defined(vax)||defined(tahoe)      /* VAX D format */
6731828Szliu #include <errno.h>
6824581Szliu     static unsigned short msign=0x7fff , mexp =0x7f80 ;
6924581Szliu     static short  prep1=57, gap=7, bias=129           ;
7024581Szliu     static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
7131856Szliu #else	/* defined(vax)||defined(tahoe) */
7224581Szliu     static unsigned short msign=0x7fff, mexp =0x7ff0  ;
7324581Szliu     static short prep1=54, gap=4, bias=1023           ;
7424581Szliu     static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
7531856Szliu #endif	/* defined(vax)||defined(tahoe) */
7624581Szliu 
7724581Szliu double scalb(x,N)
7824581Szliu double x; int N;
7924581Szliu {
8024581Szliu         int k;
8124581Szliu         double scalb();
8224581Szliu 
8331856Szliu #ifdef national
8424581Szliu         unsigned short *px=(unsigned short *) &x + 3;
8531856Szliu #else	/* national */
8624581Szliu         unsigned short *px=(unsigned short *) &x;
8731856Szliu #endif	/* national */
8824581Szliu 
8924581Szliu         if( x == zero )  return(x);
9024581Szliu 
9131856Szliu #if defined(vax)||defined(tahoe)
9224581Szliu         if( (k= *px & mexp ) != ~msign ) {
9331828Szliu             if (N < -260)
9431828Szliu 		return(nunf*nunf);
9531828Szliu 	    else if (N > 260) {
9631828Szliu 		extern double infnan(),copysign();
9731828Szliu 		return(copysign(infnan(ERANGE),x));
9831828Szliu 	    }
9931856Szliu #else	/* defined(vax)||defined(tahoe) */
10024581Szliu         if( (k= *px & mexp ) != mexp ) {
10124581Szliu             if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
10224581Szliu             if( k == 0 ) {
10324581Szliu                  x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
10431856Szliu #endif	/* defined(vax)||defined(tahoe) */
10524581Szliu 
10624581Szliu             if((k = (k>>gap)+ N) > 0 )
10724581Szliu                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
10824581Szliu                 else x=novf+novf;               /* overflow */
10924581Szliu             else
11024581Szliu                 if( k > -prep1 )
11124581Szliu                                         /* gradual underflow */
11224581Szliu                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
11324581Szliu                 else
11424581Szliu                 return(nunf*nunf);
11524581Szliu             }
11624581Szliu         return(x);
11724581Szliu }
11824581Szliu 
11924581Szliu 
12024581Szliu double copysign(x,y)
12124581Szliu double x,y;
12224581Szliu {
12331856Szliu #ifdef national
12424581Szliu         unsigned short  *px=(unsigned short *) &x+3,
12524581Szliu                         *py=(unsigned short *) &y+3;
12631856Szliu #else	/* national */
12724581Szliu         unsigned short  *px=(unsigned short *) &x,
12824581Szliu                         *py=(unsigned short *) &y;
12931856Szliu #endif	/* national */
13024581Szliu 
13131856Szliu #if defined(vax)||defined(tahoe)
13224581Szliu         if ( (*px & mexp) == 0 ) return(x);
13331856Szliu #endif	/* defined(vax)||defined(tahoe) */
13424581Szliu 
13524581Szliu         *px = ( *px & msign ) | ( *py & ~msign );
13624581Szliu         return(x);
13724581Szliu }
13824581Szliu 
13924581Szliu double logb(x)
14024581Szliu double x;
14124581Szliu {
14224581Szliu 
14331856Szliu #ifdef national
14424581Szliu         short *px=(short *) &x+3, k;
14531856Szliu #else	/* national */
14624581Szliu         short *px=(short *) &x, k;
14731856Szliu #endif	/* national */
14824581Szliu 
14931856Szliu #if defined(vax)||defined(tahoe)
15027449Szliu         return (int)(((*px&mexp)>>gap)-bias);
15131856Szliu #else	/* defined(vax)||defined(tahoe) */
15224581Szliu         if( (k= *px & mexp ) != mexp )
15324581Szliu             if ( k != 0 )
15424581Szliu                 return ( (k>>gap) - bias );
15524581Szliu             else if( x != zero)
15624581Szliu                 return ( -1022.0 );
15724581Szliu             else
15824581Szliu                 return(-(1.0/zero));
15924581Szliu         else if(x != x)
16024581Szliu             return(x);
16124581Szliu         else
16224581Szliu             {*px &= msign; return(x);}
16331856Szliu #endif	/* defined(vax)||defined(tahoe) */
16424581Szliu }
16524581Szliu 
16624581Szliu finite(x)
16724581Szliu double x;
16824581Szliu {
16931856Szliu #if defined(vax)||defined(tahoe)
17031828Szliu         return(1);
17131856Szliu #else	/* defined(vax)||defined(tahoe) */
17231856Szliu #ifdef national
17324581Szliu         return( (*((short *) &x+3 ) & mexp ) != mexp );
17431856Szliu #else	/* national */
17524581Szliu         return( (*((short *) &x ) & mexp ) != mexp );
17631856Szliu #endif	/* national */
17731856Szliu #endif	/* defined(vax)||defined(tahoe) */
17824581Szliu }
17924581Szliu 
18024581Szliu double drem(x,p)
18124581Szliu double x,p;
18224581Szliu {
18324581Szliu         short sign;
18424581Szliu         double hp,dp,tmp,drem(),scalb();
18524581Szliu         unsigned short  k;
18631856Szliu #ifdef national
18724581Szliu         unsigned short
18824581Szliu               *px=(unsigned short *) &x  +3,
18924581Szliu               *pp=(unsigned short *) &p  +3,
19024581Szliu               *pd=(unsigned short *) &dp +3,
19124581Szliu               *pt=(unsigned short *) &tmp+3;
19231856Szliu #else	/* national */
19324581Szliu         unsigned short
19424581Szliu               *px=(unsigned short *) &x  ,
19524581Szliu               *pp=(unsigned short *) &p  ,
19624581Szliu               *pd=(unsigned short *) &dp ,
19724581Szliu               *pt=(unsigned short *) &tmp;
19831856Szliu #endif	/* national */
19924581Szliu 
20024581Szliu         *pp &= msign ;
20124581Szliu 
20231856Szliu #if defined(vax)||defined(tahoe)
20331828Szliu         if( ( *px & mexp ) == ~msign )	/* is x a reserved operand? */
20431856Szliu #else	/* defined(vax)||defined(tahoe) */
20531827Szliu         if( ( *px & mexp ) == mexp )
20631856Szliu #endif	/* defined(vax)||defined(tahoe) */
20731827Szliu 		return  (x-p)-(x-p);	/* create nan if x is inf */
20831828Szliu 	if (p == zero) {
20931856Szliu #if defined(vax)||defined(tahoe)
21031828Szliu 		extern double infnan();
21131828Szliu 		return(infnan(EDOM));
21231856Szliu #else	/* defined(vax)||defined(tahoe) */
21331828Szliu 		return zero/zero;
21431856Szliu #endif	/* defined(vax)||defined(tahoe) */
21531828Szliu 	}
21631828Szliu 
21731856Szliu #if defined(vax)||defined(tahoe)
21831828Szliu         if( ( *pp & mexp ) == ~msign )	/* is p a reserved operand? */
21931856Szliu #else	/* defined(vax)||defined(tahoe) */
22031827Szliu         if( ( *pp & mexp ) == mexp )
22131856Szliu #endif	/* defined(vax)||defined(tahoe) */
22231827Szliu 		{ if (p != p) return p; else return x;}
22324581Szliu 
22424581Szliu         else  if ( ((*pp & mexp)>>gap) <= 1 )
22524581Szliu                 /* subnormal p, or almost subnormal p */
22624581Szliu             { double b; b=scalb(1.0,(int)prep1);
22724581Szliu               p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
22824581Szliu         else  if ( p >= novf/2)
22924581Szliu             { p /= 2 ; x /= 2; return(drem(x,p)*2);}
23024581Szliu         else
23124581Szliu             {
23224581Szliu                 dp=p+p; hp=p/2;
23324581Szliu                 sign= *px & ~msign ;
23424581Szliu                 *px &= msign       ;
23524581Szliu                 while ( x > dp )
23624581Szliu                     {
23724581Szliu                         k=(*px & mexp) - (*pd & mexp) ;
23824581Szliu                         tmp = dp ;
23924581Szliu                         *pt += k ;
24024581Szliu 
24131856Szliu #if defined(vax)||defined(tahoe)
24224581Szliu                         if( x < tmp ) *pt -= 128 ;
24331856Szliu #else	/* defined(vax)||defined(tahoe) */
24424581Szliu                         if( x < tmp ) *pt -= 16 ;
24531856Szliu #endif	/* defined(vax)||defined(tahoe) */
24624581Szliu 
24724581Szliu                         x -= tmp ;
24824581Szliu                     }
24924581Szliu                 if ( x > hp )
25024581Szliu                     { x -= p ;  if ( x >= hp ) x -= p ; }
25124581Szliu 
25231856Szliu #if defined(vax)||defined(tahoe)
25331825Szliu 		if (x)
25431856Szliu #endif	/* defined(vax)||defined(tahoe) */
25531825Szliu 			*px ^= sign;
25624581Szliu                 return( x);
25724581Szliu 
25824581Szliu             }
25924581Szliu }
26024581Szliu double sqrt(x)
26124581Szliu double x;
26224581Szliu {
26324581Szliu         double q,s,b,r;
26424581Szliu         double logb(),scalb();
26524581Szliu         double t,zero=0.0;
26624581Szliu         int m,n,i,finite();
26731856Szliu #if defined(vax)||defined(tahoe)
26824581Szliu         int k=54;
26931856Szliu #else	/* defined(vax)||defined(tahoe) */
27024581Szliu         int k=51;
27131856Szliu #endif	/* defined(vax)||defined(tahoe) */
27224581Szliu 
27324581Szliu     /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
27424581Szliu         if(x!=x||x==zero) return(x);
27524581Szliu 
27624581Szliu     /* sqrt(negative) is invalid */
27731828Szliu         if(x<zero) {
27831856Szliu #if defined(vax)||defined(tahoe)
27931828Szliu 		extern double infnan();
28031828Szliu 		return (infnan(EDOM));	/* NaN */
28131856Szliu #else	/* defined(vax)||defined(tahoe) */
28231828Szliu 		return(zero/zero);
28331856Szliu #endif	/* defined(vax)||defined(tahoe) */
28431828Szliu 	}
28524581Szliu 
28624581Szliu     /* sqrt(INF) is INF */
28724581Szliu         if(!finite(x)) return(x);
28824581Szliu 
28924581Szliu     /* scale x to [1,4) */
29024581Szliu         n=logb(x);
29124581Szliu         x=scalb(x,-n);
29224581Szliu         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
29324581Szliu         m += n;
29424581Szliu         n = m/2;
29524581Szliu         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
29624581Szliu 
29724581Szliu     /* generate sqrt(x) bit by bit (accumulating in q) */
29824581Szliu             q=1.0; s=4.0; x -= 1.0; r=1;
29924581Szliu             for(i=1;i<=k;i++) {
30024581Szliu                 t=s+1; x *= 4; r /= 2;
30124581Szliu                 if(t<=x) {
30224581Szliu                     s=t+t+2, x -= t; q += r;}
30324581Szliu                 else
30424581Szliu                     s *= 2;
30524581Szliu                 }
30624581Szliu 
30724581Szliu     /* generate the last bit and determine the final rounding */
30824581Szliu             r/=2; x *= 4;
30924581Szliu             if(x==zero) goto end; 100+r; /* trigger inexact flag */
31024581Szliu             if(s<x) {
31124581Szliu                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
31224581Szliu                 t = (x-s)-5;
31324581Szliu                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
31424581Szliu                 b=1.0+r/4;   if(b>1.0) t=1;	/* b>1 : Round-to-(+INF) */
31524581Szliu                 if(t>=0) q+=r; }	      /* else: Round-to-nearest */
31624581Szliu             else {
31724581Szliu                 s *= 2; x *= 4;
31824581Szliu                 t = (x-s)-1;
31924581Szliu                 b=1.0+3*r/4; if(b==1.0) goto end;
32024581Szliu                 b=1.0+r/4;   if(b>1.0) t=1;
32124581Szliu                 if(t>=0) q+=r; }
32224581Szliu 
32324581Szliu end:        return(scalb(q,n));
32424581Szliu }
32524581Szliu 
32624581Szliu #if 0
32724581Szliu /* DREM(X,Y)
32824581Szliu  * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
32924581Szliu  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
33024581Szliu  * INTENDED FOR ASSEMBLY LANGUAGE
33124581Szliu  * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
33224581Szliu  *
33324581Szliu  * Warning: this code should not get compiled in unless ALL of
33424581Szliu  * the following machine-dependent routines are supplied.
33524581Szliu  *
33624581Szliu  * Required machine dependent functions (not on a VAX):
33724581Szliu  *     swapINX(i): save inexact flag and reset it to "i"
33824581Szliu  *     swapENI(e): save inexact enable and reset it to "e"
33924581Szliu  */
34024581Szliu 
34124581Szliu double drem(x,y)
34224581Szliu double x,y;
34324581Szliu {
34424581Szliu 
34531856Szliu #ifdef national		/* order of words in floating point number */
34624581Szliu 	static n0=3,n1=2,n2=1,n3=0;
34731789Szliu #else /* VAX, SUN, ZILOG, TAHOE */
34824581Szliu 	static n0=0,n1=1,n2=2,n3=3;
34924581Szliu #endif
35024581Szliu 
35124581Szliu     	static unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
35224581Szliu 	static double zero=0.0;
35324581Szliu 	double hy,y1,t,t1;
35424581Szliu 	short k;
35524581Szliu 	long n;
35624581Szliu 	int i,e;
35724581Szliu 	unsigned short xexp,yexp, *px  =(unsigned short *) &x  ,
35824581Szliu 	      		nx,nf,	  *py  =(unsigned short *) &y  ,
35924581Szliu 	      		sign,	  *pt  =(unsigned short *) &t  ,
36024581Szliu 	      			  *pt1 =(unsigned short *) &t1 ;
36124581Szliu 
36224581Szliu 	xexp = px[n0] & mexp ;	/* exponent of x */
36324581Szliu 	yexp = py[n0] & mexp ;	/* exponent of y */
36424581Szliu 	sign = px[n0] &0x8000;	/* sign of x     */
36524581Szliu 
36624581Szliu /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
36724581Szliu 	if(x!=x) return(x); if(y!=y) return(y);	     /* x or y is NaN */
36824581Szliu 	if( xexp == mexp )   return(zero/zero);      /* x is INF */
36924581Szliu 	if(y==zero) return(y/y);
37024581Szliu 
37124581Szliu /* save the inexact flag and inexact enable in i and e respectively
37224581Szliu  * and reset them to zero
37324581Szliu  */
37424581Szliu 	i=swapINX(0);	e=swapENI(0);
37524581Szliu 
37624581Szliu /* subnormal number */
37724581Szliu 	nx=0;
37824581Szliu 	if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
37924581Szliu 
38024581Szliu /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
38124581Szliu 	if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
38224581Szliu 
38324581Szliu 	nf=nx;
38424581Szliu 	py[n0] &= 0x7fff;
38524581Szliu 	px[n0] &= 0x7fff;
38624581Szliu 
38724581Szliu /* mask off the least significant 27 bits of y */
38824581Szliu 	t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
38924581Szliu 
39024581Szliu /* LOOP: argument reduction on x whenever x > y */
39124581Szliu loop:
39224581Szliu 	while ( x > y )
39324581Szliu 	{
39424581Szliu 	    t=y;
39524581Szliu 	    t1=y1;
39624581Szliu 	    xexp=px[n0]&mexp;	  /* exponent of x */
39724581Szliu 	    k=xexp-yexp-m25;
39824581Szliu 	    if(k>0) 	/* if x/y >= 2**26, scale up y so that x/y < 2**26 */
39924581Szliu 		{pt[n0]+=k;pt1[n0]+=k;}
40024581Szliu 	    n=x/t; x=(x-n*t1)-n*(t-t1);
40124581Szliu 	}
40224581Szliu     /* end while (x > y) */
40324581Szliu 
40424581Szliu 	if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
40524581Szliu 
40624581Szliu /* final adjustment */
40724581Szliu 
40824581Szliu 	hy=y/2.0;
40924581Szliu 	if(x>hy||((x==hy)&&n%2==1)) x-=y;
41024581Szliu 	px[n0] ^= sign;
41124581Szliu 	if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
41224581Szliu 
41324581Szliu /* restore inexact flag and inexact enable */
41424581Szliu 	swapINX(i); swapENI(e);
41524581Szliu 
41624581Szliu 	return(x);
41724581Szliu }
41824581Szliu #endif
41924581Szliu 
42024581Szliu #if 0
42124581Szliu /* SQRT
42224581Szliu  * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
42324581Szliu  * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
42424581Szliu  * CODED IN C BY K.C. NG, 3/22/85.
42524581Szliu  *
42624581Szliu  * Warning: this code should not get compiled in unless ALL of
42724581Szliu  * the following machine-dependent routines are supplied.
42824581Szliu  *
42924581Szliu  * Required machine dependent functions:
43024581Szliu  *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
43124581Szliu  *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
43224581Szliu  *     swapENI(e)  ...return the status of inexact enable and reset it to "e"
43324581Szliu  *     addc(t)     ...perform t=t+1 regarding t as a 64 bit unsigned integer
43424581Szliu  *     subc(t)     ...perform t=t-1 regarding t as a 64 bit unsigned integer
43524581Szliu  */
43624581Szliu 
43724581Szliu static unsigned long table[] = {
43824581Szliu 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
43924581Szliu 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
44024581Szliu 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
44124581Szliu 
44224581Szliu double newsqrt(x)
44324581Szliu double x;
44424581Szliu {
44524581Szliu         double y,z,t,addc(),subc(),b54=134217728.*134217728.; /* b54=2**54 */
44624581Szliu         long mx,scalx,mexp=0x7ff00000;
44724581Szliu         int i,j,r,e,swapINX(),swapRM(),swapENI();
44824581Szliu         unsigned long *py=(unsigned long *) &y   ,
44924581Szliu                       *pt=(unsigned long *) &t   ,
45024581Szliu                       *px=(unsigned long *) &x   ;
45131856Szliu #ifdef national         /* ordering of word in a floating point number */
45224581Szliu         int n0=1, n1=0;
45324581Szliu #else
45424581Szliu         int n0=0, n1=1;
45524581Szliu #endif
45624581Szliu /* Rounding Mode:  RN ...round-to-nearest
45724581Szliu  *                 RZ ...round-towards 0
45824581Szliu  *                 RP ...round-towards +INF
45924581Szliu  *		   RM ...round-towards -INF
46024581Szliu  */
46124581Szliu         int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070
46224581Szliu                                  * and a National 32081 & 16081
46324581Szliu                                  */
46424581Szliu 
46524581Szliu /* exceptions */
46624581Szliu 	if(x!=x||x==0.0) return(x);  /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
46724581Szliu 	if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
46824581Szliu         if((mx=px[n0]&mexp)==mexp) return(x);  /* sqrt(+INF) is +INF */
46924581Szliu 
47024581Szliu /* save, reset, initialize */
47124581Szliu         e=swapENI(0);   /* ...save and reset the inexact enable */
47224581Szliu         i=swapINX(0);   /* ...save INEXACT flag */
47324581Szliu         r=swapRM(RN);   /* ...save and reset the Rounding Mode to RN */
47424581Szliu         scalx=0;
47524581Szliu 
47624581Szliu /* subnormal number, scale up x to x*2**54 */
47724581Szliu         if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
47824581Szliu 
47924581Szliu /* scale x to avoid intermediate over/underflow:
48024581Szliu  * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
48124581Szliu         if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
48224581Szliu         if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
48324581Szliu 
48424581Szliu /* magic initial approximation to almost 8 sig. bits */
48524581Szliu         py[n0]=(px[n0]>>1)+0x1ff80000;
48624581Szliu         py[n0]=py[n0]-table[(py[n0]>>15)&31];
48724581Szliu 
48824581Szliu /* Heron's rule once with correction to improve y to almost 18 sig. bits */
48924581Szliu         t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
49024581Szliu 
49124581Szliu /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
49224581Szliu         t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
49324581Szliu         t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
49424581Szliu 
49524581Szliu /* twiddle last bit to force y correctly rounded */
49624581Szliu         swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
49724581Szliu         swapINX(0);     /* ...clear INEXACT flag */
49824581Szliu         swapENI(e);     /* ...restore inexact enable status */
49924581Szliu         t=x/y;          /* ...chopped quotient, possibly inexact */
50024581Szliu         j=swapINX(i);   /* ...read and restore inexact flag */
50124581Szliu         if(j==0) { if(t==y) goto end; else t=subc(t); }  /* ...t=t-ulp */
50224581Szliu         b54+0.1;        /* ..trigger inexact flag, sqrt(x) is inexact */
50324581Szliu         if(r==RN) t=addc(t);            /* ...t=t+ulp */
50424581Szliu         else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
50524581Szliu         y=y+t;                          /* ...chopped sum */
50624581Szliu         py[n0]=py[n0]-0x00100000;       /* ...correctly rounded sqrt(x) */
50724581Szliu end:    py[n0]=py[n0]+scalx;            /* ...scale back y */
50824581Szliu         swapRM(r);                      /* ...restore Rounding Mode */
50924581Szliu         return(y);
51024581Szliu }
51124581Szliu #endif
512