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