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