]> git.ipfire.org Git - thirdparty/glibc.git/blame - sysdeps/ieee754/ldbl-128ibm/s_remquol.c
Update copyright dates with scripts/update-copyrights.
[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-128ibm / s_remquol.c
CommitLineData
f964490f 1/* Compute remainder and a congruent to the quotient.
b168057a 2 Copyright (C) 1997-2015 Free Software Foundation, Inc.
f964490f
RM
3 This file is part of the GNU C Library.
4 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
5 Jakub Jelinek <jj@ultra.linux.cz>, 1999.
6
7 The GNU C Library is free software; you can redistribute it and/or
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.
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
15 Lesser General Public License for more details.
16
17 You should have received a copy of the GNU Lesser General Public
59ba27a6
PE
18 License along with the GNU C Library; if not, see
19 <http://www.gnu.org/licenses/>. */
f964490f
RM
20
21#include <math.h>
22
1ed0291c 23#include <math_private.h>
f964490f
RM
24#include <math_ldbl_opt.h>
25
26
27static const long double zero = 0.0;
28
29
30long double
31__remquol (long double x, long double y, int *quo)
32{
33 int64_t hx,hy;
34 u_int64_t sx,lx,ly,qs;
35 int cquo;
765714ca
AM
36 double xhi, xlo, yhi, ylo;
37
38 ldbl_unpack (x, &xhi, &xlo);
39 EXTRACT_WORDS64 (hx, xhi);
40 EXTRACT_WORDS64 (lx, xlo);
41 ldbl_unpack (y, &yhi, &ylo);
42 EXTRACT_WORDS64 (hy, yhi);
43 EXTRACT_WORDS64 (ly, ylo);
f964490f
RM
44 sx = hx & 0x8000000000000000ULL;
45 qs = sx ^ (hy & 0x8000000000000000ULL);
46 hy &= 0x7fffffffffffffffLL;
47 hx &= 0x7fffffffffffffffLL;
48
49 /* Purge off exception values. */
765714ca 50 if (hy == 0)
f964490f
RM
51 return (x * y) / (x * y); /* y = 0 */
52 if ((hx >= 0x7ff0000000000000LL) /* x not finite */
765714ca 53 || (hy > 0x7ff0000000000000LL)) /* y is NaN */
f964490f
RM
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);
66 y = fabsl (y);
67 cquo = 0;
68
69 if (x >= 4 * y)
70 {
71 x -= 4 * y;
72 cquo += 4;
73 }
74 if (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 if (sx)
111 x = -x;
112 return x;
113}
114long_double_symbol (libm, __remquol, remquol);