]>
Commit | Line | Data |
---|---|---|
bc9f6000 | 1 | /* Compute remainder and a congruent to the quotient. |
2b778ceb | 2 | Copyright (C) 1997-2021 Free Software Foundation, Inc. |
bc9f6000 | 3 | This file is part of the GNU C Library. |
abfbdde1 UD |
4 | Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and |
5 | Jakub Jelinek <jj@ultra.linux.cz>, 1999. | |
bc9f6000 UD |
6 | |
7 | The GNU C Library is free software; you can redistribute it and/or | |
41bdb6e2 AJ |
8 | modify it under the terms of the GNU Lesser General Public |
9 | License as published by the Free Software Foundation; either | |
10 | version 2.1 of the License, or (at your option) any later version. | |
bc9f6000 UD |
11 | |
12 | The GNU C Library is distributed in the hope that it will be useful, | |
13 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
41bdb6e2 | 15 | Lesser General Public License for more details. |
bc9f6000 | 16 | |
41bdb6e2 | 17 | You should have received a copy of the GNU Lesser General Public |
59ba27a6 | 18 | License along with the GNU C Library; if not, see |
5a82c748 | 19 | <https://www.gnu.org/licenses/>. */ |
bc9f6000 UD |
20 | |
21 | #include <math.h> | |
22 | ||
1ed0291c | 23 | #include <math_private.h> |
fd3b4e7c | 24 | #include <libm-alias-ldouble.h> |
bc9f6000 UD |
25 | |
26 | ||
15089e04 | 27 | static const _Float128 zero = 0.0; |
bc9f6000 UD |
28 | |
29 | ||
15089e04 PM |
30 | _Float128 |
31 | __remquol (_Float128 x, _Float128 y, int *quo) | |
bc9f6000 | 32 | { |
abfbdde1 | 33 | int64_t hx,hy; |
24ab7723 | 34 | uint64_t sx,lx,ly,qs; |
52e1b618 | 35 | int cquo; |
bc9f6000 | 36 | |
abfbdde1 UD |
37 | GET_LDOUBLE_WORDS64 (hx, lx, x); |
38 | GET_LDOUBLE_WORDS64 (hy, ly, y); | |
39 | sx = hx & 0x8000000000000000ULL; | |
40 | qs = sx ^ (hy & 0x8000000000000000ULL); | |
41 | hy &= 0x7fffffffffffffffLL; | |
42 | hx &= 0x7fffffffffffffffLL; | |
bc9f6000 UD |
43 | |
44 | /* Purge off exception values. */ | |
34b402e5 UD |
45 | if ((hy | ly) == 0) |
46 | return (x * y) / (x * y); /* y = 0 */ | |
abfbdde1 UD |
47 | if ((hx >= 0x7fff000000000000LL) /* x not finite */ |
48 | || ((hy >= 0x7fff000000000000LL) /* y is NaN */ | |
49 | && (((hy - 0x7fff000000000000LL) | ly) != 0))) | |
34b402e5 | 50 | return (x * y) / (x * y); |
bc9f6000 | 51 | |
abfbdde1 UD |
52 | if (hy <= 0x7ffbffffffffffffLL) |
53 | x = __ieee754_fmodl (x, 8 * y); /* now x < 8y */ | |
9c84384c | 54 | |
34b402e5 | 55 | if (((hx - hy) | (lx - ly)) == 0) |
bc9f6000 UD |
56 | { |
57 | *quo = qs ? -1 : 1; | |
58 | return zero * x; | |
59 | } | |
60 | ||
abfbdde1 UD |
61 | x = fabsl (x); |
62 | y = fabsl (y); | |
bc9f6000 UD |
63 | cquo = 0; |
64 | ||
a820f9b3 | 65 | if (hy <= 0x7ffcffffffffffffLL && x >= 4 * y) |
bc9f6000 | 66 | { |
34b402e5 | 67 | x -= 4 * y; |
63551311 UD |
68 | cquo += 4; |
69 | } | |
a820f9b3 | 70 | if (hy <= 0x7ffdffffffffffffLL && x >= 2 * y) |
63551311 UD |
71 | { |
72 | x -= 2 * y; | |
bc9f6000 UD |
73 | cquo += 2; |
74 | } | |
bc9f6000 | 75 | |
abfbdde1 | 76 | if (hy < 0x0002000000000000LL) |
bc9f6000 | 77 | { |
34b402e5 | 78 | if (x + x > y) |
bc9f6000 | 79 | { |
34b402e5 | 80 | x -= y; |
d705269e | 81 | ++cquo; |
34b402e5 | 82 | if (x + x >= y) |
63551311 UD |
83 | { |
84 | x -= y; | |
85 | ++cquo; | |
86 | } | |
bc9f6000 UD |
87 | } |
88 | } | |
89 | else | |
90 | { | |
02bbfb41 | 91 | _Float128 y_half = L(0.5) * y; |
63551311 | 92 | if (x > y_half) |
bc9f6000 | 93 | { |
34b402e5 | 94 | x -= y; |
d705269e | 95 | ++cquo; |
34b402e5 | 96 | if (x >= y_half) |
63551311 UD |
97 | { |
98 | x -= y; | |
99 | ++cquo; | |
100 | } | |
bc9f6000 UD |
101 | } |
102 | } | |
103 | ||
104 | *quo = qs ? -cquo : cquo; | |
105 | ||
ce8fc784 | 106 | /* Ensure correct sign of zero result in round-downward mode. */ |
02bbfb41 PM |
107 | if (x == 0) |
108 | x = 0; | |
63551311 UD |
109 | if (sx) |
110 | x = -x; | |
bc9f6000 UD |
111 | return x; |
112 | } | |
fd3b4e7c | 113 | libm_alias_ldouble (__remquo, remquo) |