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