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{
  6	return fmod(x, y);
  7}
  8#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
  9long 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