1*f504f610SAugustin Cavalier #include <math.h> 2*f504f610SAugustin Cavalier #include <stdint.h> 3*f504f610SAugustin Cavalier fmodf(float x,float y)4*f504f610SAugustin Cavalierfloat fmodf(float x, float y) 5*f504f610SAugustin Cavalier { 6*f504f610SAugustin Cavalier union {float f; uint32_t i;} ux = {x}, uy = {y}; 7*f504f610SAugustin Cavalier int ex = ux.i>>23 & 0xff; 8*f504f610SAugustin Cavalier int ey = uy.i>>23 & 0xff; 9*f504f610SAugustin Cavalier uint32_t sx = ux.i & 0x80000000; 10*f504f610SAugustin Cavalier uint32_t i; 11*f504f610SAugustin Cavalier uint32_t uxi = ux.i; 12*f504f610SAugustin Cavalier 13*f504f610SAugustin Cavalier if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) 14*f504f610SAugustin Cavalier return (x*y)/(x*y); 15*f504f610SAugustin Cavalier if (uxi<<1 <= uy.i<<1) { 16*f504f610SAugustin Cavalier if (uxi<<1 == uy.i<<1) 17*f504f610SAugustin Cavalier return 0*x; 18*f504f610SAugustin Cavalier return x; 19*f504f610SAugustin Cavalier } 20*f504f610SAugustin Cavalier 21*f504f610SAugustin Cavalier /* normalize x and y */ 22*f504f610SAugustin Cavalier if (!ex) { 23*f504f610SAugustin Cavalier for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1); 24*f504f610SAugustin Cavalier uxi <<= -ex + 1; 25*f504f610SAugustin Cavalier } else { 26*f504f610SAugustin Cavalier uxi &= -1U >> 9; 27*f504f610SAugustin Cavalier uxi |= 1U << 23; 28*f504f610SAugustin Cavalier } 29*f504f610SAugustin Cavalier if (!ey) { 30*f504f610SAugustin Cavalier for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1); 31*f504f610SAugustin Cavalier uy.i <<= -ey + 1; 32*f504f610SAugustin Cavalier } else { 33*f504f610SAugustin Cavalier uy.i &= -1U >> 9; 34*f504f610SAugustin Cavalier uy.i |= 1U << 23; 35*f504f610SAugustin Cavalier } 36*f504f610SAugustin Cavalier 37*f504f610SAugustin Cavalier /* x mod y */ 38*f504f610SAugustin Cavalier for (; ex > ey; ex--) { 39*f504f610SAugustin Cavalier i = uxi - uy.i; 40*f504f610SAugustin Cavalier if (i >> 31 == 0) { 41*f504f610SAugustin Cavalier if (i == 0) 42*f504f610SAugustin Cavalier return 0*x; 43*f504f610SAugustin Cavalier uxi = i; 44*f504f610SAugustin Cavalier } 45*f504f610SAugustin Cavalier uxi <<= 1; 46*f504f610SAugustin Cavalier } 47*f504f610SAugustin Cavalier i = uxi - uy.i; 48*f504f610SAugustin Cavalier if (i >> 31 == 0) { 49*f504f610SAugustin Cavalier if (i == 0) 50*f504f610SAugustin Cavalier return 0*x; 51*f504f610SAugustin Cavalier uxi = i; 52*f504f610SAugustin Cavalier } 53*f504f610SAugustin Cavalier for (; uxi>>23 == 0; uxi <<= 1, ex--); 54*f504f610SAugustin Cavalier 55*f504f610SAugustin Cavalier /* scale result up */ 56*f504f610SAugustin Cavalier if (ex > 0) { 57*f504f610SAugustin Cavalier uxi -= 1U << 23; 58*f504f610SAugustin Cavalier uxi |= (uint32_t)ex << 23; 59*f504f610SAugustin Cavalier } else { 60*f504f610SAugustin Cavalier uxi >>= -ex + 1; 61*f504f610SAugustin Cavalier } 62*f504f610SAugustin Cavalier uxi |= sx; 63*f504f610SAugustin Cavalier ux.i = uxi; 64*f504f610SAugustin Cavalier return ux.f; 65*f504f610SAugustin Cavalier } 66