2002-08-26 Roland McGrath <roland@redhat.com>
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / support.c
index e976839..2651ea8 100644 (file)
  * 2. Redistributions in binary form must reproduce the above copyright
  *    notice, this list of conditions and the following disclaimer in the
  *    documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- *    must display the following acknowledgement:
- *     This product includes software developed by the University of
- *     California, Berkeley and its contributors.
  * 4. Neither the name of the University nor the names of its contributors
  *    may be used to endorse or promote products derived from this software
  *    without specific prior written permission.
 static char sccsid[] = "@(#)support.c  8.1 (Berkeley) 6/4/93";
 #endif /* not lint */
 
-/* 
- * Some IEEE standard 754 recommended functions and remainder and sqrt for 
+/*
+ * Some IEEE standard 754 recommended functions and remainder and sqrt for
  * supporting the C elementary functions.
  ******************************************************************************
  * WARNING:
  *      These codes are developed (in double) to support the C elementary
  * functions temporarily. They are not universal, and some of them are very
- * slow (in particular, drem and sqrt is extremely inefficient). Each 
- * computer system should have its implementation of these functions using 
+ * slow (in particular, drem and sqrt is extremely inefficient). Each
+ * computer system should have its implementation of these functions using
  * its own assembler.
  ******************************************************************************
  *
  * IEEE 754 required operations:
- *     drem(x,p) 
+ *     drem(x,p)
  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
  *              nearest x/y; in half way case, choose the even one.
- *     sqrt(x) 
- *              returns the square root of x correctly rounded according to 
+ *     sqrt(x)
+ *              returns the square root of x correctly rounded according to
  *             the rounding mod.
  *
  * IEEE 754 recommended functions:
- * (a) copysign(x,y) 
- *              returns x with the sign of y. 
- * (b) scalb(x,N) 
+ * (a) copysign(x,y)
+ *              returns x with the sign of y.
+ * (b) scalb(x,N)
  *              returns  x * (2**N), for integer values N.
- * (c) logb(x) 
- *              returns the unbiased exponent of x, a signed integer in 
- *              double precision, except that logb(0) is -INF, logb(INF) 
+ * (c) logb(x)
+ *              returns the unbiased exponent of x, a signed integer in
+ *              double precision, except that logb(0) is -INF, logb(INF)
  *              is +INF, and logb(NAN) is that NAN.
- * (d) finite(x) 
- *              returns the value TRUE if -INF < x < +INF and returns 
+ * (d) finite(x)
+ *              returns the value TRUE if -INF < x < +INF and returns
  *              FALSE otherwise.
  *
  *
@@ -78,7 +74,7 @@ static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93";
 #if defined(vax)||defined(tahoe)      /* VAX D format */
 #include <errno.h>
     static const unsigned short msign=0x7fff , mexp =0x7f80 ;
-    static const short  prep1=57, gap=7, bias=129           ;   
+    static const short  prep1=57, gap=7, bias=129           ;
     static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
 #else  /* defined(vax)||defined(tahoe) */
     static const unsigned short msign=0x7fff, mexp =0x7ff0  ;
@@ -97,7 +93,7 @@ double x; int N;
         unsigned short *px=(unsigned short *) &x;
 #endif /* national */
 
-        if( x == zero )  return(x); 
+        if( x == zero )  return(x);
 
 #if defined(vax)||defined(tahoe)
         if( (k= *px & mexp ) != ~msign ) {
@@ -117,7 +113,7 @@ double x; int N;
                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
                 else x=novf+novf;               /* overflow */
             else
-                if( k > -prep1 ) 
+                if( k > -prep1 )
                                         /* gradual underflow */
                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
                 else
@@ -147,7 +143,7 @@ double x,y;
 }
 
 double logb(x)
-double x; 
+double x;
 {
 
 #ifdef national
@@ -164,8 +160,8 @@ double x;
                 return ( (k>>gap) - bias );
             else if( x != zero)
                 return ( -1022.0 );
-            else        
-                return(-(1.0/zero));    
+            else
+                return(-(1.0/zero));
         else if(x != x)
             return(x);
         else
@@ -174,7 +170,7 @@ double x;
 }
 
 finite(x)
