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 | /* |
---|
17 | * Modified for ALMOS-MKH OS at UPMC, France, August 2018. (Alain Greiner) |
---|
18 | */ |
---|
19 | |
---|
20 | |
---|
21 | #include "math.h" |
---|
22 | #include "math_private.h" |
---|
23 | |
---|
24 | double fmod(double x, double y) |
---|
25 | { |
---|
26 | union { |
---|
27 | double d; |
---|
28 | uint64_t u; |
---|
29 | } ux = {x}; |
---|
30 | union { |
---|
31 | double d; |
---|
32 | uint64_t u; |
---|
33 | } uy = {y}; |
---|
34 | |
---|
35 | uint64_t absx = ux.u & ~0x8000000000000000ULL; |
---|
36 | uint64_t absy = uy.u & ~0x8000000000000000ULL; |
---|
37 | |
---|
38 | if (absx - 1ULL >= 0x7fefffffffffffffULL || absy - 1ULL >= 0x7fefffffffffffffULL) { |
---|
39 | double fabsx = fabs(x); |
---|
40 | double fabsy = fabs(y); |
---|
41 | |
---|
42 | // deal with NaN |
---|
43 | if (x != x || y != y) { |
---|
44 | return x + y; |
---|
45 | } |
---|
46 | |
---|
47 | // x = Inf or y = 0, return Invalid per IEEE-754 |
---|
48 | if (!isfinite(fabsx) || 0.0 == y) { |
---|
49 | return 0.0 / 0.0; // NaN |
---|
50 | } |
---|
51 | |
---|
52 | // handle trivial case |
---|
53 | if (!isfinite(fabsy) || 0.0 == x) { |
---|
54 | return x; |
---|
55 | } |
---|
56 | } |
---|
57 | |
---|
58 | if (absy >= absx) { |
---|
59 | if (absy == absx) { |
---|
60 | ux.u ^= absx; |
---|
61 | return ux.d; |
---|
62 | } |
---|
63 | |
---|
64 | return x; |
---|
65 | } |
---|
66 | |
---|
67 | int32_t expx = absx >> 52; |
---|
68 | int32_t expy = absy >> 52; |
---|
69 | int64_t sx = absx & 0x000fffffffffffffULL; |
---|
70 | int64_t sy = absy & 0x000fffffffffffffULL; |
---|
71 | |
---|
72 | if (0 == expx) { |
---|
73 | uint32_t shift = __builtin_clzll(absx) - (64 - 53); |
---|
74 | sx <<= shift; |
---|
75 | expx = 1 - shift; |
---|
76 | } |
---|
77 | |
---|
78 | if (0 == expy) |
---|
79 | { |
---|
80 | uint32_t shift = __builtin_clzll(absy) - (64 - 53); |
---|
81 | sy <<= shift; |
---|
82 | expy = 1 - shift; |
---|
83 | } |
---|
84 | sx |= 0x0010000000000000ULL; |
---|
85 | sy |= 0x0010000000000000ULL; |
---|
86 | |
---|
87 | |
---|
88 | int32_t idiff = expx - expy; |
---|
89 | int32_t shift = 0; |
---|
90 | int64_t mask; |
---|
91 | |
---|
92 | do { |
---|
93 | sx <<= shift; |
---|
94 | idiff += ~shift; |
---|
95 | sx -= sy; |
---|
96 | mask = sx >> 63; |
---|
97 | sx += sx; |
---|
98 | sx += sy & mask; |
---|
99 | shift = __builtin_clzll(sx) - (64 - 53); |
---|
100 | } |
---|
101 | while (idiff >= shift && sx != 0LL); |
---|
102 | |
---|
103 | if (idiff < 0) { |
---|
104 | sx += sy & mask; |
---|
105 | sx >>= 1; |
---|
106 | idiff = 0; |
---|
107 | } |
---|
108 | |
---|
109 | sx <<= idiff; |
---|
110 | |
---|
111 | if (0 == sx) { |
---|
112 | ux.u &= 0x8000000000000000ULL; |
---|
113 | return ux.d; |
---|
114 | } |
---|
115 | |
---|
116 | shift = __builtin_clzll(sx) - (64 - 53); |
---|
117 | sx <<= shift; |
---|
118 | expy -= shift; |
---|
119 | sx &= 0x000fffffffffffffULL; |
---|
120 | sx |= ux.u & 0x8000000000000000ULL; |
---|
121 | if (expy > 0) { |
---|
122 | ux.u = sx | ((int64_t) expy << 52); |
---|
123 | return ux.d; |
---|
124 | } |
---|
125 | |
---|
126 | expy += 1022; |
---|
127 | ux.u = sx | ((int64_t) expy << 52); |
---|
128 | return ux.d * 0x1.0p-1022; |
---|
129 | |
---|
130 | } |
---|
131 | |
---|
132 | |
---|