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