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