github.com/afumu/libc@v0.0.6/musl/src/math/remquol.c (about)

     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