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