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

     1  #include <math.h>
     2  #include <stdint.h>
     3  
     4  double fmod(double x, double y)
     5  {
     6  	union {double f; uint64_t i;} ux = {x}, uy = {y};
     7  	int ex = ux.i>>52 & 0x7ff;
     8  	int ey = uy.i>>52 & 0x7ff;
     9  	int sx = ux.i>>63;
    10  	uint64_t i;
    11  
    12  	/* in the followings uxi should be ux.i, but then gcc wrongly adds */
    13  	/* float load/store to inner loops ruining performance and code size */
    14  	uint64_t uxi = ux.i;
    15  
    16  	if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff)
    17  		return (x*y)/(x*y);
    18  	if (uxi<<1 <= uy.i<<1) {
    19  		if (uxi<<1 == uy.i<<1)
    20  			return 0*x;
    21  		return x;
    22  	}
    23  
    24  	/* normalize x and y */
    25  	if (!ex) {
    26  		for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1);
    27  		uxi <<= -ex + 1;
    28  	} else {
    29  		uxi &= -1ULL >> 12;
    30  		uxi |= 1ULL << 52;
    31  	}
    32  	if (!ey) {
    33  		for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1);
    34  		uy.i <<= -ey + 1;
    35  	} else {
    36  		uy.i &= -1ULL >> 12;
    37  		uy.i |= 1ULL << 52;
    38  	}
    39  
    40  	/* x mod y */
    41  	for (; ex > ey; ex--) {
    42  		i = uxi - uy.i;
    43  		if (i >> 63 == 0) {
    44  			if (i == 0)
    45  				return 0*x;
    46  			uxi = i;
    47  		}
    48  		uxi <<= 1;
    49  	}
    50  	i = uxi - uy.i;
    51  	if (i >> 63 == 0) {
    52  		if (i == 0)
    53  			return 0*x;
    54  		uxi = i;
    55  	}
    56  	for (; uxi>>52 == 0; uxi <<= 1, ex--);
    57  
    58  	/* scale result */
    59  	if (ex > 0) {
    60  		uxi -= 1ULL << 52;
    61  		uxi |= (uint64_t)ex << 52;
    62  	} else {
    63  		uxi >>= -ex + 1;
    64  	}
    65  	uxi |= (uint64_t)sx << 63;
    66  	ux.i = uxi;
    67  	return ux.f;
    68  }