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