fmodl.c raw

   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