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