(__expm1l, expm1l): Remove NO_LONG_DOUBLE aliases.
[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  * 4. Neither the name of the University nor the names of its contributors
14  *    may be used to endorse or promote products derived from this software
15  *    without specific prior written permission.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
18  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
21  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
23  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
25  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
26  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
27  * SUCH DAMAGE.
28  */
29
30 #ifndef lint
31 static char sccsid[] = "@(#)support.c   8.1 (Berkeley) 6/4/93";
32 #endif /* not lint */
33
34 /*
35  * Some IEEE standard 754 recommended functions and remainder and sqrt for
36  * supporting the C elementary functions.
37  ******************************************************************************
38  * WARNING:
39  *      These codes are developed (in double) to support the C elementary
40  * functions temporarily. They are not universal, and some of them are very
41  * slow (in particular, drem and sqrt is extremely inefficient). Each
42  * computer system should have its implementation of these functions using
43  * its own assembler.
44  ******************************************************************************
45  *
46  * IEEE 754 required operations:
47  *     drem(x,p)
48  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
49  *              nearest x/y; in half way case, choose the even one.
50  *     sqrt(x)
51  *              returns the square root of x correctly rounded according to
52  *              the rounding mod.
53  *
54  * IEEE 754 recommended functions:
55  * (a) copysign(x,y)
56  *              returns x with the sign of y.
57  * (b) scalb(x,N)
58  *              returns  x * (2**N), for integer values N.
59  * (c) logb(x)
60  *              returns the unbiased exponent of x, a signed integer in
61  *              double precision, except that logb(0) is -INF, logb(INF)
62  *              is +INF, and logb(NAN) is that NAN.
63  * (d) finite(x)
64  *              returns the value TRUE if -INF < x < +INF and returns
65  *              FALSE otherwise.
66  *
67  *
68  * CODED IN C BY K.C. NG, 11/25/84;
69  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
70  */
71
72 #include "mathimpl.h"
73
74 #if defined(vax)||defined(tahoe)      /* VAX D format */
75 #include <errno.h>
76     static const unsigned short msign=0x7fff , mexp =0x7f80 ;
77     static const short  prep1=57, gap=7, bias=129           ;
78     static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
79 #else   /* defined(vax)||defined(tahoe) */
80     static const unsigned short msign=0x7fff, mexp =0x7ff0  ;
81     static const short prep1=54, gap=4, bias=1023           ;
82     static const double novf=1.7E308, nunf=3.0E-308,zero=0.0;
83 #endif  /* defined(vax)||defined(tahoe) */
84
85 double scalb(x,N)
86 double x; int N;
87 {
88         int k;
89
90 #ifdef national
91         unsigned short *px=(unsigned short *) &x + 3;
92 #else   /* national */
93         unsigned short *px=(unsigned short *) &x;
94 #endif  /* national */
95
96         if( x == zero )  return(x);
97
98 #if defined(vax)||defined(tahoe)
99         if( (k= *px & mexp ) != ~msign ) {
100             if (N < -260)
101                 return(nunf*nunf);
102             else if (N > 260) {
103                 return(copysign(infnan(ERANGE),x));
104             }
105 #else   /* defined(vax)||defined(tahoe) */
106         if( (k= *px & mexp ) != mexp ) {
107             if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
108             if( k == 0 ) {
109                  x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
110 #endif  /* defined(vax)||defined(tahoe) */
111
112             if((k = (k>>gap)+ N) > 0 )
113                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
114                 else x=novf+novf;               /* overflow */
115             else
116                 if( k > -prep1 )
117                                         /* gradual underflow */
118                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
119                 else
120                 return(nunf*nunf);
121             }
122         return(x);
123 }
124
125
126 double copysign(x,y)
127 double x,y;
128 {
129 #ifdef national
130         unsigned short  *px=(unsigned short *) &x+3,
131                         *py=(unsigned short *) &y+3;
132 #else   /* national */
133         unsigned short  *px=(unsigned short *) &x,
134                         *py=(unsigned short *) &y;
135 #endif  /* national */
136
137 #if defined(vax)||defined(tahoe)
138         if ( (*px & mexp) == 0 ) return(x);
139 #endif  /* defined(vax)||defined(tahoe) */
140
141         *px = ( *px & msign ) | ( *py & ~msign );
142         return(x);
143 }
144
145 double logb(x)
146 double x;
147 {
148
149 #ifdef national
150         short *px=(short *) &x+3, k;
151 #else   /* national */
152         short *px=(short *) &x, k;
153 #endif  /* national */
154
155 #if defined(vax)||defined(tahoe)
156         return (int)(((*px&mexp)>>gap)-bias);
157 #else   /* defined(vax)||defined(tahoe) */
158         if( (k= *px & mexp ) != mexp )
159             if ( k != 0 )
160                 return ( (k>>gap) - bias );
161             else if( x != zero)
162                 return ( -1022.0 );
163             else
164                 return(-(1.0/zero));
165         else if(x != x)
166             return(x);
167         else
168             {*px &= msign; return(x);}
169 #endif  /* defined(vax)||defined(tahoe) */
170 }
171
172 finite(x)
173 double x;
174 {
175 #if defined(vax)||defined(tahoe)
176         return(1);
177 #else   /* defined(vax)||defined(tahoe) */
178 #ifdef national
179         return( (*((short *) &x+3 ) & mexp ) != mexp );
180 #else   /* national */
181         return( (*((short *) &x ) & mexp ) != mexp );
182 #endif  /* national */
183 #endif  /* defined(vax)||defined(tahoe) */
184 }
185
186 double drem(x,p)
187 double x,p;
188 {
189         short sign;
190         double hp,dp,tmp;
191         unsigned short  k;
192 #ifdef national
193         unsigned short
194               *px=(unsigned short *) &x  +3,
195               *pp=(unsigned short *) &p  +3,
196               *pd=(unsigned short *) &dp +3,
197               *pt=(unsigned short *) &tmp+3;
198 #else   /* national */
199         unsigned short
200               *px=(unsigned short *) &x  ,
201               *pp=(unsigned short *) &p  ,
202               *pd=(unsigned short *) &dp ,
203               *pt=(unsigned short *) &tmp;
204 #endif  /* national */
205
206         *pp &= msign ;
207
208 #if defined(vax)||defined(tahoe)
209         if( ( *px & mexp ) == ~msign )  /* is x a reserved operand? */
210 #else   /* defined(vax)||defined(tahoe) */
211         if( ( *px & mexp ) == mexp )
212 #endif  /* defined(vax)||defined(tahoe) */
213                 return  (x-p)-(x-p);    /* create nan if x is inf */
214         if (p == zero) {
215 #if defined(vax)||defined(tahoe)
216                 return(infnan(EDOM));
217 #else   /* defined(vax)||defined(tahoe) */
218                 return zero/zero;
219 #endif  /* defined(vax)||defined(tahoe) */
220         }
221
222 #if defined(vax)||defined(tahoe)
223         if( ( *pp & mexp ) == ~msign )  /* is p a reserved operand? */
224 #else   /* defined(vax)||defined(tahoe) */
225         if( ( *pp & mexp ) == mexp )
226 #endif  /* defined(vax)||defined(tahoe) */
227                 { if (p != p) return p; else return x;}
228
229         else  if ( ((*pp & mexp)>>gap) <= 1 )
230                 /* subnormal p, or almost subnormal p */
231             { double b; b=scalb(1.0,(int)prep1);
232               p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
233         else  if ( p >= novf/2)
234             { p /= 2 ; x /= 2; return(drem(x,p)*2);}
235         else
236             {
237                 dp=p+p; hp=p/2;
238                 sign= *px & ~msign ;
239                 *px &= msign       ;
240                 while ( x > dp )
241                     {
242                         k=(*px & mexp) - (*pd & mexp) ;
243                         tmp = dp ;
244                         *pt += k ;
245
246 #if defined(vax)||defined(tahoe)
247                         if( x < tmp ) *pt -= 128 ;
248 #else   /* defined(vax)||defined(tahoe) */
249                         if( x < tmp ) *pt -= 16 ;
250 #endif  /* defined(vax)||defined(tahoe) */
251
252                         x -= tmp ;
253                     }
254                 if ( x > hp )
255                     { x -= p ;  if ( x >= hp ) x -= p ; }
256
257 #if defined(vax)||defined(tahoe)
258                 if (x)
259 #endif  /* defined(vax)||defined(tahoe) */
260                         *px ^= sign;
261                 return( x);
262
263             }
264 }
265
266
267 double sqrt(x)
268 double x;
269 {
270         double q,s,b,r;
271         double t;
272         double const zero=0.0;
273         int m,n,i;
274 #if defined(vax)||defined(tahoe)
275         int k=54;
276 #else   /* defined(vax)||defined(tahoe) */
277         int k=51;
278 #endif  /* defined(vax)||defined(tahoe) */
279
280     /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
281         if(x!=x||x==zero) return(x);
282
283     /* sqrt(negative) is invalid */
284         if(x<zero) {
285 #if defined(vax)||defined(tahoe)
286                 return (infnan(EDOM));  /* NaN */
287 #else   /* defined(vax)||defined(tahoe) */
288                 return(zero/zero);
289 #endif  /* defined(vax)||defined(tahoe) */
290         }
291
292     /* sqrt(INF) is INF */
293         if(!finite(x)) return(x);
294
295     /* scale x to [1,4) */
296         n=logb(x);
297         x=scalb(x,-n);
298         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
299         m += n;
300         n = m/2;
301         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
302
303     /* generate sqrt(x) bit by bit (accumulating in q) */
304             q=1.0; s=4.0; x -= 1.0; r=1;
305             for(i=1;i<=k;i++) {
306                 t=s+1; x *= 4; r /= 2;
307                 if(t<=x) {
308                     s=t+t+2, x -= t; q += r;}
309                 else
310                     s *= 2;
311                 }
312
313     /* generate the last bit and determine the final rounding */
314             r/=2; x *= 4;
315             if(x==zero) goto end; 100+r; /* trigger inexact flag */
316             if(s<x) {
317                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
318                 t = (x-s)-5;
319                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
320                 b=1.0+r/4;   if(b>1.0) t=1;     /* b>1 : Round-to-(+INF) */
321                 if(t>=0) q+=r; }              /* else: Round-to-nearest */
322             else {
323                 s *= 2; x *= 4;
324                 t = (x-s)-1;
325                 b=1.0+3*r/4; if(b==1.0) goto end;
326                 b=1.0+r/4;   if(b>1.0) t=1;
327                 if(t>=0) q+=r; }
328
329 end:        return(scalb(q,n));
330 }
331
332 #if 0
333 /* DREM(X,Y)
334  * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
335  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
336  * INTENDED FOR ASSEMBLY LANGUAGE
337  * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
338  *
339  * Warning: this code should not get compiled in unless ALL of
340  * the following machine-dependent routines are supplied.
341  *
342  * Required machine dependent functions (not on a VAX):
343  *     swapINX(i): save inexact flag and reset it to "i"
344  *     swapENI(e): save inexact enable and reset it to "e"
345  */
346
347 double drem(x,y)
348 double x,y;
349 {
350
351 #ifdef national         /* order of words in floating point number */
352         static const n0=3,n1=2,n2=1,n3=0;
353 #else /* VAX, SUN, ZILOG, TAHOE */
354         static const n0=0,n1=1,n2=2,n3=3;
355 #endif
356
357         static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
358         static const double zero=0.0;
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 newsqrt(x)
449 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