xref: /illumos-gate/usr/src/lib/libm/common/C/fmod.c (revision b41e2fb6b2658a36accc73b3b86d99375e700dba)
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
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