1 /* @(#)e_fmod.c 1.3 95/01/18 */
2 /*-
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13
14
15
16 #define BIAS (LDBL_MAX_EXP - 1)
17
18 /*
19 * These macros add and remove an explicit integer bit in front of the
20 * fractional mantissa, if the architecture doesn't have such a bit by
21 * default already.
22 */
23 #ifdef LDBL_IMPLICIT_NBIT
24 #define LDBL_NBIT 0
25 #define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE))
26 #define HFRAC_BITS (EXT_FRACHBITS + EXT_FRACHMBITS)
27 #else
28 #define LDBL_NBIT 0x80000000
29 #define SET_NBIT(hx) (hx)
30 #define HFRAC_BITS (EXT_FRACHBITS + EXT_FRACHMBITS - 1)
31 #endif
32
33 #define MANL_SHIFT (EXT_FRACLMBITS + EXT_FRACLBITS - 1)
34
35 static const long double Zero[] = {0.0L, -0.0L};
36
37 /*
38 * Return the IEEE remainder and set *quo to the last n bits of the
39 * quotient, rounded to the nearest integer. We choose n=31 because
40 * we wind up computing all the integer bits of the quotient anyway as
41 * a side-effect of computing the remainder by the shift and subtract
42 * method. In practice, this is far more bits than are needed to use
43 * remquo in reduction algorithms.
44 *
45 * Assumptions:
46 * - The low part of the mantissa fits in a manl_t exactly.
47 * - The high part of the mantissa fits in an int64_t with enough room
48 * for an explicit integer bit in front of the fractional bits.
49 */
50 long double
remquol(long double x,long double y,int * quo)51 remquol(long double x, long double y, int *quo)
52 {
53 int64_t hx,hz,hy,_hx;
54 uint64_t lx,ly,lz;
55 uint64_t sx,sxy;
56 int ix,iy,n,q;
57
58 GET_LDOUBLE_WORDS64(hx,lx,x);
59 GET_LDOUBLE_WORDS64(hy,ly,y);
60 sx = (hx>>48)&0x8000;
61 sxy = sx ^ ((hy>>48)&0x8000);
62 hx &= 0x7fffffffffffffffLL; /* |x| */
63 hy &= 0x7fffffffffffffffLL; /* |y| */
64 SET_LDOUBLE_WORDS64(x,hx,lx);
65 SET_LDOUBLE_WORDS64(y,hy,ly);
66
67 /* purge off exception values */
68 if((hy|ly)==0 || /* y=0 */
69 ((hx>>48) == BIAS + LDBL_MAX_EXP) || /* or x not finite */
70 ((hy>>48) == BIAS + LDBL_MAX_EXP &&
71 (((hy&0x0000ffffffffffffLL)&~LDBL_NBIT)|ly)!=0)) /* or y is NaN */
72 return (x*y)/(x*y);
73 if((hx>>48)<=(hy>>48)) {
74 if(((hx>>48)<(hy>>48)) ||
75 ((hx&0x0000ffffffffffffLL)<=(hy&0x0000ffffffffffffLL) &&
76 ((hx&0x0000ffffffffffffLL)<(hy&0x0000ffffffffffffLL) ||
77 lx<ly))) {
78 q = 0;
79 goto fixup; /* |x|<|y| return x or x-y */
80 }
81 if((hx&0x0000ffffffffffffLL)==(hy&0x0000ffffffffffffLL) &&
82 lx==ly) {
83 *quo = sxy? -1 : 1;
84 return Zero[sx!=0]; /* |x|=|y| return x*0*/
85 }
86 }
87
88 /* determine ix = ilogb(x) */
89 if((hx>>48) == 0) { /* subnormal x */
90 x *= 0x1.0p512L;
91 GET_LDOUBLE_WORDS64(hx,lx,x);
92 ix = (hx>>48) - (BIAS + 512);
93 } else {
94 ix = (hx>>48) - BIAS;
95 }
96
97 /* determine iy = ilogb(y) */
98 if((hy>>48) == 0) { /* subnormal y */
99 y *= 0x1.0p512L;
100 GET_LDOUBLE_WORDS64(hy,ly,y);
101 iy = (hy>>48) - (BIAS + 512);
102 } else {
103 iy = (hy>>48) - BIAS;
104 }
105
106 /* set up {hx,lx}, {hy,ly} and align y to x */
107 _hx = SET_NBIT(hx & 0x0000ffffffffffffLL);
108 hy = SET_NBIT(hy & 0x0000ffffffffffffLL);
109
110 /* fix point fmod */
111 n = ix - iy;
112 q = 0;
113
114 while(n--) {
115 hz=_hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
116 if(hz<0){_hx = _hx+_hx+(lx>>MANL_SHIFT); lx = lx+lx;}
117 else {_hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; q++;}
118 q <<= 1;
119 }
120 hz=_hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
121 if(hz>=0) {_hx=hz;lx=lz;q++;}
122
123 /* convert back to floating value and restore the sign */
124 if((_hx|lx)==0) { /* return sign(x)*0 */
125 *quo = (sxy ? -q : q);
126 return Zero[sx!=0];
127 }
128 while(_hx<(1LL<<HFRAC_BITS)) { /* normalize x */
129 _hx = _hx+_hx+(lx>>MANL_SHIFT); lx = lx+lx;
130 iy -= 1;
131 }
132 hx = (hx&0xffff000000000000LL) | (_hx&0x0000ffffffffffffLL);
133 if (iy < LDBL_MIN_EXP) {
134 hx = (hx&0x0000ffffffffffffLL) | (uint64_t)(iy + BIAS + 512)<<48;
135 SET_LDOUBLE_WORDS64(x,hx,lx);
136 x *= 0x1p-512L;
137 GET_LDOUBLE_WORDS64(hx,lx,x);
138 } else {
139 hx = (hx&0x0000ffffffffffffLL) | (uint64_t)(iy + BIAS)<<48;
140 }
141 hx &= 0x7fffffffffffffffLL;
142 SET_LDOUBLE_WORDS64(x,hx,lx);
143 fixup:
144 y = fabsl(y);
145 if (y < LDBL_MIN * 2) {
146 if (x+x>y || (x+x==y && (q & 1))) {
147 q++;
148 x-=y;
149 }
150 } else if (x>0.5L*y || (x==0.5L*y && (q & 1))) {
151 q++;
152 x-=y;
153 }
154
155 GET_LDOUBLE_MSW64(hx,x);
156 hx ^= (sx << 48);
157 SET_LDOUBLE_MSW64(x,hx);
158
159 q &= 0x7fffffff;
160 *quo = (sxy ? -q : q);
161 return x;
162 }
163
