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, Version 1.0 only
6 * (the "License").  You may not use this file except in compliance
7 * with the License.
8 *
9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
10 * or http://www.opensolaris.org/os/licensing.
11 * See the License for the specific language governing permissions
12 * and limitations under the License.
13 *
14 * When distributing Covered Code, include this CDDL HEADER in each
15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
16 * If applicable, add the following below this CDDL HEADER, with the
17 * fields enclosed by brackets "[]" replaced with your own identifying
18 * information: Portions Copyright [yyyy] [name of copyright owner]
19 *
20 * CDDL HEADER END
21 */
22/*
23 * Copyright 1988 Sun Microsystems, Inc.  All rights reserved.
24 * Use is subject to license terms.
25 */
26
27#pragma ident	"%Z%%M%	%I%	%E% SMI"
28
29/* Special version adapted from libm for use in libc. */
30
31static int	n0 = 0, n1 = 1;
32
33static double   two52 = 4.503599627370496000E+15;
34static double   twom52 = 2.220446049250313081E-16;
35
36static double
37setexception(int n, double x)
38{
39	return (0.0);
40}
41
42double
43copysign(double x, double y)
44{
45	long           *px = (long *) &x;
46	long           *py = (long *) &y;
47	px[n0] = (px[n0] & 0x7fffffff) | (py[n0] & 0x80000000);
48	return (x);
49}
50
51static double
52fabs(double x)
53{
54	long           *px = (long *) &x;
55	px[0] &= 0x7fffffff;
56
57	return (x);
58}
59
60static int
61finite(double x)
62{
63	long           *px = (long *) &x;
64	return ((px[n0] & 0x7ff00000) != 0x7ff00000);
65}
66
67static int
68ilogb(double x)
69{
70	long           *px = (long *) &x, k;
71	k = px[n0] & 0x7ff00000;
72	if (k == 0) {
73		if ((px[n1] | (px[n0] & 0x7fffffff)) == 0)
74			return (0x80000001);
75		else {
76			x *= two52;
77			return ((px[n0] & 0x7ff00000) >> 20) - 1075;
78		}
79	} else if (k != 0x7ff00000)
80		return (k >> 20) - 1023;
81	else
82		return (0x7fffffff);
83}
84
85static double
86scalbn(double x, int n)
87{
88	long           *px = (long *) &x, k;
89	double          twom54 = twom52 * 0.25;
90	k = (px[n0] & 0x7ff00000) >> 20;
91	if (k == 0x7ff)
92		return (x + x);
93	if ((px[n1] | (px[n0] & 0x7fffffff)) == 0)
94		return (x);
95	if (k == 0) {
96		x *= two52;
97		k = ((px[n0] & 0x7ff00000) >> 20) - 52;
98	}
99	k = k + n;
100	if (n > 5000)
101		return (setexception(2, x));
102	if (n < -5000)
103		return (setexception(1, x));
104	if (k > 0x7fe)
105		return (setexception(2, x));
106	if (k <= -54)
107		return (setexception(1, x));
108	if (k > 0) {
109		px[n0] = (px[n0] & 0x800fffff) | (k << 20);
110		return (x);
111	}
112	k += 54;
113	px[n0] = (px[n0] & 0x800fffff) | (k << 20);
114	return (x * twom54);
115}
116
117double
118fmod(double x, double y)
119{
120	int             ny, nr;
121	double          r, z, w;
122
123	int finite(), ilogb();
124	double fabs(), scalbn(), copysign();
125
126	/* purge off exception values */
127	if (!finite(x) || y != y || y == 0.0) {
128		return ((x * y) / (x * y));
129	}
130	/* scale and subtract to get the remainder */
131	r = fabs(x);
132	y = fabs(y);
133	ny = ilogb(y);
134	while (r >= y) {
135		nr = ilogb(r);
136		if (nr == ny)
137			w = y;
138		else {
139			z = scalbn(y, nr - ny - 1);
140			w = z + z;
141		}
142		if (r >= w)
143			r -= w;
144		else
145			r -= z;
146	}
147
148	/* restore sign */
149	return (copysign(r, x));
150}
151