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