remquol.c raw

   1  #include "libm.h"
   2  
   3  #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
   4  long double remquol(long double x, long double y, int *quo)
   5  {
   6  	return remquo(x, y, quo);
   7  }
   8  #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
   9  long double remquol(long double x, long double y, int *quo)
  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 >> 15;
  15  	int sy = uy.i.se >> 15;
  16  	uint32_t q;
  17  
  18  	*quo = 0;
  19  	if (y == 0 || isnan(y) || ex == 0x7fff)
  20  		return (x*y)/(x*y);
  21  	if (x == 0)
  22  		return x;
  23  
  24  	/* normalize x and y */
  25  	if (!ex) {
  26  		ux.i.se = ex;
  27  		ux.f *= 0x1p120f;
  28  		ex = ux.i.se - 120;
  29  	}
  30  	if (!ey) {
  31  		uy.i.se = ey;
  32  		uy.f *= 0x1p120f;
  33  		ey = uy.i.se - 120;
  34  	}
  35  
  36  	q = 0;
  37  	if (ex >= ey) {
  38  		/* x mod y */
  39  #if LDBL_MANT_DIG == 64
  40  		uint64_t i, mx, my;
  41  		mx = ux.i.m;
  42  		my = uy.i.m;
  43  		for (; ex > ey; ex--) {
  44  			i = mx - my;
  45  			if (mx >= my) {
  46  				mx = 2*i;
  47  				q++;
  48  				q <<= 1;
  49  			} else if (2*mx < mx) {
  50  				mx = 2*mx - my;
  51  				q <<= 1;
  52  				q++;
  53  			} else {
  54  				mx = 2*mx;
  55  				q <<= 1;
  56  			}
  57  		}
  58  		i = mx - my;
  59  		if (mx >= my) {
  60  			mx = i;
  61  			q++;
  62  		}
  63  		if (mx == 0)
  64  			ex = -120;
  65  		else
  66  			for (; mx >> 63 == 0; mx *= 2, ex--);
  67  		ux.i.m = mx;
  68  #elif LDBL_MANT_DIG == 113
  69  		uint64_t hi, lo, xhi, xlo, yhi, ylo;
  70  		xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48;
  71  		yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48;
  72  		xlo = ux.i2.lo;
  73  		ylo = ux.i2.lo;
  74  		for (; ex > ey; ex--) {
  75  			hi = xhi - yhi;
  76  			lo = xlo - ylo;
  77  			if (xlo < ylo)
  78  				hi -= 1;
  79  			if (hi >> 63 == 0) {
  80  				xhi = 2*hi + (lo>>63);
  81  				xlo = 2*lo;
  82  				q++;
  83  			} else {
  84  				xhi = 2*xhi + (xlo>>63);
  85  				xlo = 2*xlo;
  86  			}
  87  			q <<= 1;
  88  		}
  89  		hi = xhi - yhi;
  90  		lo = xlo - ylo;
  91  		if (xlo < ylo)
  92  			hi -= 1;
  93  		if (hi >> 63 == 0) {
  94  			xhi = hi;
  95  			xlo = lo;
  96  			q++;
  97  		}
  98  		if ((xhi|xlo) == 0)
  99  			ex = -120;
 100  		else
 101  			for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--);
 102  		ux.i2.hi = xhi;
 103  		ux.i2.lo = xlo;
 104  #endif
 105  	}
 106  
 107  	/* scale result and decide between |x| and |x|-|y| */
 108  	if (ex <= 0) {
 109  		ux.i.se = ex + 120;
 110  		ux.f *= 0x1p-120f;
 111  	} else
 112  		ux.i.se = ex;
 113  	x = ux.f;
 114  	if (sy)
 115  		y = -y;
 116  	if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) {
 117  		x -= y;
 118  		q++;
 119  	}
 120  	q &= 0x7fffffff;
 121  	*quo = sx^sy ? -(int)q : (int)q;
 122  	return sx ? -x : x;
 123  }
 124  #endif
 125