1 #include "libm.h"
2
3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
remquol(long double x,long double y,int * quo)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
remquol(long double x,long double y,int * quo)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
125