math: use the rounding idiom consistently
[oweals/musl.git] / src / math / remquol.c
1 #include "libm.h"
2
3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
4 long double remquol(long double x, long double y, int *quo)
5 {
6         return remquo(x, y, quo);
7 }
8 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
9 long double remquol(long double x, long double y, int *quo)
10 {
11         union ldshape ux = {x}, uy = {y};
12         int ex = ux.i.se & 0x7fff;
13         int ey = uy.i.se & 0x7fff;
14         int sx = ux.i.se >> 15;
15         int sy = uy.i.se >> 15;
16         uint32_t q;
17
18         *quo = 0;
19         if (y == 0 || isnan(y) || ex == 0x7fff)
20                 return (x*y)/(x*y);
21         if (x == 0)
22                 return x;
23
24         /* normalize x and y */
25         if (!ex) {
26                 ux.i.se = ex;
27                 ux.f *= 0x1p120f;
28                 ex = ux.i.se - 120;
29         }
30         if (!ey) {
31                 uy.i.se = ey;
32                 uy.f *= 0x1p120f;
33                 ey = uy.i.se - 120;
34         }
35
36         q = 0;
37         if (ex >= ey) {
38                 /* x mod y */
39 #if LDBL_MANT_DIG == 64
40                 uint64_t i, mx, my;
41                 mx = ux.i.m;
42                 my = uy.i.m;
43                 for (; ex > ey; ex--) {
44                         i = mx - my;
45                         if (mx >= my) {
46                                 mx = 2*i;
47                                 q++;
48                                 q <<= 1;
49                         } else if (2*mx < mx) {
50                                 mx = 2*mx - my;
51                                 q <<= 1;
52                                 q++;
53                         } else {
54                                 mx = 2*mx;
55                                 q <<= 1;
56                         }
57                 }
58                 i = mx - my;
59                 if (mx >= my) {
60                         mx = i;
61                         q++;
62                 }
63                 if (mx == 0)
64                         ex = -120;
65                 else
66                         for (; mx >> 63 == 0; mx *= 2, ex--);
67                 ux.i.m = mx;
68 #elif LDBL_MANT_DIG == 113
69                 uint64_t hi, lo, xhi, xlo, yhi, ylo;
70                 xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48;
71                 yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48;
72                 xlo = ux.i2.lo;
73                 ylo = ux.i2.lo;
74                 for (; ex > ey; ex--) {
75                         hi = xhi - yhi;
76                         lo = xlo - ylo;
77                         if (xlo < ylo)
78                                 hi -= 1;
79                         if (hi >> 63 == 0) {
80                                 xhi = 2*hi + (lo>>63);
81                                 xlo = 2*lo;
82                                 q++;
83                         } else {
84                                 xhi = 2*xhi + (xlo>>63);
85                                 xlo = 2*xlo;
86                         }
87                         q <<= 1;
88                 }
89                 hi = xhi - yhi;
90                 lo = xlo - ylo;
91                 if (xlo < ylo)
92                         hi -= 1;
93                 if (hi >> 63 == 0) {
94                         xhi = hi;
95                         xlo = lo;
96                         q++;
97                 }
98                 if ((xhi|xlo) == 0)
99                         ex = -120;
100                 else
101                         for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--);
102                 ux.i2.hi = xhi;
103                 ux.i2.lo = xlo;
104 #endif
105         }
106
107         /* scale result and decide between |x| and |x|-|y| */
108         if (ex <= 0) {
109                 ux.i.se = ex + 120;
110                 ux.f *= 0x1p-120f;
111         } else
112                 ux.i.se = ex;
113         x = ux.f;
114         if (sy)
115                 y = -y;
116         if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) {
117                 x -= y;
118                 q++;
119         }
120         q &= 0x7fffffff;
121         *quo = sx^sy ? -(int)q : (int)q;
122         return sx ? -x : x;
123 }
124 #endif