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/*
23 * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
24 */
25/*
26 * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
27 * Use is subject to license terms.
28 */
29
30#pragma weak __fmodl = fmodl
31
32#include "libm.h"
33
34static const int
35	is = -0x7fffffff - 1,
36	im = 0x0000ffff,
37	iu = 0x00010000;
38
39static const long double
40	zero = 0.0L,
41	one = 1.0L;
42
43#ifdef __LITTLE_ENDIAN
44#define	__H0(x) *(3 + (int *) &x)
45#define	__H1(x) *(2 + (int *) &x)
46#define	__H2(x) *(1 + (int *) &x)
47#define	__H3(x) *(0 + (int *) &x)
48#else
49#define	__H0(x) *(0 + (int *) &x)
50#define	__H1(x) *(1 + (int *) &x)
51#define	__H2(x) *(2 + (int *) &x)
52#define	__H3(x) *(3 + (int *) &x)
53#endif
54
55long double
56fmodl(long double x, long double y) {
57	long double a, b;
58	int n, ix, iy, k, sx;
59	int hx;
60	int x0, y0, z0, carry;
61	unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
62
63	hx = __H0(x);
64	x1 = __H1(x);
65	x2 = __H2(x);
66	x3 = __H3(x);
67	y0 = __H0(y);
68	y1 = __H1(y);
69	y2 = __H2(y);
70	y3 = __H3(y);
71
72	sx = hx & 0x80000000;
73	x0 = hx ^ sx;
74	y0 &= 0x7fffffff;
75
76	/* purge off exception values */
77	if (x0 >= 0x7fff0000 ||	/* !finitel(x) */
78	    (y0 > 0x7fff0000) || (y0 == 0x7fff0000 && ((y1 | y2 | y3) != 0)) ||
79	    (y0 | y1 | y2 | y3) == 0)	/* isnanl(y) || y = 0 */
80		return ((x * y) / (x * y));
81	a = fabsl(x);
82	b = fabsl(y);
83	if (a <= b) {
84		if (a < b)
85			return (x);
86		else
87			return (zero * x);
88	}
89	/* determine ix = ilogbl(x) */
90	if (x0 < iu) {		/* subnormal x */
91		ix = -16382;
92		while (x0 == 0) {
93			ix -= 16;
94			x0 = x1 >> 16;
95			x1 = (x1 << 16) | (x2 >> 16);
96			x2 = (x2 << 16) | (x3 >> 16);
97			x3 = (x3 << 16);
98		}
99		while (x0 < iu) {
100			ix -= 1;
101			x0 = (x0 << 1) | (x1 >> 31);
102			x1 = (x1 << 1) | (x2 >> 31);
103			x2 = (x2 << 1) | (x3 >> 31);
104			x3 <<= 1;
105		}
106	} else {
107		ix = (x0 >> 16) - 16383;
108		x0 = iu | (x0 & im);
109	}
110
111	/* determine iy = ilogbl(y) */
112	if (y0 < iu) {		/* subnormal y */
113		iy = -16382;
114		while (y0 == 0) {
115			iy -= 16;
116			y0 = y1 >> 16;
117			y1 = (y1 << 16) | (y2 >> 16);
118			y2 = (y2 << 16) | (y3 >> 16);
119			y3 = (y3 << 16);
120		}
121		while (y0 < iu) {
122			iy -= 1;
123			y0 = (y0 << 1) | (y1 >> 31);
124			y1 = (y1 << 1) | (y2 >> 31);
125			y2 = (y2 << 1) | (y3 >> 31);
126			y3 <<= 1;
127		}
128	} else {
129		iy = (y0 >> 16) - 16383;
130		y0 = iu | (y0 & im);
131	}
132
133	/* fix point fmod */
134	n = ix - iy;
135	while (n--) {
136		while (x0 == 0 && n >= 16) {
137			n -= 16;
138			x0 = x1 >> 16;
139			x1 = (x1 << 16) | (x2 >> 16);
140			x2 = (x2 << 16) | (x3 >> 16);
141			x3 = (x3 << 16);
142		}
143		while (x0 < iu && n >= 1) {
144			n -= 1;
145			x0 = (x0 << 1) | (x1 >> 31);
146			x1 = (x1 << 1) | (x2 >> 31);
147			x2 = (x2 << 1) | (x3 >> 31);
148			x3 = (x3 << 1);
149		}
150		carry = 0;
151		z3 = x3 - y3;
152		carry = (z3 > x3);
153		if (carry == 0) {
154			z2 = x2 - y2;
155			carry = (z2 > x2);
156		} else {
157			z2 = x2 - y2 - 1;
158			carry = (z2 >= x2);
159		}
160		if (carry == 0) {
161			z1 = x1 - y1;
162			carry = (z1 > x1);
163		} else {
164			z1 = x1 - y1 - 1;
165			carry = (z1 >= x1);
166		}
167		z0 = x0 - y0 - carry;
168		if (z0 < 0) {	/* double x */
169			x0 = x0 + x0 + ((x1 & is) != 0);
170			x1 = x1 + x1 + ((x2 & is) != 0);
171			x2 = x2 + x2 + ((x3 & is) != 0);
172			x3 = x3 + x3;
173		} else {
174			if (z0 == 0) {
175				if ((z1 | z2 | z3) == 0) {	/* 0: done */
176					__H0(a) = hx & is;
177					__H1(a) = __H2(a) = __H3(a) = 0;
178					return (a);
179				}
180			}
181			/* x = z << 1 */
182			z0 = z0 + z0 + ((z1 & is) != 0);
183			z1 = z1 + z1 + ((z2 & is) != 0);
184			z2 = z2 + z2 + ((z3 & is) != 0);
185			z3 = z3 + z3;
186			x0 = z0;
187			x1 = z1;
188			x2 = z2;
189			x3 = z3;
190		}
191	}
192
193	carry = 0;
194	z3 = x3 - y3;
195	carry = (z3 > x3);
196	if (carry == 0) {
197		z2 = x2 - y2;
198		carry = (z2 > x2);
199	} else {
200		z2 = x2 - y2 - 1;
201		carry = (z2 >= x2);
202	}
203	if (carry == 0) {
204		z1 = x1 - y1;
205		carry = (z1 > x1);
206	} else {
207		z1 = x1 - y1 - 1;
208		carry = (z1 >= x1);
209	}
210	z0 = x0 - y0 - carry;
211	if (z0 >= 0) {
212		x0 = z0;
213		x1 = z1;
214		x2 = z2;
215		x3 = z3;
216	}
217	/* convert back to floating value and restore the sign */
218	if ((x0 | x1 | x2 | x3) == 0) {
219		__H0(a) = hx & is;
220		__H1(a) = __H2(a) = __H3(a) = 0;
221		return (a);
222	}
223	while (x0 < iu) {
224		if (x0 == 0) {
225			iy -= 16;
226			x0 = x1 >> 16;
227			x1 = (x1 << 16) | (x2 >> 16);
228			x2 = (x2 << 16) | (x3 >> 16);
229			x3 = (x3 << 16);
230		} else {
231			x0 = x0 + x0 + ((x1 & is) != 0);
232			x1 = x1 + x1 + ((x2 & is) != 0);
233			x2 = x2 + x2 + ((x3 & is) != 0);
234			x3 = x3 + x3;
235			iy -= 1;
236		}
237	}
238
239	/* normalize output */
240	if (iy >= -16382) {
241		__H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
242		__H1(a) = x1;
243		__H2(a) = x2;
244		__H3(a) = x3;
245	} else {		/* subnormal output */
246		n = -16382 - iy;
247		k = n & 31;
248		if (k != 0) {
249			if (k <= 16) {
250				x3 = (x2 << (32 - k)) | (x3 >> k);
251				x2 = (x1 << (32 - k)) | (x2 >> k);
252				x1 = (x0 << (32 - k)) | (x1 >> k);
253				x0 >>= k;
254			} else {
255				x3 = (x2 << (32 - k)) | (x3 >> k);
256				x2 = (x1 << (32 - k)) | (x2 >> k);
257				x1 = (x0 << (32 - k)) | (x1 >> k);
258				x0 = 0;
259			}
260		}
261		while (n >= 32) {
262			n -= 32;
263			x3 = x2;
264			x2 = x1;
265			x1 = x0;
266			x0 = 0;
267		}
268		__H0(a) = x0 | sx;
269		__H1(a) = x1;
270		__H2(a) = x2;
271		__H3(a) = x3;
272		a *= one;
273	}
274	return (a);
275}
276