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