1/* Compute remainder and a congruent to the quotient.
2 Copyright (C) 1997-2022 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
9
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <https://www.gnu.org/licenses/>. */
18
19#include <math.h>
20
21#include <math_private.h>
22#include <math_ldbl_opt.h>
23
24
25static const long double zero = 0.0;
26
27
28long double
29__remquol (long double x, long double y, int *quo)
30{
31 int64_t hx,hy;
32 uint64_t sx,lx,ly,qs;
33 int cquo;
34 double xhi, xlo, yhi, ylo;
35
36 ldbl_unpack (x, &xhi, &xlo);
37 EXTRACT_WORDS64 (hx, xhi);
38 EXTRACT_WORDS64 (lx, xlo);
39 ldbl_unpack (y, &yhi, &ylo);
40 EXTRACT_WORDS64 (hy, yhi);
41 EXTRACT_WORDS64 (ly, ylo);
42 sx = hx & 0x8000000000000000ULL;
43 qs = sx ^ (hy & 0x8000000000000000ULL);
44 ly ^= hy & 0x8000000000000000ULL;
45 hy &= 0x7fffffffffffffffLL;
46 lx ^= sx;
47 hx &= 0x7fffffffffffffffLL;
48
49 /* Purge off exception values. */
50 if (hy == 0)
51 return (x * y) / (x * y); /* y = 0 */
52 if ((hx >= 0x7ff0000000000000LL) /* x not finite */
53 || (hy > 0x7ff0000000000000LL)) /* y is NaN */
54 return (x * y) / (x * y);
55
56 if (hy <= 0x7fbfffffffffffffLL)
57 x = __ieee754_fmodl (x, 8 * y); /* now x < 8y */
58
59 if (((hx - hy) | (lx - ly)) == 0)
60 {
61 *quo = qs ? -1 : 1;
62 return zero * x;
63 }
64
65 x = fabsl (x: x);
66 y = fabsl (x: y);
67 cquo = 0;
68
69 if (hy <= 0x7fcfffffffffffffLL && x >= 4 * y)
70 {
71 x -= 4 * y;
72 cquo += 4;
73 }
74 if (hy <= 0x7fdfffffffffffffLL && x >= 2 * y)
75 {
76 x -= 2 * y;
77 cquo += 2;
78 }
79
80 if (hy < 0x0020000000000000LL)
81 {
82 if (x + x > y)
83 {
84 x -= y;
85 ++cquo;
86 if (x + x >= y)
87 {
88 x -= y;
89 ++cquo;
90 }
91 }
92 }
93 else
94 {
95 long double y_half = 0.5L * y;
96 if (x > y_half)
97 {
98 x -= y;
99 ++cquo;
100 if (x >= y_half)
101 {
102 x -= y;
103 ++cquo;
104 }
105 }
106 }
107
108 *quo = qs ? -cquo : cquo;
109
110 /* Ensure correct sign of zero result in round-downward mode. */
111 if (x == 0.0L)
112 x = 0.0L;
113 if (sx)
114 x = -x;
115 return x;
116}
117long_double_symbol (libm, __remquol, remquol);
118

source code of glibc/sysdeps/ieee754/ldbl-128ibm/s_remquol.c