| 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 |