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