]> git.ipfire.org Git - thirdparty/glibc.git/blame - sysdeps/ieee754/ldbl-128ibm/s_remquol.c
2.5-18.1
[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-128ibm / s_remquol.c
CommitLineData
0ecb606c
JJ
1/* Compute remainder and a congruent to the quotient.
2 Copyright (C) 1997,1999,2002,2004,2006 Free Software Foundation, Inc.
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
18 License along with the GNU C Library; if not, write to the Free
19 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20 02111-1307 USA. */
21
22#include <math.h>
23
24#include "math_private.h"
25#include <math_ldbl_opt.h>
26
27
28static const long double zero = 0.0;
29
30
31long double
32__remquol (long double x, long double y, int *quo)
33{
34 int64_t hx,hy;
35 u_int64_t sx,lx,ly,qs;
36 int cquo;
37
38 GET_LDOUBLE_WORDS64 (hx, lx, x);
39 GET_LDOUBLE_WORDS64 (hy, ly, y);
40 sx = hx & 0x8000000000000000ULL;
41 qs = sx ^ (hy & 0x8000000000000000ULL);
42 hy &= 0x7fffffffffffffffLL;
43 hx &= 0x7fffffffffffffffLL;
44
45 /* Purge off exception values. */
46 if ((hy | (ly & 0x7fffffffffffffff)) == 0)
47 return (x * y) / (x * y); /* y = 0 */
48 if ((hx >= 0x7ff0000000000000LL) /* x not finite */
49 || ((hy >= 0x7ff0000000000000LL) /* y is NaN */
50 && (((hy - 0x7ff0000000000000LL) | ly) != 0)))
51 return (x * y) / (x * y);
52
53 if (hy <= 0x7fbfffffffffffffLL)
54 x = __ieee754_fmodl (x, 8 * y); /* now x < 8y */
55
56 if (((hx - hy) | (lx - ly)) == 0)
57 {
58 *quo = qs ? -1 : 1;
59 return zero * x;
60 }
61
62 x = fabsl (x);
63 y = fabsl (y);
64 cquo = 0;
65
66 if (x >= 4 * y)
67 {
68 x -= 4 * y;
69 cquo += 4;
70 }
71 if (x >= 2 * y)
72 {
73 x -= 2 * y;
74 cquo += 2;
75 }
76
77 if (hy < 0x0020000000000000LL)
78 {
79 if (x + x > y)
80 {
81 x -= y;
82 ++cquo;
83 if (x + x >= y)
84 {
85 x -= y;
86 ++cquo;
87 }
88 }
89 }
90 else
91 {
92 long double y_half = 0.5L * y;
93 if (x > y_half)
94 {
95 x -= y;
96 ++cquo;
97 if (x >= y_half)
98 {
99 x -= y;
100 ++cquo;
101 }
102 }
103 }
104
105 *quo = qs ? -cquo : cquo;
106
107 if (sx)
108 x = -x;
109 return x;
110}
111long_double_symbol (libm, __remquol, remquol);