-double x;    
+double x;
 {
 #if defined(vax)||defined(tahoe)
         return(1);
@@ -192,16 +188,16 @@ double x,p;
 {
         short sign;
         double hp,dp,tmp;
-        unsigned short  k; 
+        unsigned short  k;
 #ifdef national
         unsigned short
-              *px=(unsigned short *) &x  +3, 
+              *px=(unsigned short *) &x  +3,
               *pp=(unsigned short *) &p  +3,
               *pd=(unsigned short *) &dp +3,
               *pt=(unsigned short *) &tmp+3;
 #else  /* national */
         unsigned short
-              *px=(unsigned short *) &x  , 
+              *px=(unsigned short *) &x  ,
               *pp=(unsigned short *) &p  ,
               *pd=(unsigned short *) &dp ,
               *pt=(unsigned short *) &tmp;
@@ -230,13 +226,13 @@ double x,p;
 #endif /* defined(vax)||defined(tahoe) */
                { if (p != p) return p; else return x;}
 
-        else  if ( ((*pp & mexp)>>gap) <= 1 ) 
+        else  if ( ((*pp & mexp)>>gap) <= 1 )
                 /* subnormal p, or almost subnormal p */
             { double b; b=scalb(1.0,(int)prep1);
               p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
         else  if ( p >= novf/2)
             { p /= 2 ; x /= 2; return(drem(x,p)*2);}
-        else 
+        else
             {
                 dp=p+p; hp=p/2;
                 sign= *px & ~msign ;
@@ -294,13 +290,13 @@ double x;
        }
 
     /* sqrt(INF) is INF */
-        if(!finite(x)) return(x);               
+        if(!finite(x)) return(x);
 
     /* scale x to [1,4) */
         n=logb(x);
         x=scalb(x,-n);
         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
-        m += n; 
+        m += n;
         n = m/2;
         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
 
@@ -313,23 +309,23 @@ double x;
                 else
                     s *= 2;
                 }
-            
+
     /* generate the last bit and determine the final rounding */
-            r/=2; x *= 4; 
+            r/=2; x *= 4;
             if(x==zero) goto end; 100+r; /* trigger inexact flag */
             if(s<x) {
                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
-                t = (x-s)-5; 
+                t = (x-s)-5;
                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
                 b=1.0+r/4;   if(b>1.0) t=1;    /* b>1 : Round-to-(+INF) */
                 if(t>=0) q+=r; }             /* else: Round-to-nearest */
-            else { 
-                s *= 2; x *= 4; 
-                t = (x-s)-1; 
+            else {
+                s *= 2; x *= 4;
+                t = (x-s)-1;
                 b=1.0+3*r/4; if(b==1.0) goto end;
                 b=1.0+r/4;   if(b>1.0) t=1;
                 if(t>=0) q+=r; }
-            
+
 end:        return(scalb(q,n));
 }
 
@@ -342,13 +338,13 @@ end:        return(scalb(q,n));
  *
  * Warning: this code should not get compiled in unless ALL of
  * the following machine-dependent routines are supplied.
- * 
+ *
  * Required machine dependent functions (not on a VAX):
  *     swapINX(i): save inexact flag and reset it to "i"
  *     swapENI(e): save inexact enable and reset it to "e"
  */
 
-double drem(x,y)       
+double drem(x,y)
 double x,y;
 {
 
@@ -363,8 +359,8 @@ double x,y;
        double hy,y1,t,t1;
        short k;
        long n;
-       int i,e; 
-       unsigned short xexp,yexp, *px  =(unsigned short *) &x  , 
+       int i,e;
+       unsigned short xexp,yexp, *px  =(unsigned short *) &x  ,
                        nx,nf,    *py  =(unsigned short *) &y  ,
                        sign,     *pt  =(unsigned short *) &t  ,
                                  *pt1 =(unsigned short *) &t1 ;
@@ -381,7 +377,7 @@ double x,y;
 /* save the inexact flag and inexact enable in i and e respectively
  * and reset them to zero
  */
-       i=swapINX(0);   e=swapENI(0);   
+       i=swapINX(0);   e=swapENI(0);
 
 /* subnormal number */
        nx=0;
@@ -391,7 +387,7 @@ double x,y;
        if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
 
        nf=nx;
-       py[n0] &= 0x7fff;       
+       py[n0] &= 0x7fff;
        px[n0] &= 0x7fff;
 
 /* mask off the least significant 27 bits of y */
@@ -408,7 +404,7 @@ loop:
            if(k>0)     /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
                {pt[n0]+=k;pt1[n0]+=k;}
            n=x/t; x=(x-n*t1)-n*(t-t1);
-       }       
+       }
     /* end while (x > y) */
 
        if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
@@ -416,14 +412,14 @@ loop:
 /* final adjustment */
 
        hy=y/2.0;
-       if(x>hy||((x==hy)&&n%2==1)) x-=y; 
+       if(x>hy||((x==hy)&&n%2==1)) x-=y;
        px[n0] ^= sign;
        if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
 
 /* restore inexact flag and inexact enable */
-       swapINX(i); swapENI(e); 
+       swapINX(i); swapENI(e);
 
-       return(x);      
+       return(x);
 }
 #endif
 
@@ -435,7 +431,7 @@ loop:
  *
  * Warning: this code should not get compiled in unless ALL of
  * the following machine-dependent routines are supplied.
- * 
+ *
  * Required machine dependent functions:
  *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
  *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
@@ -456,16 +452,16 @@ double x;
        double const b54=134217728.*134217728.; /* b54=2**54 */
         long mx,scalx;
        long const mexp=0x7ff00000;
-        int i,j,r,e,swapINX(),swapRM(),swapENI();       
+        int i,j,r,e,swapINX(),swapRM(),swapENI();
         unsigned long *py=(unsigned long *) &y   ,
                       *pt=(unsigned long *) &t   ,
                       *px=(unsigned long *) &x   ;
 #ifdef national         /* ordering of word in a floating point number */
-        const int n0=1, n1=0; 
+        const int n0=1, n1=0;
 #else
-        const int n0=0, n1=1; 
+        const int n0=0, n1=1;
 #endif
-/* Rounding Mode:  RN ...round-to-nearest 
+/* Rounding Mode:  RN ...round-to-nearest
  *                 RZ ...round-towards 0
  *                 RP ...round-towards +INF
  *                RM ...round-towards -INF
@@ -502,10 +498,10 @@ double x;
         t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
 
 /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
-        t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 
+        t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
         t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
 
-/* twiddle last bit to force y correctly rounded */ 
+/* twiddle last bit to force y correctly rounded */
         swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
         swapINX(0);     /* ...clear INEXACT flag */
         swapENI(e);     /* ...restore inexact enable status */