1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
/* e_fmodl.c -- long double version of e_fmod.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* __ieee754_fmodl(x,y)
* Return x mod y in exact arithmetic
* Method: shift and subtract
*/
#include <math.h>
#include <math_private.h>
#include <ieee754.h>
static const long double one = 1.0, Zero[] = {0.0, -0.0,};
long double
__ieee754_fmodl (long double x, long double y)
{
int64_t hx, hy, hz, sx, sy;
uint64_t lx, ly, lz;
int n, ix, iy;
double xhi, xlo, yhi, ylo;
ldbl_unpack (x, &xhi, &xlo);
EXTRACT_WORDS64 (hx, xhi);
EXTRACT_WORDS64 (lx, xlo);
ldbl_unpack (y, &yhi, &ylo);
EXTRACT_WORDS64 (hy, yhi);
EXTRACT_WORDS64 (ly, ylo);
sx = hx&0x8000000000000000ULL; /* sign of x */
hx ^= sx; /* |x| */
sy = hy&0x8000000000000000ULL; /* sign of y */
hy ^= sy; /* |y| */
/* purge off exception values */
if(__builtin_expect(hy==0 ||
(hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */
(hy>0x7ff0000000000000LL),0)) /* or y is NaN */
return (x*y)/(x*y);
if (__glibc_unlikely (hx <= hy))
{
/* If |x| < |y| return x. */
if (hx < hy)
return x;
/* At this point the absolute value of the high doubles of
x and y must be equal. */
if ((lx & 0x7fffffffffffffffLL) == 0
&& (ly & 0x7fffffffffffffffLL) == 0)
/* Both low parts are zero. The result should be an
appropriately signed zero, but the subsequent logic
could treat them as unequal, depending on the signs
of the low parts. */
return Zero[(uint64_t) sx >> 63];
/* If the low double of y is the same sign as the high
double of y (ie. the low double increases |y|)... */
if (((ly ^ sy) & 0x8000000000000000LL) == 0
/* ... then a different sign low double to high double
for x or same sign but lower magnitude... */
&& (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy))
/* ... means |x| < |y|. */
return x;
/* If the low double of x differs in sign to the high
double of x (ie. the low double decreases |x|)... */
if (((lx ^ sx) & 0x8000000000000000LL) != 0
/* ... then a different sign low double to high double
for y with lower magnitude (we've already caught
the same sign for y case above)... */
&& (int64_t) (lx ^ sx) > (int64_t) (ly ^ sy))
/* ... means |x| < |y|. */
return x;
/* If |x| == |y| return x*0. */
if ((lx ^ sx) == (ly ^ sy))
return Zero[(uint64_t) sx >> 63];
}
/* Make the IBM extended format 105 bit mantissa look like the ieee854 112
bit mantissa so the following operations will give the correct
result. */
ldbl_extract_mantissa(&hx, &lx, &ix, x);
ldbl_extract_mantissa(&hy, &ly, &iy, y);
if (__glibc_unlikely (ix == -IEEE754_DOUBLE_BIAS))
{
/* subnormal x, shift x to normal. */
while ((hx & (1LL << 48)) == 0)
{
hx = (hx << 1) | (lx >> 63);
lx = lx << 1;
ix -= 1;
}
}
if (__glibc_unlikely (iy == -IEEE754_DOUBLE_BIAS))
{
/* subnormal y, shift y to normal. */
while ((hy & (1LL << 48)) == 0)
{
hy = (hy << 1) | (ly >> 63);
ly = ly << 1;
iy -= 1;
}
}
/* fix point fmod */
n = ix - iy;
while(n--) {
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
else {
if((hz|lz)==0) /* return sign(x)*0 */
return Zero[(uint64_t)sx>>63];
hx = hz+hz+(lz>>63); lx = lz+lz;
}
}
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz>=0) {hx=hz;lx=lz;}
/* convert back to floating value and restore the sign */
if((hx|lx)==0) /* return sign(x)*0 */
return Zero[(uint64_t)sx>>63];
while(hx<0x0001000000000000LL) { /* normalize x */
hx = hx+hx+(lx>>63); lx = lx+lx;
iy -= 1;
}
if(__builtin_expect(iy>= -1022,0)) { /* normalize output */
x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
} else { /* subnormal output */
n = -1022 - iy;
/* We know 1 <= N <= 52, and that there are no nonzero
bits in places below 2^-1074. */
lx = (lx >> n) | ((uint64_t) hx << (64 - n));
hx >>= n;
x = ldbl_insert_mantissa((sx>>63), -1023, hx, lx);
x *= one; /* create necessary signal */
}
return x; /* exact output */
}
strong_alias (__ieee754_fmodl, __fmodl_finite)
|
