1 #include "FEATURE/uwin" 2 3 #if !_UWIN || (_lib__copysign||_lib_copysign) && _lib_logb && (_lib__finite||_lib_finite) && (_lib_drem||_lib_remainder) && _lib_sqrt && _lib_ilogb && (_lib__scalb||_lib_scalb) 4 5 void _STUB_support(){} 6 7 #else 8 9 /* 10 * Copyright (c) 1985, 1993 11 * The Regents of the University of California. All rights reserved. 12 * 13 * Redistribution and use in source and binary forms, with or without 14 * modification, are permitted provided that the following conditions 15 * are met: 16 * 1. Redistributions of source code must retain the above copyright 17 * notice, this list of conditions and the following disclaimer. 18 * 2. Redistributions in binary form must reproduce the above copyright 19 * notice, this list of conditions and the following disclaimer in the 20 * documentation and/or other materials provided with the distribution. 21 * 3. Neither the name of the University nor the names of its contributors 22 * may be used to endorse or promote products derived from this software 23 * without specific prior written permission. 24 * 25 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 26 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 27 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 28 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 29 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 30 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 31 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 32 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 33 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 34 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 35 * SUCH DAMAGE. 36 */ 37 38 #ifndef lint 39 static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93"; 40 #endif /* not lint */ 41 42 /* 43 * Some IEEE standard 754 recommended functions and remainder and sqrt for 44 * supporting the C elementary functions. 45 ****************************************************************************** 46 * WARNING: 47 * These codes are developed (in double) to support the C elementary 48 * functions temporarily. They are not universal, and some of them are very 49 * slow (in particular, drem and sqrt is extremely inefficient). Each 50 * computer system should have its implementation of these functions using 51 * its own assembler. 52 ****************************************************************************** 53 * 54 * IEEE 754 required operations: 55 * drem(x,p) 56 * returns x REM y = x - [x/y]*y , where [x/y] is the integer 57 * nearest x/y; in half way case, choose the even one. 58 * sqrt(x) 59 * returns the square root of x correctly rounded according to 60 * the rounding mod. 61 * 62 * IEEE 754 recommended functions: 63 * (a) copysign(x,y) 64 * returns x with the sign of y. 65 * (b) scalb(x,N) 66 * returns x * (2**N), for integer values N. 67 * (c) logb(x) 68 * returns the unbiased exponent of x, a signed integer in 69 * double precision, except that logb(0) is -INF, logb(INF) 70 * is +INF, and logb(NAN) is that NAN. 71 * (d) finite(x) 72 * returns the value TRUE if -INF < x < +INF and returns 73 * FALSE otherwise. 74 * 75 * 76 * CODED IN C BY K.C. NG, 11/25/84; 77 * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85. 78 */ 79 80 #include "mathimpl.h" 81 82 #if defined(vax)||defined(tahoe) /* VAX D format */ 83 #include <errno.h> 84 static const unsigned short msign=0x7fff , mexp =0x7f80 ; 85 static const short prep1=57, gap=7, bias=129 ; 86 static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ; 87 #else /* defined(vax)||defined(tahoe) */ 88 static const unsigned short msign=0x7fff, mexp =0x7ff0 ; 89 static const short prep1=54, gap=4, bias=1023 ; 90 static const double novf=1.7E308, nunf=3.0E-308,zero=0.0; 91 #endif /* defined(vax)||defined(tahoe) */ 92 93 #if !_lib__scalb || !_lib_scalb 94 95 extern double _scalb(x,N) 96 double x; double N; 97 { 98 int k; 99 100 #ifdef national 101 unsigned short *px=(unsigned short *) &x + 3; 102 #else /* national */ 103 unsigned short *px=(unsigned short *) &x; 104 #endif /* national */ 105 106 if( x == zero ) return(x); 107 108 #if defined(vax)||defined(tahoe) 109 if( (k= *px & mexp ) != ~msign ) { 110 if (N < -260) 111 return(nunf*nunf); 112 else if (N > 260) { 113 return(copysign(infnan(ERANGE),x)); 114 } 115 #else /* defined(vax)||defined(tahoe) */ 116 if( (k= *px & mexp ) != mexp ) { 117 if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf); 118 if( k == 0 ) { 119 x *= scalb(1.0,prep1); N -= prep1; return(scalb(x,N));} 120 #endif /* defined(vax)||defined(tahoe) */ 121 122 if((k = (k>>gap)+ N) > 0 ) 123 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap); 124 else x=novf+novf; /* overflow */ 125 else 126 if( k > -prep1 ) 127 /* gradual underflow */ 128 {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);} 129 else 130 return(nunf*nunf); 131 } 132 return(x); 133 } 134 135 #endif 136 137 #if !_lib_scalb 138 139 extern double scalb(x,N) 140 double x; double N; 141 { 142 return _scalb(x, N); 143 } 144 145 #endif 146 147 #if !_lib__copysign 148 149 extern double _copysign(x,y) 150 double x,y; 151 { 152 #ifdef national 153 unsigned short *px=(unsigned short *) &x+3, 154 *py=(unsigned short *) &y+3; 155 #else /* national */ 156 unsigned short *px=(unsigned short *) &x, 157 *py=(unsigned short *) &y; 158 #endif /* national */ 159 160 #if defined(vax)||defined(tahoe) 161 if ( (*px & mexp) == 0 ) return(x); 162 #endif /* defined(vax)||defined(tahoe) */ 163 164 *px = ( *px & msign ) | ( *py & ~msign ); 165 return(x); 166 } 167 168 #endif 169 170 #if !_lib_copysign 171 172 extern double copysign(x,y) 173 double x,y; 174 { 175 return _copysign(x,y); 176 } 177 178 #endif 179 180 #if !_lib_logb 181 182 extern double logb(x) 183 double x; 184 { 185 186 #ifdef national 187 short *px=(short *) &x+3, k; 188 #else /* national */ 189 short *px=(short *) &x, k; 190 #endif /* national */ 191 192 #if defined(vax)||defined(tahoe) 193 return (int)(((*px&mexp)>>gap)-bias); 194 #else /* defined(vax)||defined(tahoe) */ 195 if( (k= *px & mexp ) != mexp ) 196 if ( k != 0 ) 197 return ( (k>>gap) - bias ); 198 else if( x != zero) 199 return ( -1022.0 ); 200 else 201 return(-(1.0/zero)); 202 else if(x != x) 203 return(x); 204 else 205 {*px &= msign; return(x);} 206 #endif /* defined(vax)||defined(tahoe) */ 207 } 208 209 #endif 210 211 #if !_lib__finite 212 213 extern int _finite(x) 214 double x; 215 { 216 #if defined(vax)||defined(tahoe) 217 return(1); 218 #else /* defined(vax)||defined(tahoe) */ 219 #ifdef national 220 return( (*((short *) &x+3 ) & mexp ) != mexp ); 221 #else /* national */ 222 return( (*((short *) &x ) & mexp ) != mexp ); 223 #endif /* national */ 224 #endif /* defined(vax)||defined(tahoe) */ 225 } 226 227 #endif 228 229 #if !_lib_finite 230 231 extern int finite(x) 232 double x; 233 { 234 return _finite(x); 235 } 236 237 #endif 238 239 #if !_lib_drem 240 241 extern double drem(x,p) 242 double x,p; 243 { 244 #if _lib_remainder 245 return remainder(x,p); 246 #else 247 short sign; 248 double hp,dp,tmp; 249 unsigned short k; 250 #ifdef national 251 unsigned short 252 *px=(unsigned short *) &x +3, 253 *pp=(unsigned short *) &p +3, 254 *pd=(unsigned short *) &dp +3, 255 *pt=(unsigned short *) &tmp+3; 256 #else /* national */ 257 unsigned short 258 *px=(unsigned short *) &x , 259 *pp=(unsigned short *) &p , 260 *pd=(unsigned short *) &dp , 261 *pt=(unsigned short *) &tmp; 262 #endif /* national */ 263 264 *pp &= msign ; 265 266 #if defined(vax)||defined(tahoe) 267 if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */ 268 #else /* defined(vax)||defined(tahoe) */ 269 if( ( *px & mexp ) == mexp ) 270 #endif /* defined(vax)||defined(tahoe) */ 271 return (x-p)-(x-p); /* create nan if x is inf */ 272 if (p == zero) { 273 #if defined(vax)||defined(tahoe) 274 return(infnan(EDOM)); 275 #else /* defined(vax)||defined(tahoe) */ 276 return zero/zero; 277 #endif /* defined(vax)||defined(tahoe) */ 278 } 279 280 #if defined(vax)||defined(tahoe) 281 if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */ 282 #else /* defined(vax)||defined(tahoe) */ 283 if( ( *pp & mexp ) == mexp ) 284 #endif /* defined(vax)||defined(tahoe) */ 285 { if (p != p) return p; else return x;} 286 287 else if ( ((*pp & mexp)>>gap) <= 1 ) 288 /* subnormal p, or almost subnormal p */ 289 { double b; b=scalb(1.0,(int)prep1); 290 p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} 291 else if ( p >= novf/2) 292 { p /= 2 ; x /= 2; return(drem(x,p)*2);} 293 else 294 { 295 dp=p+p; hp=p/2; 296 sign= *px & ~msign ; 297 *px &= msign ; 298 while ( x > dp ) 299 { 300 k=(*px & mexp) - (*pd & mexp) ; 301 tmp = dp ; 302 *pt += k ; 303 304 #if defined(vax)||defined(tahoe) 305 if( x < tmp ) *pt -= 128 ; 306 #else /* defined(vax)||defined(tahoe) */ 307 if( x < tmp ) *pt -= 16 ; 308 #endif /* defined(vax)||defined(tahoe) */ 309 310 x -= tmp ; 311 } 312 if ( x > hp ) 313 { x -= p ; if ( x >= hp ) x -= p ; } 314 315 #if defined(vax)||defined(tahoe) 316 if (x) 317 #endif /* defined(vax)||defined(tahoe) */ 318 *px ^= sign; 319 return( x); 320 321 } 322 #endif 323 } 324 325 #endif 326 327 #if !_lib_remainder 328 329 extern double remainder(x,p) 330 double x,p; 331 { 332 return drem(x,p); 333 } 334 335 #endif 336 337 #if !_lib_sqrt 338 339 extern double sqrt(x) 340 double x; 341 { 342 double q,s,b,r; 343 double t; 344 double const zero=0.0; 345 int m,n,i; 346 #if defined(vax)||defined(tahoe) 347 int k=54; 348 #else /* defined(vax)||defined(tahoe) */ 349 int k=51; 350 #endif /* defined(vax)||defined(tahoe) */ 351 352 /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ 353 if(x!=x||x==zero) return(x); 354 355 /* sqrt(negative) is invalid */ 356 if(x<zero) { 357 #if defined(vax)||defined(tahoe) 358 return (infnan(EDOM)); /* NaN */ 359 #else /* defined(vax)||defined(tahoe) */ 360 return(zero/zero); 361 #endif /* defined(vax)||defined(tahoe) */ 362 } 363 364 /* sqrt(INF) is INF */ 365 if(!finite(x)) return(x); 366 367 /* scale x to [1,4) */ 368 n=logb(x); 369 x=scalb(x,-n); 370 if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */ 371 m += n; 372 n = m/2; 373 if((n+n)!=m) {x *= 2; m -=1; n=m/2;} 374 375 /* generate sqrt(x) bit by bit (accumulating in q) */ 376 q=1.0; s=4.0; x -= 1.0; r=1; 377 for(i=1;i<=k;i++) { 378 t=s+1; x *= 4; r /= 2; 379 if(t<=x) { 380 s=t+t+2, x -= t; q += r;} 381 else 382 s *= 2; 383 } 384 385 /* generate the last bit and determine the final rounding */ 386 r/=2; x *= 4; 387 if(x==zero) goto end; 100+r; /* trigger inexact flag */ 388 if(s<x) { 389 q+=r; x -=s; s += 2; s *= 2; x *= 4; 390 t = (x-s)-5; 391 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */ 392 b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */ 393 if(t>=0) q+=r; } /* else: Round-to-nearest */ 394 else { 395 s *= 2; x *= 4; 396 t = (x-s)-1; 397 b=1.0+3*r/4; if(b==1.0) goto end; 398 b=1.0+r/4; if(b>1.0) t=1; 399 if(t>=0) q+=r; } 400 401 end: return(scalb(q,n)); 402 } 403 404 #endif 405 406 #if 0 407 /* DREM(X,Y) 408 * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) 409 * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 410 * INTENDED FOR ASSEMBLY LANGUAGE 411 * CODED IN C BY K.C. NG, 3/23/85, 4/8/85. 412 * 413 * Warning: this code should not get compiled in unless ALL of 414 * the following machine-dependent routines are supplied. 415 * 416 * Required machine dependent functions (not on a VAX): 417 * swapINX(i): save inexact flag and reset it to "i" 418 * swapENI(e): save inexact enable and reset it to "e" 419 */ 420 421 extern double drem(x,y) 422 double x,y; 423 { 424 425 #ifdef national /* order of words in floating point number */ 426 static const n0=3,n1=2,n2=1,n3=0; 427 #else /* VAX, SUN, ZILOG, TAHOE */ 428 static const n0=0,n1=1,n2=2,n3=3; 429 #endif 430 431 static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390; 432 static const double zero=0.0; 433 double hy,y1,t,t1; 434 short k; 435 long n; 436 int i,e; 437 unsigned short xexp,yexp, *px =(unsigned short *) &x , 438 nx,nf, *py =(unsigned short *) &y , 439 sign, *pt =(unsigned short *) &t , 440 *pt1 =(unsigned short *) &t1 ; 441 442 xexp = px[n0] & mexp ; /* exponent of x */ 443 yexp = py[n0] & mexp ; /* exponent of y */ 444 sign = px[n0] &0x8000; /* sign of x */ 445 446 /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */ 447 if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */ 448 if( xexp == mexp ) return(zero/zero); /* x is INF */ 449 if(y==zero) return(y/y); 450 451 /* save the inexact flag and inexact enable in i and e respectively 452 * and reset them to zero 453 */ 454 i=swapINX(0); e=swapENI(0); 455 456 /* subnormal number */ 457 nx=0; 458 if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;} 459 460 /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */ 461 if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} 462 463 nf=nx; 464 py[n0] &= 0x7fff; 465 px[n0] &= 0x7fff; 466 467 /* mask off the least significant 27 bits of y */ 468 t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t; 469 470 /* LOOP: argument reduction on x whenever x > y */ 471 loop: 472 while ( x > y ) 473 { 474 t=y; 475 t1=y1; 476 xexp=px[n0]&mexp; /* exponent of x */ 477 k=xexp-yexp-m25; 478 if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ 479 {pt[n0]+=k;pt1[n0]+=k;} 480 n=x/t; x=(x-n*t1)-n*(t-t1); 481 } 482 /* end while (x > y) */ 483 484 if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} 485 486 /* final adjustment */ 487 488 hy=y/2.0; 489 if(x>hy||((x==hy)&&n%2==1)) x-=y; 490 px[n0] ^= sign; 491 if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} 492 493 /* restore inexact flag and inexact enable */ 494 swapINX(i); swapENI(e); 495 496 return(x); 497 } 498 #endif 499 500 #if 0 501 /* SQRT 502 * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT 503 * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE 504 * CODED IN C BY K.C. NG, 3/22/85. 505 * 506 * Warning: this code should not get compiled in unless ALL of 507 * the following machine-dependent routines are supplied. 508 * 509 * Required machine dependent functions: 510 * swapINX(i) ...return the status of INEXACT flag and reset it to "i" 511 * swapRM(r) ...return the current Rounding Mode and reset it to "r" 512 * swapENI(e) ...return the status of inexact enable and reset it to "e" 513 * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer 514 * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer 515 */ 516 517 static const unsigned long table[] = { 518 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, 519 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, 520 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; 521 522 extern double newsqrt(x) 523 double x; 524 { 525 double y,z,t,addc(),subc() 526 double const b54=134217728.*134217728.; /* b54=2**54 */ 527 long mx,scalx; 528 long const mexp=0x7ff00000; 529 int i,j,r,e,swapINX(),swapRM(),swapENI(); 530 unsigned long *py=(unsigned long *) &y , 531 *pt=(unsigned long *) &t , 532 *px=(unsigned long *) &x ; 533 #ifdef national /* ordering of word in a floating point number */ 534 const int n0=1, n1=0; 535 #else 536 const int n0=0, n1=1; 537 #endif 538 /* Rounding Mode: RN ...round-to-nearest 539 * RZ ...round-towards 0 540 * RP ...round-towards +INF 541 * RM ...round-towards -INF 542 */ 543 const int RN=0,RZ=1,RP=2,RM=3; 544 /* machine dependent: work on a Zilog Z8070 545 * and a National 32081 & 16081 546 */ 547 548 /* exceptions */ 549 if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ 550 if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */ 551 if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */ 552 553 /* save, reset, initialize */ 554 e=swapENI(0); /* ...save and reset the inexact enable */ 555 i=swapINX(0); /* ...save INEXACT flag */ 556 r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */ 557 scalx=0; 558 559 /* subnormal number, scale up x to x*2**54 */ 560 if(mx==0) {x *= b54 ; scalx-=0x01b00000;} 561 562 /* scale x to avoid intermediate over/underflow: 563 * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */ 564 if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;} 565 if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;} 566 567 /* magic initial approximation to almost 8 sig. bits */ 568 py[n0]=(px[n0]>>1)+0x1ff80000; 569 py[n0]=py[n0]-table[(py[n0]>>15)&31]; 570 571 /* Heron's rule once with correction to improve y to almost 18 sig. bits */ 572 t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; 573 574 /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ 575 t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 576 t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; 577 578 /* twiddle last bit to force y correctly rounded */ 579 swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ 580 swapINX(0); /* ...clear INEXACT flag */ 581 swapENI(e); /* ...restore inexact enable status */ 582 t=x/y; /* ...chopped quotient, possibly inexact */ 583 j=swapINX(i); /* ...read and restore inexact flag */ 584 if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */ 585 b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */ 586 if(r==RN) t=addc(t); /* ...t=t+ulp */ 587 else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */ 588 y=y+t; /* ...chopped sum */ 589 py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */ 590 end: py[n0]=py[n0]+scalx; /* ...scale back y */ 591 swapRM(r); /* ...restore Rounding Mode */ 592 return(y); 593 } 594 #endif 595 596 #if !_lib_ilogb 597 598 extern int ilogb(double x) 599 { 600 return((int)logb(x)); 601 } 602 603 #endif 604 605 #endif 606