OLD | NEW |
1 #include "libm.h" | 1 #include "libm.h" |
2 | 2 |
3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
4 long double fmodl(long double x, long double y) | 4 long double fmodl(long double x, long double y) { |
5 { | 5 return fmod(x, y); |
6 » return fmod(x, y); | |
7 } | 6 } |
8 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 | 7 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 |
9 long double fmodl(long double x, long double y) | 8 long double fmodl(long double x, long double y) { |
10 { | 9 union ldshape ux = {x}, uy = {y}; |
11 » union ldshape ux = {x}, uy = {y}; | 10 int ex = ux.i.se & 0x7fff; |
12 » int ex = ux.i.se & 0x7fff; | 11 int ey = uy.i.se & 0x7fff; |
13 » int ey = uy.i.se & 0x7fff; | 12 int sx = ux.i.se & 0x8000; |
14 » int sx = ux.i.se & 0x8000; | |
15 | 13 |
16 » if (y == 0 || isnan(y) || ex == 0x7fff) | 14 if (y == 0 || isnan(y) || ex == 0x7fff) |
17 » » return (x*y)/(x*y); | 15 return (x * y) / (x * y); |
18 » ux.i.se = ex; | 16 ux.i.se = ex; |
19 » uy.i.se = ey; | 17 uy.i.se = ey; |
20 » if (ux.f <= uy.f) { | 18 if (ux.f <= uy.f) { |
21 » » if (ux.f == uy.f) | 19 if (ux.f == uy.f) |
22 » » » return 0*x; | 20 return 0 * x; |
23 » » return x; | 21 return x; |
24 » } | 22 } |
25 | 23 |
26 » /* normalize x and y */ | 24 /* normalize x and y */ |
27 » if (!ex) { | 25 if (!ex) { |
28 » » ux.f *= 0x1p120f; | 26 ux.f *= 0x1p120f; |
29 » » ex = ux.i.se - 120; | 27 ex = ux.i.se - 120; |
30 » } | 28 } |
31 » if (!ey) { | 29 if (!ey) { |
32 » » uy.f *= 0x1p120f; | 30 uy.f *= 0x1p120f; |
33 » » ey = uy.i.se - 120; | 31 ey = uy.i.se - 120; |
34 » } | 32 } |
35 | 33 |
36 » /* x mod y */ | 34 /* x mod y */ |
37 #if LDBL_MANT_DIG == 64 | 35 #if LDBL_MANT_DIG == 64 |
38 » uint64_t i, mx, my; | 36 uint64_t i, mx, my; |
39 » mx = ux.i.m; | 37 mx = ux.i.m; |
40 » my = uy.i.m; | 38 my = uy.i.m; |
41 » for (; ex > ey; ex--) { | 39 for (; ex > ey; ex--) { |
42 » » i = mx - my; | 40 i = mx - my; |
43 » » if (mx >= my) { | 41 if (mx >= my) { |
44 » » » if (i == 0) | 42 if (i == 0) |
45 » » » » return 0*x; | 43 return 0 * x; |
46 » » » mx = 2*i; | 44 mx = 2 * i; |
47 » » } else if (2*mx < mx) { | 45 } else if (2 * mx < mx) { |
48 » » » mx = 2*mx - my; | 46 mx = 2 * mx - my; |
49 » » } else { | 47 } else { |
50 » » » mx = 2*mx; | 48 mx = 2 * mx; |
51 » » } | 49 } |
52 » } | 50 } |
53 » i = mx - my; | 51 i = mx - my; |
54 » if (mx >= my) { | 52 if (mx >= my) { |
55 » » if (i == 0) | 53 if (i == 0) |
56 » » » return 0*x; | 54 return 0 * x; |
57 » » mx = i; | 55 mx = i; |
58 » } | 56 } |
59 » for (; mx >> 63 == 0; mx *= 2, ex--); | 57 for (; mx >> 63 == 0; mx *= 2, ex--) |
60 » ux.i.m = mx; | 58 ; |
| 59 ux.i.m = mx; |
61 #elif LDBL_MANT_DIG == 113 | 60 #elif LDBL_MANT_DIG == 113 |
62 » uint64_t hi, lo, xhi, xlo, yhi, ylo; | 61 uint64_t hi, lo, xhi, xlo, yhi, ylo; |
63 » xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48; | 62 xhi = (ux.i2.hi & -1ULL >> 16) | 1ULL << 48; |
64 » yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48; | 63 yhi = (uy.i2.hi & -1ULL >> 16) | 1ULL << 48; |
65 » xlo = ux.i2.lo; | 64 xlo = ux.i2.lo; |
66 » ylo = uy.i2.lo; | 65 ylo = uy.i2.lo; |
67 » for (; ex > ey; ex--) { | 66 for (; ex > ey; ex--) { |
68 » » hi = xhi - yhi; | 67 hi = xhi - yhi; |
69 » » lo = xlo - ylo; | 68 lo = xlo - ylo; |
70 » » if (xlo < ylo) | 69 if (xlo < ylo) |
71 » » » hi -= 1; | 70 hi -= 1; |
72 » » if (hi >> 63 == 0) { | 71 if (hi >> 63 == 0) { |
73 » » » if ((hi|lo) == 0) | 72 if ((hi | lo) == 0) |
74 » » » » return 0*x; | 73 return 0 * x; |
75 » » » xhi = 2*hi + (lo>>63); | 74 xhi = 2 * hi + (lo >> 63); |
76 » » » xlo = 2*lo; | 75 xlo = 2 * lo; |
77 » » } else { | 76 } else { |
78 » » » xhi = 2*xhi + (xlo>>63); | 77 xhi = 2 * xhi + (xlo >> 63); |
79 » » » xlo = 2*xlo; | 78 xlo = 2 * xlo; |
80 » » } | 79 } |
81 » } | 80 } |
82 » hi = xhi - yhi; | 81 hi = xhi - yhi; |
83 » lo = xlo - ylo; | 82 lo = xlo - ylo; |
84 » if (xlo < ylo) | 83 if (xlo < ylo) |
85 » » hi -= 1; | 84 hi -= 1; |
86 » if (hi >> 63 == 0) { | 85 if (hi >> 63 == 0) { |
87 » » if ((hi|lo) == 0) | 86 if ((hi | lo) == 0) |
88 » » » return 0*x; | 87 return 0 * x; |
89 » » xhi = hi; | 88 xhi = hi; |
90 » » xlo = lo; | 89 xlo = lo; |
91 » } | 90 } |
92 » for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--); | 91 for (; xhi >> 48 == 0; xhi = 2 * xhi + (xlo >> 63), xlo = 2 * xlo, ex--) |
93 » ux.i2.hi = xhi; | 92 ; |
94 » ux.i2.lo = xlo; | 93 ux.i2.hi = xhi; |
| 94 ux.i2.lo = xlo; |
95 #endif | 95 #endif |
96 | 96 |
97 » /* scale result */ | 97 /* scale result */ |
98 » if (ex <= 0) { | 98 if (ex <= 0) { |
99 » » ux.i.se = (ex+120)|sx; | 99 ux.i.se = (ex + 120) | sx; |
100 » » ux.f *= 0x1p-120f; | 100 ux.f *= 0x1p-120f; |
101 » } else | 101 } else |
102 » » ux.i.se = ex|sx; | 102 ux.i.se = ex | sx; |
103 » return ux.f; | 103 return ux.f; |
104 } | 104 } |
105 #endif | 105 #endif |
OLD | NEW |