1 | /* |
2 | * ==================================================== |
3 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
4 | * |
5 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
6 | * Permission to use, copy, modify, and distribute this |
7 | * software is freely granted, provided that this notice |
8 | * is preserved. |
9 | * ==================================================== |
10 | */ |
11 | |
12 | /* |
13 | * __ieee754_fmod(x,y) |
14 | * Return x mod y in exact arithmetic |
15 | * Method: shift and subtract |
16 | */ |
17 | |
18 | #include <math.h> |
19 | #include <math_private.h> |
20 | #include <stdint.h> |
21 | #include <libm-alias-finite.h> |
22 | |
23 | static const double one = 1.0, Zero[] = {0.0, -0.0,}; |
24 | |
25 | double |
26 | __ieee754_fmod (double x, double y) |
27 | { |
28 | int32_t n,ix,iy; |
29 | int64_t hx,hy,hz,sx,i; |
30 | |
31 | EXTRACT_WORDS64(hx,x); |
32 | EXTRACT_WORDS64(hy,y); |
33 | sx = hx&UINT64_C(0x8000000000000000); /* sign of x */ |
34 | hx ^=sx; /* |x| */ |
35 | hy &= UINT64_C(0x7fffffffffffffff); /* |y| */ |
36 | |
37 | /* purge off exception values */ |
38 | if(__builtin_expect(hy==0 |
39 | || hx >= UINT64_C(0x7ff0000000000000) |
40 | || hy > UINT64_C(0x7ff0000000000000), 0)) |
41 | /* y=0,or x not finite or y is NaN */ |
42 | return (x*y)/(x*y); |
43 | if(__builtin_expect(hx<=hy, 0)) { |
44 | if(hx<hy) return x; /* |x|<|y| return x */ |
45 | return Zero[(uint64_t)sx>>63]; /* |x|=|y| return x*0*/ |
46 | } |
47 | |
48 | /* determine ix = ilogb(x) */ |
49 | if(__builtin_expect(hx<UINT64_C(0x0010000000000000), 0)) { |
50 | /* subnormal x */ |
51 | for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; |
52 | } else ix = (hx>>52)-1023; |
53 | |
54 | /* determine iy = ilogb(y) */ |
55 | if(__builtin_expect(hy<UINT64_C(0x0010000000000000), 0)) { /* subnormal y */ |
56 | for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; |
57 | } else iy = (hy>>52)-1023; |
58 | |
59 | /* set up hx, hy and align y to x */ |
60 | if(__builtin_expect(ix >= -1022, 1)) |
61 | hx = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hx); |
62 | else { /* subnormal x, shift x to normal */ |
63 | n = -1022-ix; |
64 | hx<<=n; |
65 | } |
66 | if(__builtin_expect(iy >= -1022, 1)) |
67 | hy = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hy); |
68 | else { /* subnormal y, shift y to normal */ |
69 | n = -1022-iy; |
70 | hy<<=n; |
71 | } |
72 | |
73 | /* fix point fmod */ |
74 | n = ix - iy; |
75 | while(n--) { |
76 | hz=hx-hy; |
77 | if(hz<0){hx = hx+hx;} |
78 | else { |
79 | if(hz==0) /* return sign(x)*0 */ |
80 | return Zero[(uint64_t)sx>>63]; |
81 | hx = hz+hz; |
82 | } |
83 | } |
84 | hz=hx-hy; |
85 | if(hz>=0) {hx=hz;} |
86 | |
87 | /* convert back to floating value and restore the sign */ |
88 | if(hx==0) /* return sign(x)*0 */ |
89 | return Zero[(uint64_t)sx>>63]; |
90 | while(hx<UINT64_C(0x0010000000000000)) { /* normalize x */ |
91 | hx = hx+hx; |
92 | iy -= 1; |
93 | } |
94 | if(__builtin_expect(iy>= -1022, 1)) { /* normalize output */ |
95 | hx = ((hx-UINT64_C(0x0010000000000000))|((uint64_t)(iy+1023)<<52)); |
96 | INSERT_WORDS64(x,hx|sx); |
97 | } else { /* subnormal output */ |
98 | n = -1022 - iy; |
99 | hx>>=n; |
100 | INSERT_WORDS64(x,hx|sx); |
101 | x *= one; /* create necessary signal */ |
102 | } |
103 | return x; /* exact output */ |
104 | } |
105 | libm_alias_finite (__ieee754_fmod, __fmod) |
106 | |