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

     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