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