overby.me
this repo has no description
1/*
2 * fmod.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 * Modified for armLibm, removed code to set flags as this is a soft-float fallback only.
11 *
12 * Copyright 2007 Apple Inc. All rights reserved.
13 *
14 */
15
16#if defined(__SOFTFP__)
17
18 #include <math.h>
19 #include <float.h>
20 #include <stdint.h>
21
22static inline int local_ilogb( double x ) __attribute__ ((always_inline));
23static inline int local_ilogb( double x )
24{
25 union{ double d; uint64_t u;}u = {x};
26
27 u.u &= 0x7fffffffffffffffULL;
28 int32_t exp = (int32_t) (u.u >> 52);
29
30 if( __builtin_expect( (uint32_t) exp - 1U >= 2046, 0 ) )
31 { // +-denorm, the other possibilities: +0 +-inf, NaN are screend out by caller
32 u.u |= 0x3ff0000000000000ULL;
33 u.d -= 1.0;
34 exp = (int32_t) (u.u >> 52);
35
36 return exp - (1023+1022);
37 }
38
39 return exp - 1023;
40}
41
42static inline double local_scalbn( double x, int n );
43static inline double local_scalbn( double x, int n )
44{
45 union{ uint64_t u; double d;} u;
46
47 uint32_t absn = n >> (8*sizeof(n)-1);
48 absn = (n ^ absn) - absn;
49
50 if( absn > 1022 )
51 {
52 // step = copysign( 1022, n )
53 int signMask = n >> (8 * sizeof( int ) - 1);
54 int step = (1022 ^ signMask) - signMask;
55
56 u.u = ( (int64_t) step + 1023ULL) << 52;
57
58 if( n < 0 )
59 {
60 do
61 {
62 x *= u.d;
63 n -= step;
64 }while( n < -1022 );
65 }
66 else
67 {
68 do
69 {
70 x *= u.d;
71 n -= step;
72 }while( n > 1022 );
73 }
74 }
75
76 //create 2**n in double
77 u.u = ( (int64_t) n + 1023ULL) << 52;
78
79 //scale by appropriate power of 2
80 x *= u.d;
81
82 //return result
83 return x;
84}
85
86
87double fmod( double x, double y )
88{
89 union{ double d; uint64_t u;}ux = {x};
90 union{ double d; uint64_t u;}uy = {y};
91
92 uint64_t absx = ux.u & ~0x8000000000000000ULL;
93 uint64_t absy = uy.u & ~0x8000000000000000ULL;
94
95 if( absx-1ULL >= 0x7fefffffffffffffULL || absy-1ULL >= 0x7fefffffffffffffULL )
96 {
97 double fabsx = __builtin_fabs(x);
98 double fabsy = __builtin_fabs(y);
99
100 // deal with NaN
101 if( x != x || y != y )
102 return x + y;
103
104 //x = Inf or y = 0, return Invalid per IEEE-754
105 if( fabsx == __builtin_inf() || 0.0 == y )
106 {
107 return __builtin_nan("");
108 }
109
110 //handle trivial case
111 if( fabsy == __builtin_inf() || 0.0 == x )
112 {
113 return x;
114 }
115 }
116
117 if( absy >= absx )
118 {
119 if( absy == absx )
120 {
121 ux.u ^= absx;
122 return ux.d;
123 }
124
125 return x;
126 }
127
128 int32_t expx = absx >> 52;
129 int32_t expy = absy >> 52;
130 int64_t sx = absx & 0x000fffffffffffff;
131 int64_t sy = absy & 0x000fffffffffffff;
132
133 if( 0 == expx )
134 {
135 uint32_t shift = __builtin_clzll( absx ) - (64-53);
136 sx <<= shift;
137 expx = 1 - shift;
138 }
139
140 if( 0 == expy )
141 {
142 uint32_t shift = __builtin_clzll( absy ) - (64-53);
143 sy <<= shift;
144 expy = 1 - shift;
145 }
146 sx |= 0x0010000000000000ULL;
147 sy |= 0x0010000000000000ULL;
148
149
150 int32_t idiff = expx - expy;
151 int32_t shift = 0;
152 int64_t mask;
153
154 do
155 {
156 sx <<= shift;
157 idiff += ~shift;
158 sx -= sy;
159 mask = sx >> 63;
160 sx += sx;
161 sx += sy & mask;
162 shift = __builtin_clzll( sx ) - (64-53);
163 }
164 while( idiff >= shift && sx != 0LL );
165
166 if( idiff < 0 )
167 {
168 sx += sy & mask;
169 sx >>= 1;
170 idiff = 0;
171 }
172
173 sx <<= idiff;
174
175 if( 0 == sx )
176 {
177 ux.u &= 0x8000000000000000;
178 return ux.d;
179 }
180
181 shift = __builtin_clzll( sx ) - (64-53);
182 sx <<= shift;
183 expy -= shift;
184 sx &= 0x000fffffffffffffULL;
185 sx |= ux.u & 0x8000000000000000ULL;
186 if( expy > 0 )
187 {
188 ux.u = sx | ((int64_t) expy << 52);
189 return ux.d;
190 }
191
192 expy += 1022;
193 ux.u = sx | ((int64_t) expy << 52);
194 return ux.d * 0x1.0p-1022;
195
196}
197
198#endif
199