lib/libm/noieee_src/n_support.c

   1 /*      $NetBSD: n_support.c,v 1.5 2003/08/07 16:44:52 agc Exp $ */
   2 /*
   3  * Copyright (c) 1985, 1993
   4  *      The Regents of the University of California.  All rights reserved.
   5  *
   6  * Redistribution and use in source and binary forms, with or without
   7  * modification, are permitted provided that the following conditions
   8  * are met:
   9  * 1. Redistributions of source code must retain the above copyright
  10  *    notice, this list of conditions and the following disclaimer.
  11  * 2. Redistributions in binary form must reproduce the above copyright
  12  *    notice, this list of conditions and the following disclaimer in the
  13  *    documentation and/or other materials provided with the distribution.
  14  * 3. Neither the name of the University nor the names of its contributors
  15  *    may be used to endorse or promote products derived from this software
  16  *    without specific prior written permission.
  17  *
  18  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  19  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  22  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  23  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  24  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  25  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  26  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  27  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  28  * SUCH DAMAGE.
  29  */
  30
  31 #ifndef lint
  32 static char sccsid[] = "@(#)support.c   8.1 (Berkeley) 6/4/93";
  33 #endif /* not lint */
  34
  35 /*
  36  * Some IEEE standard 754 recommended functions and remainder and sqrt for
  37  * supporting the C elementary functions.
  38  ******************************************************************************
  39  * WARNING:
  40  *      These codes are developed (in double) to support the C elementary
  41  * functions temporarily. They are not universal, and some of them are very
  42  * slow (in particular, drem and sqrt is extremely inefficient). Each
  43  * computer system should have its implementation of these functions using
  44  * its own assembler.
  45  ******************************************************************************
  46  *
  47  * IEEE 754 required operations:
  48  *     drem(x,p)
  49  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
  50  *              nearest x/y; in half way case, choose the even one.
  51  *     sqrt(x)
  52  *              returns the square root of x correctly rounded according to
  53  *              the rounding mod.
  54  *
  55  * IEEE 754 recommended functions:
  56  * (a) copysign(x,y)
  57  *              returns x with the sign of y.
  58  * (b) scalb(x,N)
  59  *              returns  x * (2**N), for integer values N.
  60  * (c) logb(x)
  61  *              returns the unbiased exponent of x, a signed integer in
  62  *              double precision, except that logb(0) is -INF, logb(INF)
  63  *              is +INF, and logb(NAN) is that NAN.
  64  * (d) finite(x)
  65  *              returns the value TRUE if -INF < x < +INF and returns
  66  *              FALSE otherwise.
  67  *
  68  *
  69  * CODED IN C BY K.C. NG, 11/25/84;
  70  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
  71  */
  72
  73 #include "mathimpl.h"
  74 #include "trig.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 ;
  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;
  85 #endif  /* defined(__vax__)||defined(tahoe) */
  86
  87 double
  88 scalb(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
 129 copysign(double x, double 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
 148 logb(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 int
 175 finite(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
 189 drem(double x, double 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
 270 sqrt(double x)
 271 {
 272         double q,s,b,r;
 273         double t;
 274         int m,n,i;
 275 #if defined(__vax__)||defined(tahoe)
 276         int k=54;
 277 #else   /* defined(__vax__)||defined(tahoe) */
 278         int k=51;
 279 #endif  /* defined(__vax__)||defined(tahoe) */
 280
 281     /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
 282         if(x!=x||x==__zero) return(x);
 283
 284     /* sqrt(negative) is invalid */
 285         if(x<__zero) {
 286 #if defined(__vax__)||defined(tahoe)
 287                 return (infnan(EDOM));  /* NaN */
 288 #else   /* defined(__vax__)||defined(tahoe) */
 289                 return(__zero/__zero);
 290 #endif  /* defined(__vax__)||defined(tahoe) */
 291         }
 292
 293     /* sqrt(INF) is INF */
 294         if(!finite(x)) return(x);
 295
 296     /* scale x to [1,4) */
 297         n=logb(x);
 298         x=scalb(x,-n);
 299         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
 300         m += n;
 301         n = m/2;
 302         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
 303
 304     /* generate sqrt(x) bit by bit (accumulating in q) */
 305             q=1.0; s=4.0; x -= 1.0; r=1;
 306             for(i=1;i<=k;i++) {
 307                 t=s+1; x *= 4; r /= 2;
 308                 if(t<=x) {
 309                     s=t+t+2, x -= t; q += r;}
 310                 else
 311                     s *= 2;
 312                 }
 313
 314     /* generate the last bit and determine the final rounding */
 315             r/=2; x *= 4;
 316             if(x==__zero) goto end; 100+r; /* trigger inexact flag */
 317             if(s<x) {
 318                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
 319                 t = (x-s)-5;
 320                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
 321                 b=1.0+r/4;   if(b>1.0) t=1;     /* b>1 : Round-to-(+INF) */
 322                 if(t>=0) q+=r; }              /* else: Round-to-nearest */
 323             else {
 324                 s *= 2; x *= 4;
 325                 t = (x-s)-1;
 326                 b=1.0+3*r/4; if(b==1.0) goto end;
 327                 b=1.0+r/4;   if(b>1.0) t=1;
 328                 if(t>=0) q+=r; }
 329
 330 end:        return(scalb(q,n));
 331 }
 332
 333 #if 0
 334 /* DREM(X,Y)
 335  * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
 336  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
 337  * INTENDED FOR ASSEMBLY LANGUAGE
 338  * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
 339  *
 340  * Warning: this code should not get compiled in unless ALL of
 341  * the following machine-dependent routines are supplied.
 342  *
 343  * Required machine dependent functions (not on a VAX):
 344  *     swapINX(i): save inexact flag and reset it to "i"
 345  *     swapENI(e): save inexact enable and reset it to "e"
 346  */
 347
 348 double
 349 drem(double x, double y)
 350 {
 351
 352 #ifdef national         /* order of words in floating point number */
 353         static const n0=3,n1=2,n2=1,n3=0;
 354 #else /* VAX, SUN, ZILOG, TAHOE */
 355         static const n0=0,n1=1,n2=2,n3=3;
 356 #endif
 357
 358         static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
 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
 449 newsqrt(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