1 /* 2 * Copyright (c) 1985 Regents of the University of California. 3 * 4 * Use and reproduction of this software are granted in accordance with 5 * the terms and conditions specified in the Berkeley Software License 6 * Agreement (in particular, this entails acknowledgement of the programs' 7 * source, and inclusion of this notice) with the additional understanding 8 * that all recipients should regard themselves as participants in an 9 * ongoing research project and hence should feel obligated to report 10 * their experiences (good or bad) with these elementary function codes, 11 * using "sendbug 4bsd-bugs@BERKELEY", to the authors. 12 */ 13 14 #ifndef lint 15 static char sccsid[] = 16 "@(#)support.c 1.1 (Berkeley) 5/23/85; 1.8 (ucb.elefunt) 07/11/87"; 17 #endif not lint 18 19 /* 20 * Some IEEE standard p754 recommended functions and remainder and sqrt for 21 * supporting the C elementary functions. 22 ****************************************************************************** 23 * WARNING: 24 * These codes are developed (in double) to support the C elementary 25 * functions temporarily. They are not universal, and some of them are very 26 * slow (in particular, drem and sqrt is extremely inefficient). Each 27 * computer system should have its implementation of these functions using 28 * its own assembler. 29 ****************************************************************************** 30 * 31 * IEEE p754 required operations: 32 * drem(x,p) 33 * returns x REM y = x - [x/y]*y , where [x/y] is the integer 34 * nearest x/y; in half way case, choose the even one. 35 * sqrt(x) 36 * returns the square root of x correctly rounded according to 37 * the rounding mod. 38 * 39 * IEEE p754 recommended functions: 40 * (a) copysign(x,y) 41 * returns x with the sign of y. 42 * (b) scalb(x,N) 43 * returns x * (2**N), for integer values N. 44 * (c) logb(x) 45 * returns the unbiased exponent of x, a signed integer in 46 * double precision, except that logb(0) is -INF, logb(INF) 47 * is +INF, and logb(NAN) is that NAN. 48 * (d) finite(x) 49 * returns the value TRUE if -INF < x < +INF and returns 50 * FALSE otherwise. 51 * 52 * 53 * CODED IN C BY K.C. NG, 11/25/84; 54 * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85. 55 */ 56 57 58 #if (defined(VAX)||defined(TAHOE)) /* VAX D format */ 59 static unsigned short msign=0x7fff , mexp =0x7f80 ; 60 static short prep1=57, gap=7, bias=129 ; 61 static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ; 62 #else /*IEEE double format */ 63 static unsigned short msign=0x7fff, mexp =0x7ff0 ; 64 static short prep1=54, gap=4, bias=1023 ; 65 static double novf=1.7E308, nunf=3.0E-308,zero=0.0; 66 #endif 67 68 double scalb(x,N) 69 double x; int N; 70 { 71 int k; 72 double scalb(); 73 74 #ifdef NATIONAL 75 unsigned short *px=(unsigned short *) &x + 3; 76 #else /* VAX, SUN, ZILOG, TAHOE */ 77 unsigned short *px=(unsigned short *) &x; 78 #endif 79 80 if( x == zero ) return(x); 81 82 #if (defined(VAX)||defined(TAHOE)) 83 if( (k= *px & mexp ) != ~msign ) { 84 if( N<-260) return(nunf*nunf); else if(N>260) return(novf+novf); 85 #else /* IEEE */ 86 if( (k= *px & mexp ) != mexp ) { 87 if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf); 88 if( k == 0 ) { 89 x *= scalb(1.0,(int)prep1); N -= prep1; return(scalb(x,N));} 90 #endif 91 92 if((k = (k>>gap)+ N) > 0 ) 93 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap); 94 else x=novf+novf; /* overflow */ 95 else 96 if( k > -prep1 ) 97 /* gradual underflow */ 98 {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);} 99 else 100 return(nunf*nunf); 101 } 102 return(x); 103 } 104 105 106 double copysign(x,y) 107 double x,y; 108 { 109 #ifdef NATIONAL 110 unsigned short *px=(unsigned short *) &x+3, 111 *py=(unsigned short *) &y+3; 112 #else /* VAX, SUN, ZILOG,TAHOE */ 113 unsigned short *px=(unsigned short *) &x, 114 *py=(unsigned short *) &y; 115 #endif 116 117 #if (defined(VAX)||defined(TAHOE)) 118 if ( (*px & mexp) == 0 ) return(x); 119 #endif 120 121 *px = ( *px & msign ) | ( *py & ~msign ); 122 return(x); 123 } 124 125 double logb(x) 126 double x; 127 { 128 129 #ifdef NATIONAL 130 short *px=(short *) &x+3, k; 131 #else /* VAX, SUN, ZILOG, TAHOE */ 132 short *px=(short *) &x, k; 133 #endif 134 135 #if (defined(VAX)||defined(TAHOE)) 136 return (int)(((*px&mexp)>>gap)-bias); 137 #else /* IEEE */ 138 if( (k= *px & mexp ) != mexp ) 139 if ( k != 0 ) 140 return ( (k>>gap) - bias ); 141 else if( x != zero) 142 return ( -1022.0 ); 143 else 144 return(-(1.0/zero)); 145 else if(x != x) 146 return(x); 147 else 148 {*px &= msign; return(x);} 149 #endif 150 } 151 152 finite(x) 153 double x; 154 { 155 #if (defined(VAX)||defined(TAHOE)) 156 return(1.0); 157 #else /* IEEE */ 158 #ifdef NATIONAL 159 return( (*((short *) &x+3 ) & mexp ) != mexp ); 160 #else /* SUN, ZILOG */ 161 return( (*((short *) &x ) & mexp ) != mexp ); 162 #endif 163 #endif 164 } 165 166 double drem(x,p) 167 double x,p; 168 { 169 short sign; 170 double hp,dp,tmp,drem(),scalb(); 171 unsigned short k; 172 #ifdef NATIONAL 173 unsigned short 174 *px=(unsigned short *) &x +3, 175 *pp=(unsigned short *) &p +3, 176 *pd=(unsigned short *) &dp +3, 177 *pt=(unsigned short *) &tmp+3; 178 #else /* VAX, SUN, ZILOG, TAHOE */ 179 unsigned short 180 *px=(unsigned short *) &x , 181 *pp=(unsigned short *) &p , 182 *pd=(unsigned short *) &dp , 183 *pt=(unsigned short *) &tmp; 184 #endif 185 186 *pp &= msign ; 187 188 #if (defined(VAX)||defined(TAHOE)) 189 if( ( *px & mexp ) == ~msign ) 190 #else /* IEEE */ 191 if( ( *px & mexp ) == mexp ) 192 #endif 193 return (x-p)-(x-p); /* create nan if x is inf */ 194 if(p==zero) return zero/zero; 195 #if (defined(VAX)||defined(TAHOE)) 196 if( ( *pp & mexp ) == ~msign ) 197 #else /* IEEE */ 198 if( ( *pp & mexp ) == mexp ) 199 #endif 200 { if (p != p) return p; else return x;} 201 202 else if ( ((*pp & mexp)>>gap) <= 1 ) 203 /* subnormal p, or almost subnormal p */ 204 { double b; b=scalb(1.0,(int)prep1); 205 p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} 206 else if ( p >= novf/2) 207 { p /= 2 ; x /= 2; return(drem(x,p)*2);} 208 else 209 { 210 dp=p+p; hp=p/2; 211 sign= *px & ~msign ; 212 *px &= msign ; 213 while ( x > dp ) 214 { 215 k=(*px & mexp) - (*pd & mexp) ; 216 tmp = dp ; 217 *pt += k ; 218 219 #if (defined(VAX)||defined(TAHOE)) 220 if( x < tmp ) *pt -= 128 ; 221 #else /* IEEE */ 222 if( x < tmp ) *pt -= 16 ; 223 #endif 224 225 x -= tmp ; 226 } 227 if ( x > hp ) 228 { x -= p ; if ( x >= hp ) x -= p ; } 229 230 #if (defined(VAX)||defined(TAHOE)) 231 if (x) 232 #endif 233 *px ^= sign; 234 return( x); 235 236 } 237 } 238 double sqrt(x) 239 double x; 240 { 241 double q,s,b,r; 242 double logb(),scalb(); 243 double t,zero=0.0; 244 int m,n,i,finite(); 245 #if (defined(VAX)||defined(TAHOE)) 246 int k=54; 247 #else /* IEEE */ 248 int k=51; 249 #endif 250 251 /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ 252 if(x!=x||x==zero) return(x); 253 254 /* sqrt(negative) is invalid */ 255 if(x<zero) return(zero/zero); 256 257 /* sqrt(INF) is INF */ 258 if(!finite(x)) return(x); 259 260 /* scale x to [1,4) */ 261 n=logb(x); 262 x=scalb(x,-n); 263 if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */ 264 m += n; 265 n = m/2; 266 if((n+n)!=m) {x *= 2; m -=1; n=m/2;} 267 268 /* generate sqrt(x) bit by bit (accumulating in q) */ 269 q=1.0; s=4.0; x -= 1.0; r=1; 270 for(i=1;i<=k;i++) { 271 t=s+1; x *= 4; r /= 2; 272 if(t<=x) { 273 s=t+t+2, x -= t; q += r;} 274 else 275 s *= 2; 276 } 277 278 /* generate the last bit and determine the final rounding */ 279 r/=2; x *= 4; 280 if(x==zero) goto end; 100+r; /* trigger inexact flag */ 281 if(s<x) { 282 q+=r; x -=s; s += 2; s *= 2; x *= 4; 283 t = (x-s)-5; 284 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */ 285 b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */ 286 if(t>=0) q+=r; } /* else: Round-to-nearest */ 287 else { 288 s *= 2; x *= 4; 289 t = (x-s)-1; 290 b=1.0+3*r/4; if(b==1.0) goto end; 291 b=1.0+r/4; if(b>1.0) t=1; 292 if(t>=0) q+=r; } 293 294 end: return(scalb(q,n)); 295 } 296 297 #if 0 298 /* DREM(X,Y) 299 * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) 300 * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 301 * INTENDED FOR ASSEMBLY LANGUAGE 302 * CODED IN C BY K.C. NG, 3/23/85, 4/8/85. 303 * 304 * Warning: this code should not get compiled in unless ALL of 305 * the following machine-dependent routines are supplied. 306 * 307 * Required machine dependent functions (not on a VAX): 308 * swapINX(i): save inexact flag and reset it to "i" 309 * swapENI(e): save inexact enable and reset it to "e" 310 */ 311 312 double drem(x,y) 313 double x,y; 314 { 315 316 #ifdef NATIONAL /* order of words in floating point number */ 317 static n0=3,n1=2,n2=1,n3=0; 318 #else /* VAX, SUN, ZILOG, TAHOE */ 319 static n0=0,n1=1,n2=2,n3=3; 320 #endif 321 322 static unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390; 323 static double zero=0.0; 324 double hy,y1,t,t1; 325 short k; 326 long n; 327 int i,e; 328 unsigned short xexp,yexp, *px =(unsigned short *) &x , 329 nx,nf, *py =(unsigned short *) &y , 330 sign, *pt =(unsigned short *) &t , 331 *pt1 =(unsigned short *) &t1 ; 332 333 xexp = px[n0] & mexp ; /* exponent of x */ 334 yexp = py[n0] & mexp ; /* exponent of y */ 335 sign = px[n0] &0x8000; /* sign of x */ 336 337 /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */ 338 if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */ 339 if( xexp == mexp ) return(zero/zero); /* x is INF */ 340 if(y==zero) return(y/y); 341 342 /* save the inexact flag and inexact enable in i and e respectively 343 * and reset them to zero 344 */ 345 i=swapINX(0); e=swapENI(0); 346 347 /* subnormal number */ 348 nx=0; 349 if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;} 350 351 /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */ 352 if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} 353 354 nf=nx; 355 py[n0] &= 0x7fff; 356 px[n0] &= 0x7fff; 357 358 /* mask off the least significant 27 bits of y */ 359 t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t; 360 361 /* LOOP: argument reduction on x whenever x > y */ 362 loop: 363 while ( x > y ) 364 { 365 t=y; 366 t1=y1; 367 xexp=px[n0]&mexp; /* exponent of x */ 368 k=xexp-yexp-m25; 369 if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ 370 {pt[n0]+=k;pt1[n0]+=k;} 371 n=x/t; x=(x-n*t1)-n*(t-t1); 372 } 373 /* end while (x > y) */ 374 375 if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} 376 377 /* final adjustment */ 378 379 hy=y/2.0; 380 if(x>hy||((x==hy)&&n%2==1)) x-=y; 381 px[n0] ^= sign; 382 if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} 383 384 /* restore inexact flag and inexact enable */ 385 swapINX(i); swapENI(e); 386 387 return(x); 388 } 389 #endif 390 391 #if 0 392 /* SQRT 393 * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT 394 * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE 395 * CODED IN C BY K.C. NG, 3/22/85. 396 * 397 * Warning: this code should not get compiled in unless ALL of 398 * the following machine-dependent routines are supplied. 399 * 400 * Required machine dependent functions: 401 * swapINX(i) ...return the status of INEXACT flag and reset it to "i" 402 * swapRM(r) ...return the current Rounding Mode and reset it to "r" 403 * swapENI(e) ...return the status of inexact enable and reset it to "e" 404 * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer 405 * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer 406 */ 407 408 static unsigned long table[] = { 409 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, 410 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, 411 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; 412 413 double newsqrt(x) 414 double x; 415 { 416 double y,z,t,addc(),subc(),b54=134217728.*134217728.; /* b54=2**54 */ 417 long mx,scalx,mexp=0x7ff00000; 418 int i,j,r,e,swapINX(),swapRM(),swapENI(); 419 unsigned long *py=(unsigned long *) &y , 420 *pt=(unsigned long *) &t , 421 *px=(unsigned long *) &x ; 422 #ifdef NATIONAL /* ordering of word in a floating point number */ 423 int n0=1, n1=0; 424 #else 425 int n0=0, n1=1; 426 #endif 427 /* Rounding Mode: RN ...round-to-nearest 428 * RZ ...round-towards 0 429 * RP ...round-towards +INF 430 * RM ...round-towards -INF 431 */ 432 int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070 433 * and a National 32081 & 16081 434 */ 435 436 /* exceptions */ 437 if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ 438 if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */ 439 if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */ 440 441 /* save, reset, initialize */ 442 e=swapENI(0); /* ...save and reset the inexact enable */ 443 i=swapINX(0); /* ...save INEXACT flag */ 444 r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */ 445 scalx=0; 446 447 /* subnormal number, scale up x to x*2**54 */ 448 if(mx==0) {x *= b54 ; scalx-=0x01b00000;} 449 450 /* scale x to avoid intermediate over/underflow: 451 * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */ 452 if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;} 453 if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;} 454 455 /* magic initial approximation to almost 8 sig. bits */ 456 py[n0]=(px[n0]>>1)+0x1ff80000; 457 py[n0]=py[n0]-table[(py[n0]>>15)&31]; 458 459 /* Heron's rule once with correction to improve y to almost 18 sig. bits */ 460 t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; 461 462 /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ 463 t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 464 t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; 465 466 /* twiddle last bit to force y correctly rounded */ 467 swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ 468 swapINX(0); /* ...clear INEXACT flag */ 469 swapENI(e); /* ...restore inexact enable status */ 470 t=x/y; /* ...chopped quotient, possibly inexact */ 471 j=swapINX(i); /* ...read and restore inexact flag */ 472 if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */ 473 b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */ 474 if(r==RN) t=addc(t); /* ...t=t+ulp */ 475 else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */ 476 y=y+t; /* ...chopped sum */ 477 py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */ 478 end: py[n0]=py[n0]+scalx; /* ...scale back y */ 479 swapRM(r); /* ...restore Rounding Mode */ 480 return(y); 481 } 482 #endif 483