1 /* $NetBSD: n_support.c,v 1.3 1999/07/02 15:37:37 simonb 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 79 #if defined(__vax__)||defined(tahoe) /* VAX D format */ 80 #include <errno.h> 81 static const unsigned short msign=0x7fff , mexp =0x7f80 ; 82 static const short prep1=57, gap=7, bias=129 ; 83 static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ; 84 #else /* defined(__vax__)||defined(tahoe) */ 85 static const unsigned short msign=0x7fff, mexp =0x7ff0 ; 86 static const short prep1=54, gap=4, bias=1023 ; 87 static const double novf=1.7E308, nunf=3.0E-308,zero=0.0; 88 #endif /* defined(__vax__)||defined(tahoe) */ 89 90 double scalb(x,N) 91 double x; int N; 92 { 93 int k; 94 95 #ifdef national 96 unsigned short *px=(unsigned short *) &x + 3; 97 #else /* national */ 98 unsigned short *px=(unsigned short *) &x; 99 #endif /* national */ 100 101 if( x == zero ) return(x); 102 103 #if defined(__vax__)||defined(tahoe) 104 if( (k= *px & mexp ) != ~msign ) { 105 if (N < -260) 106 return(nunf*nunf); 107 else if (N > 260) { 108 return(copysign(infnan(ERANGE),x)); 109 } 110 #else /* defined(__vax__)||defined(tahoe) */ 111 if( (k= *px & mexp ) != mexp ) { 112 if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf); 113 if( k == 0 ) { 114 x *= scalb(1.0,(int)prep1); N -= prep1; return(scalb(x,N));} 115 #endif /* defined(__vax__)||defined(tahoe) */ 116 117 if((k = (k>>gap)+ N) > 0 ) 118 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap); 119 else x=novf+novf; /* overflow */ 120 else 121 if( k > -prep1 ) 122 /* gradual underflow */ 123 {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);} 124 else 125 return(nunf*nunf); 126 } 127 return(x); 128 } 129 130 131 double copysign(x,y) 132 double x,y; 133 { 134 #ifdef national 135 unsigned short *px=(unsigned short *) &x+3, 136 *py=(unsigned short *) &y+3; 137 #else /* national */ 138 unsigned short *px=(unsigned short *) &x, 139 *py=(unsigned short *) &y; 140 #endif /* national */ 141 142 #if defined(__vax__)||defined(tahoe) 143 if ( (*px & mexp) == 0 ) return(x); 144 #endif /* defined(__vax__)||defined(tahoe) */ 145 146 *px = ( *px & msign ) | ( *py & ~msign ); 147 return(x); 148 } 149 150 double logb(x) 151 double x; 152 { 153 154 #ifdef national 155 short *px=(short *) &x+3, k; 156 #else /* national */ 157 short *px=(short *) &x, k; 158 #endif /* national */ 159 160 #if defined(__vax__)||defined(tahoe) 161 return (int)(((*px&mexp)>>gap)-bias); 162 #else /* defined(__vax__)||defined(tahoe) */ 163 if( (k= *px & mexp ) != mexp ) 164 if ( k != 0 ) 165 return ( (k>>gap) - bias ); 166 else if( x != zero) 167 return ( -1022.0 ); 168 else 169 return(-(1.0/zero)); 170 else if(x != x) 171 return(x); 172 else 173 {*px &= msign; return(x);} 174 #endif /* defined(__vax__)||defined(tahoe) */ 175 } 176 177 finite(x) 178 double x; 179 { 180 #if defined(__vax__)||defined(tahoe) 181 return(1); 182 #else /* defined(__vax__)||defined(tahoe) */ 183 #ifdef national 184 return( (*((short *) &x+3 ) & mexp ) != mexp ); 185 #else /* national */ 186 return( (*((short *) &x ) & mexp ) != mexp ); 187 #endif /* national */ 188 #endif /* defined(__vax__)||defined(tahoe) */ 189 } 190 191 double drem(x,p) 192 double x,p; 193 { 194 short sign; 195 double hp,dp,tmp; 196 unsigned short k; 197 #ifdef national 198 unsigned short 199 *px=(unsigned short *) &x +3, 200 *pp=(unsigned short *) &p +3, 201 *pd=(unsigned short *) &dp +3, 202 *pt=(unsigned short *) &tmp+3; 203 #else /* national */ 204 unsigned short 205 *px=(unsigned short *) &x , 206 *pp=(unsigned short *) &p , 207 *pd=(unsigned short *) &dp , 208 *pt=(unsigned short *) &tmp; 209 #endif /* national */ 210 211 *pp &= msign ; 212 213 #if defined(__vax__)||defined(tahoe) 214 if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */ 215 #else /* defined(__vax__)||defined(tahoe) */ 216 if( ( *px & mexp ) == mexp ) 217 #endif /* defined(__vax__)||defined(tahoe) */ 218 return (x-p)-(x-p); /* create nan if x is inf */ 219 if (p == zero) { 220 #if defined(__vax__)||defined(tahoe) 221 return(infnan(EDOM)); 222 #else /* defined(__vax__)||defined(tahoe) */ 223 return zero/zero; 224 #endif /* defined(__vax__)||defined(tahoe) */ 225 } 226 227 #if defined(__vax__)||defined(tahoe) 228 if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */ 229 #else /* defined(__vax__)||defined(tahoe) */ 230 if( ( *pp & mexp ) == mexp ) 231 #endif /* defined(__vax__)||defined(tahoe) */ 232 { if (p != p) return p; else return x;} 233 234 else if ( ((*pp & mexp)>>gap) <= 1 ) 235 /* subnormal p, or almost subnormal p */ 236 { double b; b=scalb(1.0,(int)prep1); 237 p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} 238 else if ( p >= novf/2) 239 { p /= 2 ; x /= 2; return(drem(x,p)*2);} 240 else 241 { 242 dp=p+p; hp=p/2; 243 sign= *px & ~msign ; 244 *px &= msign ; 245 while ( x > dp ) 246 { 247 k=(*px & mexp) - (*pd & mexp) ; 248 tmp = dp ; 249 *pt += k ; 250 251 #if defined(__vax__)||defined(tahoe) 252 if( x < tmp ) *pt -= 128 ; 253 #else /* defined(__vax__)||defined(tahoe) */ 254 if( x < tmp ) *pt -= 16 ; 255 #endif /* defined(__vax__)||defined(tahoe) */ 256 257 x -= tmp ; 258 } 259 if ( x > hp ) 260 { x -= p ; if ( x >= hp ) x -= p ; } 261 262 #if defined(__vax__)||defined(tahoe) 263 if (x) 264 #endif /* defined(__vax__)||defined(tahoe) */ 265 *px ^= sign; 266 return( x); 267 268 } 269 } 270 271 272 double sqrt(x) 273 double x; 274 { 275 double q,s,b,r; 276 double t; 277 double const zero=0.0; 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 drem(x,y) 353 double x,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 static const double zero=0.0; 364 double hy,y1,t,t1; 365 short k; 366 long n; 367 int i,e; 368 unsigned short xexp,yexp, *px =(unsigned short *) &x , 369 nx,nf, *py =(unsigned short *) &y , 370 sign, *pt =(unsigned short *) &t , 371 *pt1 =(unsigned short *) &t1 ; 372 373 xexp = px[n0] & mexp ; /* exponent of x */ 374 yexp = py[n0] & mexp ; /* exponent of y */ 375 sign = px[n0] &0x8000; /* sign of x */ 376 377 /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */ 378 if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */ 379 if( xexp == mexp ) return(zero/zero); /* x is INF */ 380 if(y==zero) return(y/y); 381 382 /* save the inexact flag and inexact enable in i and e respectively 383 * and reset them to zero 384 */ 385 i=swapINX(0); e=swapENI(0); 386 387 /* subnormal number */ 388 nx=0; 389 if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;} 390 391 /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */ 392 if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} 393 394 nf=nx; 395 py[n0] &= 0x7fff; 396 px[n0] &= 0x7fff; 397 398 /* mask off the least significant 27 bits of y */ 399 t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t; 400 401 /* LOOP: argument reduction on x whenever x > y */ 402 loop: 403 while ( x > y ) 404 { 405 t=y; 406 t1=y1; 407 xexp=px[n0]&mexp; /* exponent of x */ 408 k=xexp-yexp-m25; 409 if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ 410 {pt[n0]+=k;pt1[n0]+=k;} 411 n=x/t; x=(x-n*t1)-n*(t-t1); 412 } 413 /* end while (x > y) */ 414 415 if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} 416 417 /* final adjustment */ 418 419 hy=y/2.0; 420 if(x>hy||((x==hy)&&n%2==1)) x-=y; 421 px[n0] ^= sign; 422 if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} 423 424 /* restore inexact flag and inexact enable */ 425 swapINX(i); swapENI(e); 426 427 return(x); 428 } 429 #endif 430 431 #if 0 432 /* SQRT 433 * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT 434 * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE 435 * CODED IN C BY K.C. NG, 3/22/85. 436 * 437 * Warning: this code should not get compiled in unless ALL of 438 * the following machine-dependent routines are supplied. 439 * 440 * Required machine dependent functions: 441 * swapINX(i) ...return the status of INEXACT flag and reset it to "i" 442 * swapRM(r) ...return the current Rounding Mode and reset it to "r" 443 * swapENI(e) ...return the status of inexact enable and reset it to "e" 444 * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer 445 * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer 446 */ 447 448 static const unsigned long table[] = { 449 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, 450 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, 451 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; 452 453 double newsqrt(x) 454 double x; 455 { 456 double y,z,t,addc(),subc() 457 double const b54=134217728.*134217728.; /* b54=2**54 */ 458 long mx,scalx; 459 long const mexp=0x7ff00000; 460 int i,j,r,e,swapINX(),swapRM(),swapENI(); 461 unsigned long *py=(unsigned long *) &y , 462 *pt=(unsigned long *) &t , 463 *px=(unsigned long *) &x ; 464 #ifdef national /* ordering of word in a floating point number */ 465 const int n0=1, n1=0; 466 #else 467 const int n0=0, n1=1; 468 #endif 469 /* Rounding Mode: RN ...round-to-nearest 470 * RZ ...round-towards 0 471 * RP ...round-towards +INF 472 * RM ...round-towards -INF 473 */ 474 const int RN=0,RZ=1,RP=2,RM=3; 475 /* machine dependent: work on a Zilog Z8070 476 * and a National 32081 & 16081 477 */ 478 479 /* exceptions */ 480 if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ 481 if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */ 482 if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */ 483 484 /* save, reset, initialize */ 485 e=swapENI(0); /* ...save and reset the inexact enable */ 486 i=swapINX(0); /* ...save INEXACT flag */ 487 r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */ 488 scalx=0; 489 490 /* subnormal number, scale up x to x*2**54 */ 491 if(mx==0) {x *= b54 ; scalx-=0x01b00000;} 492 493 /* scale x to avoid intermediate over/underflow: 494 * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */ 495 if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;} 496 if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;} 497 498 /* magic initial approximation to almost 8 sig. bits */ 499 py[n0]=(px[n0]>>1)+0x1ff80000; 500 py[n0]=py[n0]-table[(py[n0]>>15)&31]; 501 502 /* Heron's rule once with correction to improve y to almost 18 sig. bits */ 503 t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; 504 505 /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ 506 t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 507 t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; 508 509 /* twiddle last bit to force y correctly rounded */ 510 swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ 511 swapINX(0); /* ...clear INEXACT flag */ 512 swapENI(e); /* ...restore inexact enable status */ 513 t=x/y; /* ...chopped quotient, possibly inexact */ 514 j=swapINX(i); /* ...read and restore inexact flag */ 515 if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */ 516 b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */ 517 if(r==RN) t=addc(t); /* ...t=t+ulp */ 518 else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */ 519 y=y+t; /* ...chopped sum */ 520 py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */ 521 end: py[n0]=py[n0]+scalx; /* ...scale back y */ 522 swapRM(r); /* ...restore Rounding Mode */ 523 return(y); 524 } 525 #endif 526