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