1 /* $NetBSD: n_support.c,v 1.4 2002/06/15 00:10:18 matt 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. All advertising materials mentioning features or use of this software 15 * must display the following acknowledgement: 16 * This product includes software developed by the University of 17 * California, Berkeley and its contributors. 18 * 4. Neither the name of the University nor the names of its contributors 19 * may be used to endorse or promote products derived from this software 20 * without specific prior written permission. 21 * 22 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 23 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 25 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 26 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 27 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 28 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 29 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 31 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 32 * SUCH DAMAGE. 33 */ 34 35 #ifndef lint 36 static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93"; 37 #endif /* not lint */ 38 39 /* 40 * Some IEEE standard 754 recommended functions and remainder and sqrt for 41 * supporting the C elementary functions. 42 ****************************************************************************** 43 * WARNING: 44 * These codes are developed (in double) to support the C elementary 45 * functions temporarily. They are not universal, and some of them are very 46 * slow (in particular, drem and sqrt is extremely inefficient). Each 47 * computer system should have its implementation of these functions using 48 * its own assembler. 49 ****************************************************************************** 50 * 51 * IEEE 754 required operations: 52 * drem(x,p) 53 * returns x REM y = x - [x/y]*y , where [x/y] is the integer 54 * nearest x/y; in half way case, choose the even one. 55 * sqrt(x) 56 * returns the square root of x correctly rounded according to 57 * the rounding mod. 58 * 59 * IEEE 754 recommended functions: 60 * (a) copysign(x,y) 61 * returns x with the sign of y. 62 * (b) scalb(x,N) 63 * returns x * (2**N), for integer values N. 64 * (c) logb(x) 65 * returns the unbiased exponent of x, a signed integer in 66 * double precision, except that logb(0) is -INF, logb(INF) 67 * is +INF, and logb(NAN) is that NAN. 68 * (d) finite(x) 69 * returns the value TRUE if -INF < x < +INF and returns 70 * FALSE otherwise. 71 * 72 * 73 * CODED IN C BY K.C. NG, 11/25/84; 74 * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85. 75 */ 76 77 #include "mathimpl.h" 78 #include "trig.h" 79 80 #if defined(__vax__)||defined(tahoe) /* VAX D format */ 81 #include <errno.h> 82 static const unsigned short msign=0x7fff , mexp =0x7f80 ; 83 static const short prep1=57, gap=7, bias=129 ; 84 static const double novf=1.7E38, nunf=3.0E-39 ; 85 #else /* defined(__vax__)||defined(tahoe) */ 86 static const unsigned short msign=0x7fff, mexp =0x7ff0 ; 87 static const short prep1=54, gap=4, bias=1023 ; 88 static const double novf=1.7E308, nunf=3.0E-308; 89 #endif /* defined(__vax__)||defined(tahoe) */ 90 91 double 92 scalb(double x, int N) 93 { 94 int k; 95 96 #ifdef national 97 unsigned short *px=(unsigned short *) &x + 3; 98 #else /* national */ 99 unsigned short *px=(unsigned short *) &x; 100 #endif /* national */ 101 102 if( x == __zero ) return(x); 103 104 #if defined(__vax__)||defined(tahoe) 105 if( (k= *px & mexp ) != ~msign ) { 106 if (N < -260) 107 return(nunf*nunf); 108 else if (N > 260) { 109 return(copysign(infnan(ERANGE),x)); 110 } 111 #else /* defined(__vax__)||defined(tahoe) */ 112 if( (k= *px & mexp ) != mexp ) { 113 if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf); 114 if( k == 0 ) { 115 x *= scalb(1.0,(int)prep1); N -= prep1; return(scalb(x,N));} 116 #endif /* defined(__vax__)||defined(tahoe) */ 117 118 if((k = (k>>gap)+ N) > 0 ) 119 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap); 120 else x=novf+novf; /* overflow */ 121 else 122 if( k > -prep1 ) 123 /* gradual underflow */ 124 {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);} 125 else 126 return(nunf*nunf); 127 } 128 return(x); 129 } 130 131 132 double 133 copysign(double x, double y) 134 { 135 #ifdef national 136 unsigned short *px=(unsigned short *) &x+3, 137 *py=(unsigned short *) &y+3; 138 #else /* national */ 139 unsigned short *px=(unsigned short *) &x, 140 *py=(unsigned short *) &y; 141 #endif /* national */ 142 143 #if defined(__vax__)||defined(tahoe) 144 if ( (*px & mexp) == 0 ) return(x); 145 #endif /* defined(__vax__)||defined(tahoe) */ 146 147 *px = ( *px & msign ) | ( *py & ~msign ); 148 return(x); 149 } 150 151 double 152 logb(double x) 153 { 154 155 #ifdef national 156 short *px=(short *) &x+3, k; 157 #else /* national */ 158 short *px=(short *) &x, k; 159 #endif /* national */ 160 161 #if defined(__vax__)||defined(tahoe) 162 return (int)(((*px&mexp)>>gap)-bias); 163 #else /* defined(__vax__)||defined(tahoe) */ 164 if( (k= *px & mexp ) != mexp ) 165 if ( k != 0 ) 166 return ( (k>>gap) - bias ); 167 else if( x != __zero) 168 return ( -1022.0 ); 169 else 170 return(-(1.0/__zero)); 171 else if(x != x) 172 return(x); 173 else 174 {*px &= msign; return(x);} 175 #endif /* defined(__vax__)||defined(tahoe) */ 176 } 177 178 int 179 finite(double x) 180 { 181 #if defined(__vax__)||defined(tahoe) 182 return(1); 183 #else /* defined(__vax__)||defined(tahoe) */ 184 #ifdef national 185 return( (*((short *) &x+3 ) & mexp ) != mexp ); 186 #else /* national */ 187 return( (*((short *) &x ) & mexp ) != mexp ); 188 #endif /* national */ 189 #endif /* defined(__vax__)||defined(tahoe) */ 190 } 191 192 double 193 drem(double x, double p) 194 { 195 short sign; 196 double hp,dp,tmp; 197 unsigned short k; 198 #ifdef national 199 unsigned short 200 *px=(unsigned short *) &x +3, 201 *pp=(unsigned short *) &p +3, 202 *pd=(unsigned short *) &dp +3, 203 *pt=(unsigned short *) &tmp+3; 204 #else /* national */ 205 unsigned short 206 *px=(unsigned short *) &x , 207 *pp=(unsigned short *) &p , 208 *pd=(unsigned short *) &dp , 209 *pt=(unsigned short *) &tmp; 210 #endif /* national */ 211 212 *pp &= msign ; 213 214 #if defined(__vax__)||defined(tahoe) 215 if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */ 216 #else /* defined(__vax__)||defined(tahoe) */ 217 if( ( *px & mexp ) == mexp ) 218 #endif /* defined(__vax__)||defined(tahoe) */ 219 return (x-p)-(x-p); /* create nan if x is inf */ 220 if (p == __zero) { 221 #if defined(__vax__)||defined(tahoe) 222 return(infnan(EDOM)); 223 #else /* defined(__vax__)||defined(tahoe) */ 224 return __zero/__zero; 225 #endif /* defined(__vax__)||defined(tahoe) */ 226 } 227 228 #if defined(__vax__)||defined(tahoe) 229 if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */ 230 #else /* defined(__vax__)||defined(tahoe) */ 231 if( ( *pp & mexp ) == mexp ) 232 #endif /* defined(__vax__)||defined(tahoe) */ 233 { if (p != p) return p; else return x;} 234 235 else if ( ((*pp & mexp)>>gap) <= 1 ) 236 /* subnormal p, or almost subnormal p */ 237 { double b; b=scalb(1.0,(int)prep1); 238 p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} 239 else if ( p >= novf/2) 240 { p /= 2 ; x /= 2; return(drem(x,p)*2);} 241 else 242 { 243 dp=p+p; hp=p/2; 244 sign= *px & ~msign ; 245 *px &= msign ; 246 while ( x > dp ) 247 { 248 k=(*px & mexp) - (*pd & mexp) ; 249 tmp = dp ; 250 *pt += k ; 251 252 #if defined(__vax__)||defined(tahoe) 253 if( x < tmp ) *pt -= 128 ; 254 #else /* defined(__vax__)||defined(tahoe) */ 255 if( x < tmp ) *pt -= 16 ; 256 #endif /* defined(__vax__)||defined(tahoe) */ 257 258 x -= tmp ; 259 } 260 if ( x > hp ) 261 { x -= p ; if ( x >= hp ) x -= p ; } 262 263 #if defined(__vax__)||defined(tahoe) 264 if (x) 265 #endif /* defined(__vax__)||defined(tahoe) */ 266 *px ^= sign; 267 return( x); 268 269 } 270 } 271 272 273 double 274 sqrt(double x) 275 { 276 double q,s,b,r; 277 double t; 278 int m,n,i; 279 #if defined(__vax__)||defined(tahoe) 280 int k=54; 281 #else /* defined(__vax__)||defined(tahoe) */ 282 int k=51; 283 #endif /* defined(__vax__)||defined(tahoe) */ 284 285 /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ 286 if(x!=x||x==__zero) return(x); 287 288 /* sqrt(negative) is invalid */ 289 if(x<__zero) { 290 #if defined(__vax__)||defined(tahoe) 291 return (infnan(EDOM)); /* NaN */ 292 #else /* defined(__vax__)||defined(tahoe) */ 293 return(__zero/__zero); 294 #endif /* defined(__vax__)||defined(tahoe) */ 295 } 296 297 /* sqrt(INF) is INF */ 298 if(!finite(x)) return(x); 299 300 /* scale x to [1,4) */ 301 n=logb(x); 302 x=scalb(x,-n); 303 if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */ 304 m += n; 305 n = m/2; 306 if((n+n)!=m) {x *= 2; m -=1; n=m/2;} 307 308 /* generate sqrt(x) bit by bit (accumulating in q) */ 309 q=1.0; s=4.0; x -= 1.0; r=1; 310 for(i=1;i<=k;i++) { 311 t=s+1; x *= 4; r /= 2; 312 if(t<=x) { 313 s=t+t+2, x -= t; q += r;} 314 else 315 s *= 2; 316 } 317 318 /* generate the last bit and determine the final rounding */ 319 r/=2; x *= 4; 320 if(x==__zero) goto end; 100+r; /* trigger inexact flag */ 321 if(s<x) { 322 q+=r; x -=s; s += 2; s *= 2; x *= 4; 323 t = (x-s)-5; 324 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */ 325 b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */ 326 if(t>=0) q+=r; } /* else: Round-to-nearest */ 327 else { 328 s *= 2; x *= 4; 329 t = (x-s)-1; 330 b=1.0+3*r/4; if(b==1.0) goto end; 331 b=1.0+r/4; if(b>1.0) t=1; 332 if(t>=0) q+=r; } 333 334 end: return(scalb(q,n)); 335 } 336 337 #if 0 338 /* DREM(X,Y) 339 * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) 340 * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 341 * INTENDED FOR ASSEMBLY LANGUAGE 342 * CODED IN C BY K.C. NG, 3/23/85, 4/8/85. 343 * 344 * Warning: this code should not get compiled in unless ALL of 345 * the following machine-dependent routines are supplied. 346 * 347 * Required machine dependent functions (not on a VAX): 348 * swapINX(i): save inexact flag and reset it to "i" 349 * swapENI(e): save inexact enable and reset it to "e" 350 */ 351 352 double 353 drem(double x, double y) 354 { 355 356 #ifdef national /* order of words in floating point number */ 357 static const n0=3,n1=2,n2=1,n3=0; 358 #else /* VAX, SUN, ZILOG, TAHOE */ 359 static const n0=0,n1=1,n2=2,n3=3; 360 #endif 361 362 static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390; 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 453 newsqrt(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 525