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