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