1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21 /*
22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 */
24 /*
25 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
26 * Use is subject to license terms.
27 */
28
29 #pragma weak __fmod = fmod
30
31 #include "libm.h"
32
33 static const double zero = 0.0;
34
35 /*
36 * The following implementation assumes fast 64-bit integer arith-
37 * metic. This is fine for sparc because we build libm in v8plus
38 * mode. It's also fine for sparcv9 and amd64, although we have
39 * assembly code on amd64. For x86, it would be better to use
40 * 32-bit code, but we have assembly for x86, too.
41 */
42 double
fmod(double x,double y)43 fmod(double x, double y) {
44 double w;
45 long long hx, ix, iy, iz;
46 int nd, k, ny;
47
48 hx = *(long long *)&x;
49 ix = hx & ~0x8000000000000000ull;
50 iy = *(long long *)&y & ~0x8000000000000000ull;
51
52 /* handle special cases */
53 if (iy == 0ll)
54 return (_SVID_libm_err(x, y, 27));
55
56 if (ix >= 0x7ff0000000000000ll || iy > 0x7ff0000000000000ll)
57 return ((x * y) * zero);
58
59 if (ix <= iy)
60 return ((ix < iy)? x : x * zero);
61
62 /*
63 * Set:
64 * ny = true exponent of y
65 * nd = true exponent of x minus true exponent of y
66 * ix = normalized significand of x
67 * iy = normalized significand of y
68 */
69 ny = iy >> 52;
70 k = ix >> 52;
71 if (ny == 0) {
72 /* y is subnormal, x could be normal or subnormal */
73 ny = 1;
74 while (iy < 0x0010000000000000ll) {
75 ny -= 1;
76 iy += iy;
77 }
78 nd = k - ny;
79 if (k == 0) {
80 nd += 1;
81 while (ix < 0x0010000000000000ll) {
82 nd -= 1;
83 ix += ix;
84 }
85 } else {
86 ix = 0x0010000000000000ll | (ix & 0x000fffffffffffffll);
87 }
88 } else {
89 /* both x and y are normal */
90 nd = k - ny;
91 ix = 0x0010000000000000ll | (ix & 0x000fffffffffffffll);
92 iy = 0x0010000000000000ll | (iy & 0x000fffffffffffffll);
93 }
94
95 /* perform fixed point mod */
96 while (nd--) {
97 iz = ix - iy;
98 if (iz >= 0)
99 ix = iz;
100 ix += ix;
101 }
102 iz = ix - iy;
103 if (iz >= 0)
104 ix = iz;
105
106 /* convert back to floating point and restore the sign */
107 if (ix == 0ll)
108 return (x * zero);
109 while (ix < 0x0010000000000000ll) {
110 ix += ix;
111 ny -= 1;
112 }
113 while (ix > 0x0020000000000000ll) { /* XXX can this ever happen? */
114 ny += 1;
115 ix >>= 1;
116 }
117 if (ny <= 0) {
118 /* result is subnormal */
119 k = -ny + 1;
120 ix >>= k;
121 *(long long *)&w = (hx & 0x8000000000000000ull) | ix;
122 return (w);
123 }
124 *(long long *)&w = (hx & 0x8000000000000000ull) |
125 ((long long)ny << 52) | (ix & 0x000fffffffffffffll);
126 return (w);
127 }
128