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