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 __remquol = remquol
31
32#include "libm.h"
33#if defined(__SUNPRO_C)
34#include <sunmath.h>			/* fabsl */
35#endif
36/* INDENT OFF */
37static const int
38	is = -0x7fffffff - 1,
39	im = 0x0000ffff,
40	iu = 0x00010000;
41
42static const long double zero = 0.0L, one = 1.0L;
43/* INDENT ON */
44
45#if defined(__sparc)
46#define	__H0(x)	((int *) &x)[0]
47#define	__H1(x)	((int *) &x)[1]
48#define	__H2(x)	((int *) &x)[2]
49#define	__H3(x)	((int *) &x)[3]
50#else
51#error Unsupported architecture
52#endif
53
54/*
55 * On entrance: *quo is initialized to 0, x finite and y non-zero & ordered
56 */
57static long double
58fmodquol(long double x, long double y, int *quo) {
59	long double a, b;
60	int n, ix, iy, k, sx, sq, m;
61	int hx;
62	int x0, y0, z0, carry;
63	unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
64
65	hx = __H0(x);
66	x1 = __H1(x);
67	x2 = __H2(x);
68	x3 = __H3(x);
69	y0 = __H0(y);
70	y1 = __H1(y);
71	y2 = __H2(y);
72	y3 = __H3(y);
73
74	sx = hx & is;
75	sq = (hx ^ y0) & is;
76	x0 = hx ^ sx;
77	y0 &= ~0x80000000;
78
79	a = fabsl(x);
80	b = fabsl(y);
81	if (a <= b) {
82		if (a < b)
83			return (x);
84		else {
85			*quo = 1 + (sq >> 30);
86			return (zero * x);
87		}
88	}
89	/* determine ix = ilogbl(x) */
90	if (x0 < iu) {		/* subnormal x */
91		ix = 0;
92		ix = -16382;
93		while (x0 == 0) {
94			ix -= 16;
95			x0 = x1 >> 16;
96			x1 = (x1 << 16) | (x2 >> 16);
97			x2 = (x2 << 16) | (x3 >> 16);
98			x3 = (x3 << 16);
99		}
100		while (x0 < iu) {
101			ix -= 1;
102			x0 = (x0 << 1) | (x1 >> 31);
103			x1 = (x1 << 1) | (x2 >> 31);
104			x2 = (x2 << 1) | (x3 >> 31);
105			x3 <<= 1;
106		}
107	} else {
108		ix = (x0 >> 16) - 16383;
109		x0 = iu | (x0 & im);
110	}
111
112	/* determine iy = ilogbl(y) */
113	if (y0 < iu) {		/* subnormal y */
114		iy = -16382;
115		while (y0 == 0) {
116			iy -= 16;
117			y0 = y1 >> 16;
118			y1 = (y1 << 16) | (y2 >> 16);
119			y2 = (y2 << 16) | (y3 >> 16);
120			y3 = (y3 << 16);
121		}
122		while (y0 < iu) {
123			iy -= 1;
124			y0 = (y0 << 1) | (y1 >> 31);
125			y1 = (y1 << 1) | (y2 >> 31);
126			y2 = (y2 << 1) | (y3 >> 31);
127			y3 <<= 1;
128		}
129	} else {
130		iy = (y0 >> 16) - 16383;
131		y0 = iu | (y0 & im);
132	}
133
134
135	/* fix point fmod */
136	n = ix - iy;
137	m = 0;
138	while (n--) {
139		while (x0 == 0 && n >= 16) {
140			m <<= 16;
141			n -= 16;
142			x0 = x1 >> 16;
143			x1 = (x1 << 16) | (x2 >> 16);
144			x2 = (x2 << 16) | (x3 >> 16);
145			x3 = (x3 << 16);
146		}
147		while (x0 < iu && n >= 1) {
148			m += m;
149			n -= 1;
150			x0 = (x0 << 1) | (x1 >> 31);
151			x1 = (x1 << 1) | (x2 >> 31);
152			x2 = (x2 << 1) | (x3 >> 31);
153			x3 = (x3 << 1);
154		}
155		carry = 0;
156		z3 = x3 - y3;
157		carry = z3 > x3;
158		if (carry == 0) {
159			z2 = x2 - y2;
160			carry = z2 > x2;
161		} else {
162			z2 = x2 - y2 - 1;
163			carry = z2 >= x2;
164		}
165		if (carry == 0) {
166			z1 = x1 - y1;
167			carry = z1 > x1;
168		} else {
169			z1 = x1 - y1 - 1;
170			carry = z1 >= x1;
171		}
172		z0 = x0 - y0 - carry;
173		if (z0 < 0) {	/* double x */
174			x0 = x0 + x0 + ((x1 & is) != 0);
175			x1 = x1 + x1 + ((x2 & is) != 0);
176			x2 = x2 + x2 + ((x3 & is) != 0);
177			x3 = x3 + x3;
178			m += m;
179		} else {
180			m += 1;
181			if (z0 == 0) {
182				if ((z1 | z2 | z3) == 0) {
183					/* 0: we are done */
184					if (n < 31)
185						m <<= (1 + n);
186					else
187						m = 0;
188					m &= ~0x80000000;
189					*quo = sq >= 0 ? m : -m;
190					__H0(a) = hx & is;
191					__H1(a) = __H2(a) = __H3(a) = 0;
192					return (a);
193				}
194			}
195			/* x = z << 1 */
196			z0 = z0 + z0 + ((z1 & is) != 0);
197			z1 = z1 + z1 + ((z2 & is) != 0);
198			z2 = z2 + z2 + ((z3 & is) != 0);
199			z3 = z3 + z3;
200			x0 = z0;
201			x1 = z1;
202			x2 = z2;
203			x3 = z3;
204			m += m;
205		}
206	}
207	carry = 0;
208	z3 = x3 - y3;
209	carry = z3 > x3;
210	if (carry == 0) {
211		z2 = x2 - y2;
212		carry = z2 > x2;
213	} else {
214		z2 = x2 - y2 - 1;
215		carry = z2 >= x2;
216	}
217	if (carry == 0) {
218		z1 = x1 - y1;
219		carry = z1 > x1;
220	} else {
221		z1 = x1 - y1 - 1;
222		carry = z1 >= x1;
223	}
224	z0 = x0 - y0 - carry;
225	if (z0 >= 0) {
226		x0 = z0;
227		x1 = z1;
228		x2 = z2;
229		x3 = z3;
230		m += 1;
231	}
232	m &= ~0x80000000;
233	*quo = sq >= 0 ? m : -m;
234
235	/* convert back to floating value and restore the sign */
236	if ((x0 | x1 | x2 | x3) == 0) {
237		__H0(a) = hx & is;
238		__H1(a) = __H2(a) = __H3(a) = 0;
239		return (a);
240	}
241	while (x0 < iu) {
242		if (x0 == 0) {
243			iy -= 16;
244			x0 = x1 >> 16;
245			x1 = (x1 << 16) | (x2 >> 16);
246			x2 = (x2 << 16) | (x3 >> 16);
247			x3 = (x3 << 16);
248		} else {
249			x0 = x0 + x0 + ((x1 & is) != 0);
250			x1 = x1 + x1 + ((x2 & is) != 0);
251			x2 = x2 + x2 + ((x3 & is) != 0);
252			x3 = x3 + x3;
253			iy -= 1;
254		}
255	}
256
257	/* normalize output */
258	if (iy >= -16382) {
259		__H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
260		__H1(a) = x1;
261		__H2(a) = x2;
262		__H3(a) = x3;
263	} else {		/* subnormal output */
264		n = -16382 - iy;
265		k = n & 31;
266		if (k <= 16) {
267			x3 = (x2 << (32 - k)) | (x3 >> k);
268			x2 = (x1 << (32 - k)) | (x2 >> k);
269			x1 = (x0 << (32 - k)) | (x1 >> k);
270			x0 >>= k;
271		} else {
272			x3 = (x2 << (32 - k)) | (x3 >> k);
273			x2 = (x1 << (32 - k)) | (x2 >> k);
274			x1 = (x0 << (32 - k)) | (x1 >> k);
275			x0 = 0;
276		}
277		while (n >= 32) {
278			n -= 32;
279			x3 = x2;
280			x2 = x1;
281			x1 = x0;
282			x0 = 0;
283		}
284		__H0(a) = x0 | sx;
285		__H1(a) = x1;
286		__H2(a) = x2;
287		__H3(a) = x3;
288		a *= one;
289	}
290	return (a);
291}
292
293long double
294remquol(long double x, long double y, int *quo) {
295	int hx, hy, sx, sq;
296	long double v;
297
298	hx = __H0(x);		/* high word of x */
299	hy = __H0(y);		/* high word of y */
300	sx = hx & is;		/* sign of x */
301	sq = (hx ^ hy) & is;	/* sign of x/y */
302	hx ^= sx;		/* |x| */
303	hy &= ~0x80000000;
304
305	/* purge off exception values */
306	*quo = 0;
307	/* y=0, y is NaN, x is NaN or inf */
308	if (y == 0.0L || y != y || hx >= 0x7fff0000)
309		return ((x * y) / (x * y));
310
311	y = fabsl(y);
312	x = fabsl(x);
313	if (hy <= 0x7ffdffff) {
314		x = fmodquol(x, y + y, quo);
315		*quo = ((*quo) & 0x3fffffff) << 1;
316	}
317	if (hy < 0x00020000) {
318		if (x + x > y) {
319			*quo += 1;
320			if (x == y)
321				x = zero;
322			else
323				x -= y;
324			if (x + x >= y) {
325				x -= y;
326				*quo += 1;
327			}
328		}
329	} else {
330		v = 0.5L * y;
331		if (x > v) {
332			*quo += 1;
333			if (x == y)
334				x = zero;
335			else
336				x -= y;
337			if (x >= v) {
338				x -= y;
339				*quo += 1;
340			}
341		}
342	}
343	if (sq != 0)
344		*quo = -(*quo);
345	return (sx == 0 ? x : -x);
346}
347