1*f504f610SAugustin Cavalier #include <math.h> 2*f504f610SAugustin Cavalier #include <stdint.h> 3*f504f610SAugustin Cavalier remquof(float x,float y,int * quo)4*f504f610SAugustin Cavalierfloat remquof(float x, float y, int *quo) 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 int sx = ux.i>>31; 10*f504f610SAugustin Cavalier int sy = uy.i>>31; 11*f504f610SAugustin Cavalier uint32_t q; 12*f504f610SAugustin Cavalier uint32_t i; 13*f504f610SAugustin Cavalier uint32_t uxi = ux.i; 14*f504f610SAugustin Cavalier 15*f504f610SAugustin Cavalier *quo = 0; 16*f504f610SAugustin Cavalier if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) 17*f504f610SAugustin Cavalier return (x*y)/(x*y); 18*f504f610SAugustin Cavalier if (ux.i<<1 == 0) 19*f504f610SAugustin Cavalier return x; 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 q = 0; 38*f504f610SAugustin Cavalier if (ex < ey) { 39*f504f610SAugustin Cavalier if (ex+1 == ey) 40*f504f610SAugustin Cavalier goto end; 41*f504f610SAugustin Cavalier return x; 42*f504f610SAugustin Cavalier } 43*f504f610SAugustin Cavalier 44*f504f610SAugustin Cavalier /* x mod y */ 45*f504f610SAugustin Cavalier for (; ex > ey; ex--) { 46*f504f610SAugustin Cavalier i = uxi - uy.i; 47*f504f610SAugustin Cavalier if (i >> 31 == 0) { 48*f504f610SAugustin Cavalier uxi = i; 49*f504f610SAugustin Cavalier q++; 50*f504f610SAugustin Cavalier } 51*f504f610SAugustin Cavalier uxi <<= 1; 52*f504f610SAugustin Cavalier q <<= 1; 53*f504f610SAugustin Cavalier } 54*f504f610SAugustin Cavalier i = uxi - uy.i; 55*f504f610SAugustin Cavalier if (i >> 31 == 0) { 56*f504f610SAugustin Cavalier uxi = i; 57*f504f610SAugustin Cavalier q++; 58*f504f610SAugustin Cavalier } 59*f504f610SAugustin Cavalier if (uxi == 0) 60*f504f610SAugustin Cavalier ex = -30; 61*f504f610SAugustin Cavalier else 62*f504f610SAugustin Cavalier for (; uxi>>23 == 0; uxi <<= 1, ex--); 63*f504f610SAugustin Cavalier end: 64*f504f610SAugustin Cavalier /* scale result and decide between |x| and |x|-|y| */ 65*f504f610SAugustin Cavalier if (ex > 0) { 66*f504f610SAugustin Cavalier uxi -= 1U << 23; 67*f504f610SAugustin Cavalier uxi |= (uint32_t)ex << 23; 68*f504f610SAugustin Cavalier } else { 69*f504f610SAugustin Cavalier uxi >>= -ex + 1; 70*f504f610SAugustin Cavalier } 71*f504f610SAugustin Cavalier ux.i = uxi; 72*f504f610SAugustin Cavalier x = ux.f; 73*f504f610SAugustin Cavalier if (sy) 74*f504f610SAugustin Cavalier y = -y; 75*f504f610SAugustin Cavalier if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) { 76*f504f610SAugustin Cavalier x -= y; 77*f504f610SAugustin Cavalier q++; 78*f504f610SAugustin Cavalier } 79*f504f610SAugustin Cavalier q &= 0x7fffffff; 80*f504f610SAugustin Cavalier *quo = sx^sy ? -(int)q : (int)q; 81*f504f610SAugustin Cavalier return sx ? -x : x; 82*f504f610SAugustin Cavalier } 83