Back to index

glibc  2.9
support.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 1985, 1993
00003  *     The Regents of the University of California.  All rights reserved.
00004  *
00005  * Redistribution and use in source and binary forms, with or without
00006  * modification, are permitted provided that the following conditions
00007  * are met:
00008  * 1. Redistributions of source code must retain the above copyright
00009  *    notice, this list of conditions and the following disclaimer.
00010  * 2. Redistributions in binary form must reproduce the above copyright
00011  *    notice, this list of conditions and the following disclaimer in the
00012  *    documentation and/or other materials provided with the distribution.
00013  * 4. Neither the name of the University nor the names of its contributors
00014  *    may be used to endorse or promote products derived from this software
00015  *    without specific prior written permission.
00016  *
00017  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
00018  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00019  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00020  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
00021  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
00022  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
00023  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00024  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00025  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
00026  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
00027  * SUCH DAMAGE.
00028  */
00029 
00030 #ifndef lint
00031 static char sccsid[] = "@(#)support.c     8.1 (Berkeley) 6/4/93";
00032 #endif /* not lint */
00033 
00034 /*
00035  * Some IEEE standard 754 recommended functions and remainder and sqrt for
00036  * supporting the C elementary functions.
00037  ******************************************************************************
00038  * WARNING:
00039  *      These codes are developed (in double) to support the C elementary
00040  * functions temporarily. They are not universal, and some of them are very
00041  * slow (in particular, drem and sqrt is extremely inefficient). Each
00042  * computer system should have its implementation of these functions using
00043  * its own assembler.
00044  ******************************************************************************
00045  *
00046  * IEEE 754 required operations:
00047  *     drem(x,p)
00048  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
00049  *              nearest x/y; in half way case, choose the even one.
00050  *     sqrt(x)
00051  *              returns the square root of x correctly rounded according to
00052  *            the rounding mod.
00053  *
00054  * IEEE 754 recommended functions:
00055  * (a) copysign(x,y)
00056  *              returns x with the sign of y.
00057  * (b) scalb(x,N)
00058  *              returns  x * (2**N), for integer values N.
00059  * (c) logb(x)
00060  *              returns the unbiased exponent of x, a signed integer in
00061  *              double precision, except that logb(0) is -INF, logb(INF)
00062  *              is +INF, and logb(NAN) is that NAN.
00063  * (d) finite(x)
00064  *              returns the value TRUE if -INF < x < +INF and returns
00065  *              FALSE otherwise.
00066  *
00067  *
00068  * CODED IN C BY K.C. NG, 11/25/84;
00069  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
00070  */
00071 
00072 #include "mathimpl.h"
00073 
00074 #if defined(vax)||defined(tahoe)      /* VAX D format */
00075 #include <errno.h>
00076     static const unsigned short msign=0x7fff , mexp =0x7f80 ;
00077     static const short  prep1=57, gap=7, bias=129           ;
00078     static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
00079 #else  /* defined(vax)||defined(tahoe) */
00080     static const unsigned short msign=0x7fff, mexp =0x7ff0  ;
00081     static const short prep1=54, gap=4, bias=1023           ;
00082     static const double novf=1.7E308, nunf=3.0E-308,zero=0.0;
00083 #endif /* defined(vax)||defined(tahoe) */
00084 
00085 double scalb(x,N)
00086 double x; int N;
00087 {
00088         int k;
00089 
00090 #ifdef national
00091         unsigned short *px=(unsigned short *) &x + 3;
00092 #else  /* national */
00093         unsigned short *px=(unsigned short *) &x;
00094 #endif /* national */
00095 
00096         if( x == zero )  return(x);
00097 
00098 #if defined(vax)||defined(tahoe)
00099         if( (k= *px & mexp ) != ~msign ) {
00100             if (N < -260)
00101               return(nunf*nunf);
00102            else if (N > 260) {
00103               return(copysign(infnan(ERANGE),x));
00104            }
00105 #else  /* defined(vax)||defined(tahoe) */
00106         if( (k= *px & mexp ) != mexp ) {
00107             if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
00108             if( k == 0 ) {
00109                  x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
00110 #endif /* defined(vax)||defined(tahoe) */
00111 
00112             if((k = (k>>gap)+ N) > 0 )
00113                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
00114                 else x=novf+novf;               /* overflow */
00115             else
00116                 if( k > -prep1 )
00117                                         /* gradual underflow */
00118                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
00119                 else
00120                 return(nunf*nunf);
00121             }
00122         return(x);
00123 }
00124 
00125 
00126 double copysign(x,y)
00127 double x,y;
00128 {
00129 #ifdef national
00130         unsigned short  *px=(unsigned short *) &x+3,
00131                         *py=(unsigned short *) &y+3;
00132 #else  /* national */
00133         unsigned short  *px=(unsigned short *) &x,
00134                         *py=(unsigned short *) &y;
00135 #endif /* national */
00136 
00137 #if defined(vax)||defined(tahoe)
00138         if ( (*px & mexp) == 0 ) return(x);
00139 #endif /* defined(vax)||defined(tahoe) */
00140 
00141         *px = ( *px & msign ) | ( *py & ~msign );
00142         return(x);
00143 }
00144 
00145 double logb(x)
00146 double x;
00147 {
00148 
00149 #ifdef national
00150         short *px=(short *) &x+3, k;
00151 #else  /* national */
00152         short *px=(short *) &x, k;
00153 #endif /* national */
00154 
00155 #if defined(vax)||defined(tahoe)
00156         return (int)(((*px&mexp)>>gap)-bias);
00157 #else  /* defined(vax)||defined(tahoe) */
00158         if( (k= *px & mexp ) != mexp )
00159             if ( k != 0 )
00160                 return ( (k>>gap) - bias );
00161             else if( x != zero)
00162                 return ( -1022.0 );
00163             else
00164                 return(-(1.0/zero));
00165         else if(x != x)
00166             return(x);
00167         else
00168             {*px &= msign; return(x);}
00169 #endif /* defined(vax)||defined(tahoe) */
00170 }
00171 
00172 finite(x)
00173 double x;
00174 {
00175 #if defined(vax)||defined(tahoe)
00176         return(1);
00177 #else  /* defined(vax)||defined(tahoe) */
00178 #ifdef national
00179         return( (*((short *) &x+3 ) & mexp ) != mexp );
00180 #else  /* national */
00181         return( (*((short *) &x ) & mexp ) != mexp );
00182 #endif /* national */
00183 #endif /* defined(vax)||defined(tahoe) */
00184 }
00185 
00186 double drem(x,p)
00187 double x,p;
00188 {
00189         short sign;
00190         double hp,dp,tmp;
00191         unsigned short  k;
00192 #ifdef national
00193         unsigned short
00194               *px=(unsigned short *) &x  +3,
00195               *pp=(unsigned short *) &p  +3,
00196               *pd=(unsigned short *) &dp +3,
00197               *pt=(unsigned short *) &tmp+3;
00198 #else  /* national */
00199         unsigned short
00200               *px=(unsigned short *) &x  ,
00201               *pp=(unsigned short *) &p  ,
00202               *pd=(unsigned short *) &dp ,
00203               *pt=(unsigned short *) &tmp;
00204 #endif /* national */
00205 
00206         *pp &= msign ;
00207 
00208 #if defined(vax)||defined(tahoe)
00209         if( ( *px & mexp ) == ~msign )    /* is x a reserved operand? */
00210 #else  /* defined(vax)||defined(tahoe) */
00211         if( ( *px & mexp ) == mexp )
00212 #endif /* defined(vax)||defined(tahoe) */
00213               return  (x-p)-(x-p); /* create nan if x is inf */
00214        if (p == zero) {
00215 #if defined(vax)||defined(tahoe)
00216               return(infnan(EDOM));
00217 #else  /* defined(vax)||defined(tahoe) */
00218               return zero/zero;
00219 #endif /* defined(vax)||defined(tahoe) */
00220        }
00221 
00222 #if defined(vax)||defined(tahoe)
00223         if( ( *pp & mexp ) == ~msign )    /* is p a reserved operand? */
00224 #else  /* defined(vax)||defined(tahoe) */
00225         if( ( *pp & mexp ) == mexp )
00226 #endif /* defined(vax)||defined(tahoe) */
00227               { if (p != p) return p; else return x;}
00228 
00229         else  if ( ((*pp & mexp)>>gap) <= 1 )
00230                 /* subnormal p, or almost subnormal p */
00231             { double b; b=scalb(1.0,(int)prep1);
00232               p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
00233         else  if ( p >= novf/2)
00234             { p /= 2 ; x /= 2; return(drem(x,p)*2);}
00235         else
00236             {
00237                 dp=p+p; hp=p/2;
00238                 sign= *px & ~msign ;
00239                 *px &= msign       ;
00240                 while ( x > dp )
00241                     {
00242                         k=(*px & mexp) - (*pd & mexp) ;
00243                         tmp = dp ;
00244                         *pt += k ;
00245 
00246 #if defined(vax)||defined(tahoe)
00247                         if( x < tmp ) *pt -= 128 ;
00248 #else  /* defined(vax)||defined(tahoe) */
00249                         if( x < tmp ) *pt -= 16 ;
00250 #endif /* defined(vax)||defined(tahoe) */
00251 
00252                         x -= tmp ;
00253                     }
00254                 if ( x > hp )
00255                     { x -= p ;  if ( x >= hp ) x -= p ; }
00256 
00257 #if defined(vax)||defined(tahoe)
00258               if (x)
00259 #endif /* defined(vax)||defined(tahoe) */
00260                      *px ^= sign;
00261                 return( x);
00262 
00263             }
00264 }
00265 
00266 
00267 double sqrt(x)
00268 double x;
00269 {
00270         double q,s,b,r;
00271         double t;
00272        double const zero=0.0;
00273         int m,n,i;
00274 #if defined(vax)||defined(tahoe)
00275         int k=54;
00276 #else  /* defined(vax)||defined(tahoe) */
00277         int k=51;
00278 #endif /* defined(vax)||defined(tahoe) */
00279 
00280     /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
00281         if(x!=x||x==zero) return(x);
00282 
00283     /* sqrt(negative) is invalid */
00284         if(x<zero) {
00285 #if defined(vax)||defined(tahoe)
00286               return (infnan(EDOM));      /* NaN */
00287 #else  /* defined(vax)||defined(tahoe) */
00288               return(zero/zero);
00289 #endif /* defined(vax)||defined(tahoe) */
00290        }
00291 
00292     /* sqrt(INF) is INF */
00293         if(!finite(x)) return(x);
00294 
00295     /* scale x to [1,4) */
00296         n=logb(x);
00297         x=scalb(x,-n);
00298         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
00299         m += n;
00300         n = m/2;
00301         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
00302 
00303     /* generate sqrt(x) bit by bit (accumulating in q) */
00304             q=1.0; s=4.0; x -= 1.0; r=1;
00305             for(i=1;i<=k;i++) {
00306                 t=s+1; x *= 4; r /= 2;
00307                 if(t<=x) {
00308                     s=t+t+2, x -= t; q += r;}
00309                 else
00310                     s *= 2;
00311                 }
00312 
00313     /* generate the last bit and determine the final rounding */
00314             r/=2; x *= 4;
00315             if(x==zero) goto end; 100+r; /* trigger inexact flag */
00316             if(s<x) {
00317                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
00318                 t = (x-s)-5;
00319                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
00320                 b=1.0+r/4;   if(b>1.0) t=1;      /* b>1 : Round-to-(+INF) */
00321                 if(t>=0) q+=r; }         /* else: Round-to-nearest */
00322             else {
00323                 s *= 2; x *= 4;
00324                 t = (x-s)-1;
00325                 b=1.0+3*r/4; if(b==1.0) goto end;
00326                 b=1.0+r/4;   if(b>1.0) t=1;
00327                 if(t>=0) q+=r; }
00328 
00329 end:        return(scalb(q,n));
00330 }
00331 
00332 #if 0
00333 /* DREM(X,Y)
00334  * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
00335  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
00336  * INTENDED FOR ASSEMBLY LANGUAGE
00337  * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
00338  *
00339  * Warning: this code should not get compiled in unless ALL of
00340  * the following machine-dependent routines are supplied.
00341  *
00342  * Required machine dependent functions (not on a VAX):
00343  *     swapINX(i): save inexact flag and reset it to "i"
00344  *     swapENI(e): save inexact enable and reset it to "e"
00345  */
00346 
00347 double drem(x,y)
00348 double x,y;
00349 {
00350 
00351 #ifdef national             /* order of words in floating point number */
00352        static const n0=3,n1=2,n2=1,n3=0;
00353 #else /* VAX, SUN, ZILOG, TAHOE */
00354        static const n0=0,n1=1,n2=2,n3=3;
00355 #endif
00356 
00357        static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
00358        static const double zero=0.0;
00359        double hy,y1,t,t1;
00360        short k;
00361        long n;
00362        int i,e;
00363        unsigned short xexp,yexp, *px  =(unsigned short *) &x  ,
00364                      nx,nf,   *py  =(unsigned short *) &y  ,
00365                      sign,    *pt  =(unsigned short *) &t  ,
00366                               *pt1 =(unsigned short *) &t1 ;
00367 
00368        xexp = px[n0] & mexp ;      /* exponent of x */
00369        yexp = py[n0] & mexp ;      /* exponent of y */
00370        sign = px[n0] &0x8000;      /* sign of x     */
00371 
00372 /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
00373        if(x!=x) return(x); if(y!=y) return(y);        /* x or y is NaN */
00374        if( xexp == mexp )   return(zero/zero);      /* x is INF */
00375        if(y==zero) return(y/y);
00376 
00377 /* save the inexact flag and inexact enable in i and e respectively
00378  * and reset them to zero
00379  */
00380        i=swapINX(0); e=swapENI(0);
00381 
00382 /* subnormal number */
00383        nx=0;
00384        if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
00385 
00386 /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
00387        if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
00388 
00389        nf=nx;
00390        py[n0] &= 0x7fff;
00391        px[n0] &= 0x7fff;
00392 
00393 /* mask off the least significant 27 bits of y */
00394        t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
00395 
00396 /* LOOP: argument reduction on x whenever x > y */
00397 loop:
00398        while ( x > y )
00399        {
00400            t=y;
00401            t1=y1;
00402            xexp=px[n0]&mexp;         /* exponent of x */
00403            k=xexp-yexp-m25;
00404            if(k>0)   /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
00405               {pt[n0]+=k;pt1[n0]+=k;}
00406            n=x/t; x=(x-n*t1)-n*(t-t1);
00407        }
00408     /* end while (x > y) */
00409 
00410        if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
00411 
00412 /* final adjustment */
00413 
00414        hy=y/2.0;
00415        if(x>hy||((x==hy)&&n%2==1)) x-=y;
00416        px[n0] ^= sign;
00417        if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
00418 
00419 /* restore inexact flag and inexact enable */
00420        swapINX(i); swapENI(e);
00421 
00422        return(x);
00423 }
00424 #endif
00425 
00426 #if 0
00427 /* SQRT
00428  * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
00429  * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
00430  * CODED IN C BY K.C. NG, 3/22/85.
00431  *
00432  * Warning: this code should not get compiled in unless ALL of
00433  * the following machine-dependent routines are supplied.
00434  *
00435  * Required machine dependent functions:
00436  *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
00437  *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
00438  *     swapENI(e)  ...return the status of inexact enable and reset it to "e"
00439  *     addc(t)     ...perform t=t+1 regarding t as a 64 bit unsigned integer
00440  *     subc(t)     ...perform t=t-1 regarding t as a 64 bit unsigned integer
00441  */
00442 
00443 static const unsigned long table[] = {
00444 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
00445 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
00446 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
00447 
00448 double newsqrt(x)
00449 double x;
00450 {
00451         double y,z,t,addc(),subc()
00452        double const b54=134217728.*134217728.; /* b54=2**54 */
00453         long mx,scalx;
00454        long const mexp=0x7ff00000;
00455         int i,j,r,e,swapINX(),swapRM(),swapENI();
00456         unsigned long *py=(unsigned long *) &y   ,
00457                       *pt=(unsigned long *) &t   ,
00458                       *px=(unsigned long *) &x   ;
00459 #ifdef national         /* ordering of word in a floating point number */
00460         const int n0=1, n1=0;
00461 #else
00462         const int n0=0, n1=1;
00463 #endif
00464 /* Rounding Mode:  RN ...round-to-nearest
00465  *                 RZ ...round-towards 0
00466  *                 RP ...round-towards +INF
00467  *               RM ...round-towards -INF
00468  */
00469         const int RN=0,RZ=1,RP=2,RM=3;
00470                             /* machine dependent: work on a Zilog Z8070
00471                                  * and a National 32081 & 16081
00472                                  */
00473 
00474 /* exceptions */
00475        if(x!=x||x==0.0) return(x);  /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
00476        if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
00477         if((mx=px[n0]&mexp)==mexp) return(x);  /* sqrt(+INF) is +INF */
00478 
00479 /* save, reset, initialize */
00480         e=swapENI(0);   /* ...save and reset the inexact enable */
00481         i=swapINX(0);   /* ...save INEXACT flag */
00482         r=swapRM(RN);   /* ...save and reset the Rounding Mode to RN */
00483         scalx=0;
00484 
00485 /* subnormal number, scale up x to x*2**54 */
00486         if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
00487 
00488 /* scale x to avoid intermediate over/underflow:
00489  * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
00490         if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
00491         if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
00492 
00493 /* magic initial approximation to almost 8 sig. bits */
00494         py[n0]=(px[n0]>>1)+0x1ff80000;
00495         py[n0]=py[n0]-table[(py[n0]>>15)&31];
00496 
00497 /* Heron's rule once with correction to improve y to almost 18 sig. bits */
00498         t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
00499 
00500 /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
00501         t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
00502         t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
00503 
00504 /* twiddle last bit to force y correctly rounded */
00505         swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
00506         swapINX(0);     /* ...clear INEXACT flag */
00507         swapENI(e);     /* ...restore inexact enable status */
00508         t=x/y;          /* ...chopped quotient, possibly inexact */
00509         j=swapINX(i);   /* ...read and restore inexact flag */
00510         if(j==0) { if(t==y) goto end; else t=subc(t); }  /* ...t=t-ulp */
00511         b54+0.1;        /* ..trigger inexact flag, sqrt(x) is inexact */
00512         if(r==RN) t=addc(t);            /* ...t=t+ulp */
00513         else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
00514         y=y+t;                          /* ...chopped sum */
00515         py[n0]=py[n0]-0x00100000;       /* ...correctly rounded sqrt(x) */
00516 end:    py[n0]=py[n0]+scalx;            /* ...scale back y */
00517         swapRM(r);                      /* ...restore Rounding Mode */
00518         return(y);
00519 }
00520 #endif