![[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.