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.7 (ucb.elefunt) 07/10/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 || p == zero ) 190 #else /* IEEE */ 191 if( ( *px & mexp ) == mexp || p == zero ) 192 #endif 193 return( (x != x)? x:zero/zero ); 194 195 else if ( ((*pp & mexp)>>gap) <= 1 ) 196 /* subnormal p, or almost subnormal p */ 197 { double b; b=scalb(1.0,(int)prep1); 198 p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} 199 else if ( p >= novf/2) 200 { p /= 2 ; x /= 2; return(drem(x,p)*2);} 201 else 202 { 203 dp=p+p; hp=p/2; 204 sign= *px & ~msign ; 205 *px &= msign ; 206 while ( x > dp ) 207 { 208 k=(*px & mexp) - (*pd & mexp) ; 209 tmp = dp ; 210 *pt += k ; 211 212 #if (defined(VAX)||defined(TAHOE)) 213 if( x < tmp ) *pt -= 128 ; 214 #else /* IEEE */ 215 if( x < tmp ) *pt -= 16 ; 216 #endif 217 218 x -= tmp ; 219 } 220 if ( x > hp ) 221 { x -= p ; if ( x >= hp ) x -= p ; } 222 223 #if (defined(VAX)||defined(TAHOE)) 224 if (x) 225 #endif 226 *px ^= sign; 227 return( x); 228 229 } 230 } 231 double sqrt(x) 232 double x; 233 { 234 double q,s,b,r; 235 double logb(),scalb(); 236 double t,zero=0.0; 237 int m,n,i,finite(); 238 #if (defined(VAX)||defined(TAHOE)) 239 int k=54; 240 #else /* IEEE */ 241 int k=51; 242 #endif 243 244 /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ 245 if(x!=x||x==zero) return(x); 246 247 /* sqrt(negative) is invalid */ 248 if(x<zero) return(zero/zero); 249 250 /* sqrt(INF) is INF */ 251 if(!finite(x)) return(x); 252 253 /* scale x to [1,4) */ 254 n=logb(x); 255 x=scalb(x,-n); 256 if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */ 257 m += n; 258 n = m/2; 259 if((n+n)!=m) {x *= 2; m -=1; n=m/2;} 260 261 /* generate sqrt(x) bit by bit (accumulating in q) */ 262 q=1.0; s=4.0; x -= 1.0; r=1; 263 for(i=1;i<=k;i++) { 264 t=s+1; x *= 4; r /= 2; 265 if(t<=x) { 266 s=t+t+2, x -= t; q += r;} 267 else 268 s *= 2; 269 } 270 271 /* generate the last bit and determine the final rounding */ 272 r/=2; x *= 4; 273 if(x==zero) goto end; 100+r; /* trigger inexact flag */ 274 if(s<x) { 275 q+=r; x -=s; s += 2; s *= 2; x *= 4; 276 t = (x-s)-5; 277 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */ 278 b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */ 279 if(t>=0) q+=r; } /* else: Round-to-nearest */ 280 else { 281 s *= 2; x *= 4; 282 t = (x-s)-1; 283 b=1.0+3*r/4; if(b==1.0) goto end; 284 b=1.0+r/4; if(b>1.0) t=1; 285 if(t>=0) q+=r; } 286 287 end: return(scalb(q,n)); 288 } 289 290 #if 0 291 /* DREM(X,Y) 292 * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) 293 * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 294 * INTENDED FOR ASSEMBLY LANGUAGE 295 * CODED IN C BY K.C. NG, 3/23/85, 4/8/85. 296 * 297 * Warning: this code should not get compiled in unless ALL of 298 * the following machine-dependent routines are supplied. 299 * 300 * Required machine dependent functions (not on a VAX): 301 * swapINX(i): save inexact flag and reset it to "i" 302 * swapENI(e): save inexact enable and reset it to "e" 303 */ 304 305 double drem(x,y) 306 double x,y; 307 { 308 309 #ifdef NATIONAL /* order of words in floating point number */ 310 static n0=3,n1=2,n2=1,n3=0; 311 #else /* VAX, SUN, ZILOG, TAHOE */ 312 static n0=0,n1=1,n2=2,n3=3; 313 #endif 314 315 static unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390; 316 static double zero=0.0; 317 double hy,y1,t,t1; 318 short k; 319 long n; 320 int i,e; 321 unsigned short xexp,yexp, *px =(unsigned short *) &x , 322 nx,nf, *py =(unsigned short *) &y , 323 sign, *pt =(unsigned short *) &t , 324 *pt1 =(unsigned short *) &t1 ; 325 326 xexp = px[n0] & mexp ; /* exponent of x */ 327 yexp = py[n0] & mexp ; /* exponent of y */ 328 sign = px[n0] &0x8000; /* sign of x */ 329 330 /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */ 331 if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */ 332 if( xexp == mexp ) return(zero/zero); /* x is INF */ 333 if(y==zero) return(y/y); 334 335 /* save the inexact flag and inexact enable in i and e respectively 336 * and reset them to zero 337 */ 338 i=swapINX(0); e=swapENI(0); 339 340 /* subnormal number */ 341 nx=0; 342 if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;} 343 344 /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */ 345 if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} 346 347 nf=nx; 348 py[n0] &= 0x7fff; 349 px[n0] &= 0x7fff; 350 351 /* mask off the least significant 27 bits of y */ 352 t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t; 353 354 /* LOOP: argument reduction on x whenever x > y */ 355 loop: 356 while ( x > y ) 357 { 358 t=y; 359 t1=y1; 360 xexp=px[n0]&mexp; /* exponent of x */ 361 k=xexp-yexp-m25; 362 if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ 363 {pt[n0]+=k;pt1[n0]+=k;} 364 n=x/t; x=(x-n*t1)-n*(t-t1); 365 } 366 /* end while (x > y) */ 367 368 if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} 369 370 /* final adjustment */ 371 372 hy=y/2.0; 373 if(x>hy||((x==hy)&&n%2==1)) x-=y; 374 px[n0] ^= sign; 375 if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} 376 377 /* restore inexact flag and inexact enable */ 378 swapINX(i); swapENI(e); 379 380 return(x); 381 } 382 #endif 383 384 #if 0 385 /* SQRT 386 * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT 387 * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE 388 * CODED IN C BY K.C. NG, 3/22/85. 389 * 390 * Warning: this code should not get compiled in unless ALL of 391 * the following machine-dependent routines are supplied. 392 * 393 * Required machine dependent functions: 394 * swapINX(i) ...return the status of INEXACT flag and reset it to "i" 395 * swapRM(r) ...return the current Rounding Mode and reset it to "r" 396 * swapENI(e) ...return the status of inexact enable and reset it to "e" 397 * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer 398 * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer 399 */ 400 401 static unsigned long table[] = { 402 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, 403 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, 404 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; 405 406 double newsqrt(x) 407 double x; 408 { 409 double y,z,t,addc(),subc(),b54=134217728.*134217728.; /* b54=2**54 */ 410 long mx,scalx,mexp=0x7ff00000; 411 int i,j,r,e,swapINX(),swapRM(),swapENI(); 412 unsigned long *py=(unsigned long *) &y , 413 *pt=(unsigned long *) &t , 414 *px=(unsigned long *) &x ; 415 #ifdef NATIONAL /* ordering of word in a floating point number */ 416 int n0=1, n1=0; 417 #else 418 int n0=0, n1=1; 419 #endif 420 /* Rounding Mode: RN ...round-to-nearest 421 * RZ ...round-towards 0 422 * RP ...round-towards +INF 423 * RM ...round-towards -INF 424 */ 425 int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070 426 * and a National 32081 & 16081 427 */ 428 429 /* exceptions */ 430 if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ 431 if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */ 432 if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */ 433 434 /* save, reset, initialize */ 435 e=swapENI(0); /* ...save and reset the inexact enable */ 436 i=swapINX(0); /* ...save INEXACT flag */ 437 r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */ 438 scalx=0; 439 440 /* subnormal number, scale up x to x*2**54 */ 441 if(mx==0) {x *= b54 ; scalx-=0x01b00000;} 442 443 /* scale x to avoid intermediate over/underflow: 444 * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */ 445 if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;} 446 if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;} 447 448 /* magic initial approximation to almost 8 sig. bits */ 449 py[n0]=(px[n0]>>1)+0x1ff80000; 450 py[n0]=py[n0]-table[(py[n0]>>15)&31]; 451 452 /* Heron's rule once with correction to improve y to almost 18 sig. bits */ 453 t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; 454 455 /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ 456 t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 457 t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; 458 459 /* twiddle last bit to force y correctly rounded */ 460 swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ 461 swapINX(0); /* ...clear INEXACT flag */ 462 swapENI(e); /* ...restore inexact enable status */ 463 t=x/y; /* ...chopped quotient, possibly inexact */ 464 j=swapINX(i); /* ...read and restore inexact flag */ 465 if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */ 466 b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */ 467 if(r==RN) t=addc(t); /* ...t=t+ulp */ 468 else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */ 469 y=y+t; /* ...chopped sum */ 470 py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */ 471 end: py[n0]=py[n0]+scalx; /* ...scale back y */ 472 swapRM(r); /* ...restore Rounding Mode */ 473 return(y); 474 } 475 #endif 476