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