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
_STUB_support()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
ilogb(double x)598 extern int ilogb(double x)
599 {
600 return((int)logb(x));
601 }
602
603 #endif
604
605 #endif
606