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