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