overby.me
this repo has no description
1/*
2 * remquo.c
3 * cLibm
4 *
5 * Written by Jon Okada, started on December 7th, 1992.
6 * Modified by Paul Finlayson (PAF) for MathLib v2.
7 * Modified by A. Sazegari (ali) for MathLib v3.
8 * Modified and ported by Robert A. Murley (ram) for Mac OS X.
9 * Modified for cLibm, fixed a few edge cases, rewrote local_ funcs by Ian Ollmann.
10 *
11 * Copyright 2007 Apple Inc. All rights reserved.
12 *
13 */
14
15 #include <math.h>
16 #include <float.h>
17 #include <stdint.h>
18
19static inline int local_ilogb( double x ) __attribute__ ((always_inline));
20static inline int local_ilogb( double x )
21{
22 union{ double d; uint64_t u;}u = {x};
23
24 u.u &= 0x7fffffffffffffffULL;
25 int32_t exp = (int32_t) (u.u >> 52);
26
27 if( __builtin_expect( (uint32_t) exp - 1U >= 2046, 0 ) )
28 { // +-denorm, the other possibilities: +0 +-inf, NaN are screend out by caller
29 u.u |= 0x3ff0000000000000ULL;
30 u.d -= 1.0;
31 exp = (int32_t) (u.u >> 52);
32
33 return exp - (1023+1022);
34 }
35
36 return exp - 1023;
37}
38
39static inline double local_scalbn( double x, int n );
40static inline double local_scalbn( double x, int n )
41{
42 union{ uint64_t u; double d;} u;
43
44 uint32_t absn = n >> (8*sizeof(n)-1);
45 absn = (n ^ absn) - absn;
46
47 if( absn > 1022 )
48 {
49 // step = copysign( 1022, n )
50 int signMask = n >> (8 * sizeof( int ) - 1);
51 int step = (1022 ^ signMask) - signMask;
52
53 u.u = ( (int64_t) step + 1023ULL) << 52;
54
55 if( n < 0 )
56 {
57 do
58 {
59 x *= u.d;
60 n -= step;
61 }while( n < -1022 );
62 }
63 else
64 {
65 do
66 {
67 x *= u.d;
68 n -= step;
69 }while( n > 1022 );
70 }
71 }
72
73 //create 2**n in double
74 u.u = ( (int64_t) n + 1023ULL) << 52;
75
76 //scale by appropriate power of 2
77 x *= u.d;
78
79 //return result
80 return x;
81}
82
83#ifdef ARMLIBM_SET_FLAGS
84#include "required_arithmetic.h"
85
86double remquo( double x, double y, int *quo )
87{
88
89 *quo = 0; // Init quotient result
90
91 // deal with NaN
92 if( x != x || y != y )
93 return x + y;
94
95 //x = Inf or y = 0, return Invalid per IEEE-754
96 double fabsx = __builtin_fabs(x);
97 if( fabsx == __builtin_inf() || 0.0 == y )
98 {
99 return required_add_double( __builtin_inf(), -__builtin_inf() ); //set invalid
100 }
101
102 //handle trivial case
103 double fabsy = __builtin_fabs(y);
104 if( fabsy == __builtin_inf() || 0.0 == x )
105 {
106 required_add_double( fabsx, 0.0 ); //signal underflow (that is, no flag set, but exception occurs if unmasked) if x is denorm
107 return x;
108 }
109
110 //we have to work
111 int32_t iquo = 0;
112 int32_t iscx = local_ilogb(fabsx);
113 int32_t iscy = local_ilogb(fabsy);
114 int32_t idiff = iscx - iscy;
115 double x1, y1, z;
116 x1 = fabsx;
117 y1 = fabsy;
118 if( idiff >= 0 )
119 {
120 int32_t i;
121
122 if( idiff )
123 {
124 y1 = local_scalbn( y1, -iscy );
125 x1 = local_scalbn( x1, -iscx );
126 for( i = idiff; i != 0; i-- )
127 {
128 if( x1 >= y1 )
129 {
130 x1 -= y1;
131 iquo += 1;
132 }
133 iquo += iquo;
134 x1 += x1;
135 }
136 x1 = local_scalbn( x1, iscy );
137 }
138 if( x1 >= fabsy )
139 {
140 x1 -= fabsy;
141 iquo += 1;
142 }
143 }
144
145 if( x1 < 0x1.0p1022 )
146 {
147 z = x1 + x1;
148 if( (z > fabsy) || ((z == fabsy) & ((iquo & 1) != 0 )))
149 {
150 x1 -= fabsy;
151 iquo += 1;
152 }
153 }
154 else
155 {
156 // x1 is only this large if y is even bigger, so we can safely divide y by 2
157 double halfy = 0.5*fabsy;
158 if( (x1 > halfy) || ((x1 == halfy) & ((iquo & 1) != 0 )))
159 {
160 x1 -= fabsy;
161 iquo += 1;
162 }
163 }
164
165
166 if( x < 0 )
167 x1 = -x1;
168 iquo &= 0x7f;
169
170 union{ double d[2]; uint64_t u[2];} u = { { x, y } };
171 int32_t sign = (int32_t) ( (u.u[0] ^ u.u[1]) >> 32 );
172 sign = sign >> 31;
173 *quo = (iquo ^ sign) - sign;
174
175 return x1;
176}
177
178#else
179
180double remquo( double x, double y, int *quo )
181{
182
183 *quo = 0; // Init quotient result
184
185 // deal with NaN
186 if( x != x || y != y )
187 return x + y;
188
189 //x = Inf or y = 0, return Invalid per IEEE-754
190 double fabsx = __builtin_fabs(x);
191 if( fabsx == __builtin_inf() || 0.0 == y )
192 {
193 return __builtin_nan("");
194 }
195
196 //handle trivial case
197 double fabsy = __builtin_fabs(y);
198 if( fabsy == __builtin_inf() || 0.0 == x )
199 {
200 return x;
201 }
202
203 //we have to work
204 int32_t iquo = 0;
205 int32_t iscx = local_ilogb(fabsx);
206 int32_t iscy = local_ilogb(fabsy);
207 int32_t idiff = iscx - iscy;
208 double x1, y1, z;
209 x1 = fabsx;
210 y1 = fabsy;
211 if( idiff >= 0 )
212 {
213 int32_t i;
214
215 if( idiff )
216 {
217 y1 = local_scalbn( y1, -iscy );
218 x1 = local_scalbn( x1, -iscx );
219 for( i = idiff; i != 0; i-- )
220 {
221 if( x1 >= y1 )
222 {
223 x1 -= y1;
224 iquo += 1;
225 }
226 iquo += iquo;
227 x1 += x1;
228 }
229 x1 = local_scalbn( x1, iscy );
230 }
231 if( x1 >= fabsy )
232 {
233 x1 -= fabsy;
234 iquo += 1;
235 }
236 }
237
238 if( x1 < 0x1.0p1022 )
239 {
240 z = x1 + x1;
241 if( (z > fabsy) || ((z == fabsy) & ((iquo & 1) != 0 )))
242 {
243 x1 -= fabsy;
244 iquo += 1;
245 }
246 }
247 else
248 {
249 // x1 is only this large if y is even bigger, so we can safely divide y by 2
250 double halfy = 0.5*fabsy;
251 if( (x1 > halfy) || ((x1 == halfy) & ((iquo & 1) != 0 )))
252 {
253 x1 -= fabsy;
254 iquo += 1;
255 }
256 }
257
258 if( x < 0 )
259 x1 = -x1;
260 iquo &= 0x7f;
261
262 union{ double d[2]; uint64_t u[2];} u = { { x, y } };
263 int32_t sign = (int32_t) ( (u.u[0] ^ u.u[1]) >> 32 );
264 sign = sign >> 31;
265 *quo = (iquo ^ sign) - sign;
266
267 return x1;
268}
269
270
271#endif // ARMLIBM_SET_FLAGS