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