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