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