1 /*
2 * Copyright (c) 1985, 1993
3 * The Regents of the University of California. All rights reserved.
4 *
5 * %sccs.include.redist.c%
6 */
7
8 #ifndef lint
9 static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 06/04/93";
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
scalb(x,N)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
copysign(x,y)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
logb(x)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
finite(x)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
drem(x,p)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
sqrt(x)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