--- /dev/null	Sat May 10 12:08:49 2014
+++ new/usr/src/lib/libm/common/C/__rem_pio2.c	Sat May 10 12:08:49 2014
@@ -0,0 +1,168 @@
+/*
+ * CDDL HEADER START
+ *
+ * The contents of this file are subject to the terms of the
+ * Common Development and Distribution License (the "License").
+ * You may not use this file except in compliance with the License.
+ *
+ * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
+ * or http://www.opensolaris.org/os/licensing.
+ * See the License for the specific language governing permissions
+ * and limitations under the License.
+ *
+ * When distributing Covered Code, include this CDDL HEADER in each
+ * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
+ * If applicable, add the following below this CDDL HEADER, with the
+ * fields enclosed by brackets "[]" replaced with your own identifying
+ * information: Portions Copyright [yyyy] [name of copyright owner]
+ *
+ * CDDL HEADER END
+ */
+/*
+ * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
+ */
+/*
+ * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
+ * Use is subject to license terms.
+ */
+
+/*
+ * __rem_pio2(x, y) passes back a better-than-double-precision
+ * approximation to x mod pi/2 in y[0]+y[1] and returns an integer
+ * congruent mod 8 to the integer part of x/(pi/2).
+ *
+ * This implementation tacitly assumes that x is finite and at
+ * least about pi/4 in magnitude.
+ */
+
+#include "libm.h"
+
+extern const int _TBL_ipio2_inf[];
+
+/* INDENT OFF */
+/*
+ * invpio2:  53 bits of 2/pi
+ * pio2_1:   first  33 bit of pi/2
+ * pio2_1t:  pi/2 - pio2_1
+ * pio2_2:   second 33 bit of pi/2
+ * pio2_2t:  pi/2 - pio2_2
+ * pio2_3:   third  33 bit of pi/2
+ * pio2_3t:  pi/2 - pio2_3
+ */
+static const double
+	half	= 0.5,
+	invpio2	= 0.636619772367581343075535,	/* 2^ -1  * 1.45F306DC9C883 */
+	pio2_1	= 1.570796326734125614166,	/* 2^  0  * 1.921FB54400000 */
+	pio2_1t	= 6.077100506506192601475e-11,	/* 2^-34  * 1.0B4611A626331 */
+	pio2_2	= 6.077100506303965976596e-11,	/* 2^-34  * 1.0B4611A600000 */
+	pio2_2t	= 2.022266248795950732400e-21,	/* 2^-69  * 1.3198A2E037073 */
+	pio2_3	= 2.022266248711166455796e-21,	/* 2^-69  * 1.3198A2E000000 */
+	pio2_3t	= 8.478427660368899643959e-32;	/* 2^-104 * 1.B839A252049C1 */
+/* INDENT ON */
+
+int
+__rem_pio2(double x, double *y) {
+	double	w, t, r, fn;
+	double	tx[3];
+	int	e0, i, j, nx, n, ix, hx, lx;
+
+	hx = ((int *)&x)[HIWORD];
+	ix = hx & 0x7fffffff;
+
+	if (ix < 0x4002d97c) {
+		/* |x| < 3pi/4, special case with n=1 */
+		t = fabs(x) - pio2_1;
+		if (ix != 0x3ff921fb) {	/* 33+53 bit pi is good enough */
+			y[0] = t - pio2_1t;
+			y[1] = (t - y[0]) - pio2_1t;
+		} else {		/* near pi/2, use 33+33+53 bit pi */
+			t -= pio2_2;
+			y[0] = t - pio2_2t;
+			y[1] = (t - y[0]) - pio2_2t;
+		}
+		if (hx < 0) {
+			y[0] = -y[0];
+			y[1] = -y[1];
+			return (-1);
+		}
+		return (1);
+	}
+
+	if (ix <= 0x413921fb) {
+		/* |x| <= 2^19 pi */
+		t = fabs(x);
+		n = (int)(t * invpio2 + half);
+		fn = (double)n;
+		r = t - fn * pio2_1;
+		j = ix >> 20;
+		w = fn * pio2_1t;	/* 1st round good to 85 bit */
+		y[0] = r - w;
+		i = j - ((((int *)y)[HIWORD] >> 20) & 0x7ff);
+		if (i > 16) {	/* 2nd iteration (rare) */
+			/* 2nd round good to 118 bit */
+			if (i < 35) {
+				t = r;	/* r-fn*pio2_2 may not be exact */
+				w = fn * pio2_2;
+				r = t - w;
+				w = fn * pio2_2t - ((t - r) - w);
+				y[0] = r - w;
+			} else {
+				r -= fn * pio2_2;
+				w = fn * pio2_2t;
+				y[0] = r - w;
+				i = j - ((((int *)y)[HIWORD] >> 20) & 0x7ff);
+				if (i > 49) {
+					/* 3rd iteration (extremely rare) */
+					if (i < 68) {
+						t = r;
+						w = fn * pio2_3;
+						r = t - w;
+						w = fn * pio2_3t -
+						    ((t - r) - w);
+						y[0] = r - w;
+					} else {
+						/*
+						 * 3rd round good to 151 bits;
+						 * covered all possible cases
+						 */
+						r -= fn * pio2_3;
+						w = fn * pio2_3t;
+						y[0] = r - w;
+					}
+				}
+			}
+		}
+		y[1] = (r - y[0]) - w;
+		if (hx < 0) {
+			y[0] = -y[0];
+			y[1] = -y[1];
+			return (-n);
+		}
+		return (n);
+	}
+
+	e0 = (ix >> 20) - 1046;	/* e0 = ilogb(x)-23; */
+
+	/* break x into three 24 bit pieces */
+	lx = ((int *)&x)[LOWORD];
+	i = (lx & 0x1f) << 19;
+	tx[2] = (double)i;
+	j = (lx >> 5) & 0xffffff;
+	tx[1] = (double)j;
+	tx[0] = (double)((((ix & 0xfffff) | 0x100000) << 3) |
+	    ((unsigned)lx >> 29));
+	nx = 3;
+	if (i == 0) {
+		/* skip zero term */
+		nx--;
+		if (j == 0)
+			nx--;
+	}
+	n = __rem_pio2m(tx, y, e0, nx, 2, _TBL_ipio2_inf);
+	if (hx < 0) {
+		y[0] = -y[0];
+		y[1] = -y[1];
+		return (-n);
+	}
+	return (n);
+}
