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. 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 32static 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 87double 88scalb(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 128double 129copysign(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 147double 148logb(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 174int 175finite(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 188double 189drem(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 269double 270sqrt(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 330end: 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 348double 349drem(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 */ 397loop: 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 443static const unsigned long table[] = { 4440, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, 44558733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, 44621581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; 447 448double 449newsqrt(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) */ 516end: py[n0]=py[n0]+scalx; /* ...scale back y */ 517 swapRM(r); /* ...restore Rounding Mode */ 518 return(y); 519} 520#endif 521