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