Implement and document remquo() and remquof().

12 files changed, 583 insertions, 20 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 7c0f56b..63b2792 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -50,7 +50,8 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
 	s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_lrint.c s_lrintf.c \
 	s_lround.c s_lroundf.c s_modff.c \
 	s_nearbyint.c s_nextafter.c s_nextafterf.c \
-	s_nexttowardf.c s_rint.c s_rintf.c s_round.c s_roundf.c \
+	s_nexttowardf.c s_remquo.c s_remquof.c \
+	s_rint.c s_rintf.c s_round.c s_roundf.c \
 	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
 	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c s_tan.c \
 	s_tanf.c s_tanh.c s_tanhf.c s_trunc.c s_truncf.c \
@@ -140,6 +141,7 @@ MLINKS+=nextafter.3 nextafterf.3 nextafter.3 nextafterl.3
 MLINKS+=nextafter.3 nexttoward.3 nextafter.3 nexttowardf.3
 MLINKS+=nextafter.3 nexttowardl.3
 MLINKS+=remainder.3 remainderf.3
+MLINKS+=remainder.3 remquo.3 remainder.3 remquof.3
 MLINKS+=rint.3 rintf.3 rint.3 nearbyint.3 rint.3 nearbyintf.3
 MLINKS+=round.3 roundf.3
 MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3
diff --git a/lib/msun/amd64/Makefile.inc b/lib/msun/amd64/Makefile.inc
index 8f75738..a747b67 100644
--- a/lib/msun/amd64/Makefile.inc
+++ b/lib/msun/amd64/Makefile.inc
@@ -1,4 +1,4 @@
 # $FreeBSD$
 
-ARCH_SRCS = e_sqrt.S s_lrint.S s_llrint.S
+ARCH_SRCS = e_sqrt.S s_lrint.S s_llrint.S s_remquo.S s_remquof.S
 LDBL_PREC = 64
diff --git a/lib/msun/amd64/s_remquo.S b/lib/msun/amd64/s_remquo.S
new file mode 100644
index 0000000..eb113d7
--- /dev/null
+++ b/lib/msun/amd64/s_remquo.S
@@ -0,0 +1,65 @@
+/*-
+ * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/*
+ * Based on public-domain remainder routine by J.T. Conklin <jtc@NetBSD.org>.
+ */
+
+#include <machine/asm.h>
+__FBSDID("$FreeBSD$"); 
+	
+ENTRY(remquo)
+	movsd	%xmm0,-8(%rsp)
+	movsd	%xmm1,-16(%rsp)
+	fldl	-16(%rsp)
+	fldl	-8(%rsp)
+1:	fprem1
+	fstsw	%ax
+	btw	$10,%ax
+	jc	1b
+	fstp	%st(1)
+/* Extract the three low-order bits of the quotient from C0,C3,C1. */
+	shrl	$6,%eax
+	movl	%eax,%ecx
+	andl	$0x108,%eax
+	rorl	$7,%eax
+	orl	%eax,%ecx
+	roll	$4,%eax
+	orl	%ecx,%eax
+	andl	$7,%eax
+/* Negate the quotient bits if x*y<0.  Avoid using an unpredictable branch. */
+	movl	-12(%rsp),%ecx
+	xorl	-4(%rsp),%ecx
+	sarl	$16,%ecx
+	sarl	$16,%ecx
+	xorl	%ecx,%eax
+	andl	$1,%ecx
+	addl	%ecx,%eax
+/* Store the quotient and return. */
+	movl	%eax,(%rdi)
+	fstpl	-8(%rsp)
+	movsd	-8(%rsp),%xmm0
+	ret
diff --git a/lib/msun/amd64/s_remquof.S b/lib/msun/amd64/s_remquof.S
new file mode 100644
index 0000000..0833f5b
--- /dev/null
+++ b/lib/msun/amd64/s_remquof.S
@@ -0,0 +1,65 @@
+/*-
+ * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/*
+ * Based on public-domain remainder routine by J.T. Conklin <jtc@NetBSD.org>.
+ */
+
+#include <machine/asm.h>
+__FBSDID("$FreeBSD$"); 
+	
+ENTRY(remquof)
+	movss	%xmm0,-4(%rsp)
+	movss	%xmm1,-8(%rsp)
+	flds	-8(%rsp)
+	flds	-4(%rsp)
+1:	fprem1
+	fstsw	%ax
+	btw	$10,%ax
+	jc	1b
+	fstp	%st(1)
+/* Extract the three low-order bits of the quotient from C0,C3,C1. */
+	shrl	$6,%eax
+	movl	%eax,%ecx
+	andl	$0x108,%eax
+	rorl	$7,%eax
+	orl	%eax,%ecx
+	roll	$4,%eax
+	orl	%ecx,%eax
+	andl	$7,%eax
+/* Negate the quotient bits if x*y<0.  Avoid using an unpredictable branch. */
+	movl	-8(%rsp),%ecx
+	xorl	-4(%rsp),%ecx
+	sarl	$16,%ecx
+	sarl	$16,%ecx
+	xorl	%ecx,%eax
+	andl	$1,%ecx
+	addl	%ecx,%eax
+/* Store the quotient and return. */
+	movl	%eax,(%rdi)
+	fstps	-4(%rsp)
+	movss	-4(%rsp),%xmm0
+	ret
diff --git a/lib/msun/i387/Makefile.inc b/lib/msun/i387/Makefile.inc
index 3be260e..38cb0fc 100644
--- a/lib/msun/i387/Makefile.inc
+++ b/lib/msun/i387/Makefile.inc
@@ -3,12 +3,12 @@
 ARCH_SRCS = e_exp.S e_fmod.S e_log.S e_log10.S \
 	    e_remainder.S e_scalb.S e_sqrt.S s_ceil.S s_copysign.S \
 	    s_cos.S s_finite.S s_floor.S s_llrint.S s_logb.S s_lrint.S \
-	    s_rint.S s_scalbn.S s_significand.S s_sin.S s_tan.S
+	    s_remquo.S s_rint.S s_scalbn.S s_significand.S s_sin.S s_tan.S
 
 # float counterparts
 ARCH_SRCS+= e_log10f.S e_logf.S e_remainderf.S e_scalbf.S \
 	    e_sqrtf.S s_ceilf.S s_copysignf.S s_floorf.S s_logbf.S \
-	    s_rintf.S s_scalbnf.S s_significandf.S
+	    s_remquof.S s_rintf.S s_scalbnf.S s_significandf.S
 
 # long double counterparts
 ARCH_SRCS+= s_scalbnl.S
diff --git a/lib/msun/i387/s_remquo.S b/lib/msun/i387/s_remquo.S
new file mode 100644
index 0000000..fd8af5e
--- /dev/null
+++ b/lib/msun/i387/s_remquo.S
@@ -0,0 +1,62 @@
+/*-
+ * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/*
+ * Based on public-domain remainder routine by J.T. Conklin <jtc@NetBSD.org>.
+ */
+
+#include <machine/asm.h>
+__FBSDID("$FreeBSD$"); 
+	
+ENTRY(remquo)
+	fldl	12(%esp)
+	fldl	4(%esp)
+1:	fprem1
+	fstsw	%ax
+	sahf
+	jp	1b
+	fstp	%st(1)
+/* Extract the three low-order bits of the quotient from C0,C3,C1. */
+	shrl	$6,%eax
+	movl	%eax,%ecx
+	andl	$0x108,%eax
+	rorl	$7,%eax
+	orl	%eax,%ecx
+	roll	$4,%eax
+	orl	%ecx,%eax
+	andl	$7,%eax
+/* Negate the quotient bits if x*y<0.  Avoid using an unpredictable branch. */
+	movl	16(%esp),%ecx
+	xorl	8(%esp),%ecx
+	sarl	$16,%ecx
+	sarl	$16,%ecx
+	xorl	%ecx,%eax
+	andl	$1,%ecx
+	addl	%ecx,%eax
+/* Store the quotient and return. */
+	movl	20(%esp),%ecx
+	movl	%eax,(%ecx)
+	ret
diff --git a/lib/msun/i387/s_remquof.S b/lib/msun/i387/s_remquof.S
new file mode 100644
index 0000000..70a43f4
--- /dev/null
+++ b/lib/msun/i387/s_remquof.S
@@ -0,0 +1,62 @@
+/*-
+ * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/*
+ * Based on public-domain remainder routine by J.T. Conklin <jtc@NetBSD.org>.
+ */
+
+#include <machine/asm.h>
+__FBSDID("$FreeBSD$"); 
+	
+ENTRY(remquof)
+	flds	8(%esp)
+	flds	4(%esp)
+1:	fprem1
+	fstsw	%ax
+	sahf
+	jp	1b
+	fstp	%st(1)
+/* Extract the three low-order bits of the quotient from C0,C3,C1. */
+	shrl	$6,%eax
+	movl	%eax,%ecx
+	andl	$0x108,%eax
+	rorl	$7,%eax
+	orl	%eax,%ecx
+	roll	$4,%eax
+	orl	%ecx,%eax
+	andl	$7,%eax
+/* Negate the quotient bits if x*y<0.  Avoid using an unpredictable branch. */
+	movl	8(%esp),%ecx
+	xorl	4(%esp),%ecx
+	sarl	$16,%ecx
+	sarl	$16,%ecx
+	xorl	%ecx,%eax
+	andl	$1,%ecx
+	addl	%ecx,%eax
+/* Store the quotient and return. */
+	movl	12(%esp),%ecx
+	movl	%eax,(%ecx)
+	ret
diff --git a/lib/msun/man/math.3 b/lib/msun/man/math.3
index 0a5ad41..142f75e 100644
--- a/lib/msun/man/math.3
+++ b/lib/msun/man/math.3
@@ -32,7 +32,7 @@
 .\"	from: @(#)math.3	6.10 (Berkeley) 5/6/91
 .\" $FreeBSD$
 .\"
-.Dd January 11, 2005
+.Dd March 24, 2005
 .Dt MATH 3
 .Os
 .if n \{\
@@ -127,7 +127,7 @@ nearbyint	round to integer (silent)
 nextafter	next representable value
 nexttoward	next representable value
 remainder	remainder
-.\" remquo	remainder with partial quotient
+remquo	remainder with partial quotient
 rint	round to integer
 round	round to nearest integer
 trunc	integer no greater in magnitude than
@@ -223,9 +223,8 @@ values, were written for or imported into subsequent versions of FreeBSD.
 The
 .Fn exp2 ,
 .Fn log2 ,
-.Fn nan ,
 and
-.Fn remquo
+.Fn nan
 functions are missing, and many functions are not available in their
 .Vt "long double"
 variants.
diff --git a/lib/msun/man/remainder.3 b/lib/msun/man/remainder.3
index 4b5b877..10ac4de 100644
--- a/lib/msun/man/remainder.3
+++ b/lib/msun/man/remainder.3
@@ -32,12 +32,14 @@
 .\"     from: @(#)ieee.3	6.4 (Berkeley) 5/6/91
 .\" $FreeBSD$
 .\"
-.Dd January 26, 2005
+.Dd March 24, 2005
 .Dt REMAINDER 3
 .Os
 .Sh NAME
 .Nm remainder ,
-.Nm remainderf
+.Nm remainderf ,
+.Nm remquo ,
+.Nm remquof
 .Nd minimal residue functions
 .Sh LIBRARY
 .Lb libm
@@ -47,10 +49,16 @@
 .Fn remainder "double x" "double y"
 .Ft float
 .Fn remainderf "float x" "float y"
+.Ft double
+.Fn remquo "double x" "double y" "int *quo"
+.Ft float
+.Fn remquo "float x" "float y" "int *quo"
 .Sh DESCRIPTION
-.Fn remainder
+.Fn remainder ,
+.Fn remainderf ,
+.Fn remquo ,
 and
-.Fn remainderf
+.Fn remquof
 return the remainder
 .Fa r
 :=
@@ -83,23 +91,42 @@ the remainder is computed exactly and
 .Sm off
 .Pf \\*(Ba Fa y No \\*(Ba/2 .
 .Sm on
-But
-.Fn remainder x 0
+But attempting to take the remainder when
+.Fa y
+is 0 or
+.Fa x
+is \*(Pm\*(If is an invalid operation that produces a \*(Na.
+.Pp
+The
+.Fn remquo
 and
-.Fn remainder \*(If 0
-are invalid operations that produce a \*(Na.
+.Fn remquof
+functions also store the last
+.Va k
+bits of
+.Fa n
+in the location pointed to by
+.Fa quo ,
+provided that
+.Fa n
+exists.
+The number of bits
+.Va k
+is platform-specific, but is guaranteed to be at least 3.
 .Sh SEE ALSO
 .Xr fmod 3 ,
 .Xr ieee 3 ,
 .Xr math 3
 .Sh STANDARDS
 The
-.Fn remainder
+.Fn remainder ,
+.Fn remainderf ,
+.Fn remquo ,
 and
-.Fn remainderf
+.Fn remquof
 routines conform to
-.St -isoC-99
-and
+.St -isoC-99 .
+The remainder is as defined in
 .St -ieee754 .
 .Sh HISTORY
 The
@@ -111,3 +138,9 @@ functions appeared in
 and
 .Fx 2.0 ,
 respectively.
+The
+.Fn remquo
+and
+.Fn remquof
+functions were added in
+.Fx 5.5 .
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index 9b8f73c..b29ace2 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -239,6 +239,7 @@ long	lrint(double);
 long	lround(double);
 double	nextafter(double, double);
 double	remainder(double, double);
+double	remquo(double, double, int *);
 double	rint(double);
 #endif /* __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999 || __XSI_VISIBLE */
 
@@ -341,6 +342,7 @@ long	lroundf(float);
 float	nearbyintf(float);
 float	nextafterf(float, float);
 float	remainderf(float, float);
+float	remquof(float, float, int *);
 float	rintf(float);
 float	scalblnf(float, long);
 float	scalbnf(float, int);
diff --git a/lib/msun/src/s_remquo.c b/lib/msun/src/s_remquo.c
new file mode 100644
index 0000000..8756c2a
--- /dev/null
+++ b/lib/msun/src/s_remquo.c
@@ -0,0 +1,152 @@
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include "math.h"
+#include "math_private.h"
+
+static const double Zero[] = {0.0, -0.0,};
+
+/*
+ * Return the IEEE remainder and set *quo to the last n bits of the
+ * quotient, rounded to the nearest integer.  We choose n=31 because
+ * we wind up computing all the integer bits of the quotient anyway as
+ * a side-effect of computing the remainder by the shift and subtract
+ * method.  In practice, this is far more bits than are needed to use
+ * remquo in reduction algorithms.
+ */
+double
+remquo(double x, double y, int *quo)
+{
+	int32_t n,hx,hy,hz,ix,iy,sx,i;
+	u_int32_t lx,ly,lz,q,sxy;
+
+	EXTRACT_WORDS(hx,lx,x);
+	EXTRACT_WORDS(hy,ly,y);
+	sxy = (hx ^ hy) & 0x80000000;
+	sx = hx&0x80000000;		/* sign of x */
+	hx ^=sx;		/* |x| */
+	hy &= 0x7fffffff;	/* |y| */
+
+    /* purge off exception values */
+	if((hy|ly)==0||(hx>=0x7ff00000)||	/* y=0,or x not finite */
+	  ((hy|((ly|-ly)>>31))>0x7ff00000))	/* or y is NaN */
+	    return (x*y)/(x*y);
+	if(hx<=hy) {
+	    if((hx<hy)||(lx<ly)) {
+		q = 0;
+		goto fixup;	/* |x|<|y| return x or x-y */
+	    }
+	    if(lx==ly) {
+		*quo = 1;
+		return Zero[(u_int32_t)sx>>31];	/* |x|=|y| return x*0*/
+	    }
+	}
+
+    /* determine ix = ilogb(x) */
+	if(hx<0x00100000) {	/* subnormal x */
+	    if(hx==0) {
+		for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
+	    } else {
+		for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
+	    }
+	} else ix = (hx>>20)-1023;
+
+    /* determine iy = ilogb(y) */
+	if(hy<0x00100000) {	/* subnormal y */
+	    if(hy==0) {
+		for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
+	    } else {
+		for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
+	    }
+	} else iy = (hy>>20)-1023;
+
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+	if(ix >= -1022) 
+	    hx = 0x00100000|(0x000fffff&hx);
+	else {		/* subnormal x, shift x to normal */
+	    n = -1022-ix;
+	    if(n<=31) {
+	        hx = (hx<<n)|(lx>>(32-n));
+	        lx <<= n;
+	    } else {
+		hx = lx<<(n-32);
+		lx = 0;
+	    }
+	}
+	if(iy >= -1022) 
+	    hy = 0x00100000|(0x000fffff&hy);
+	else {		/* subnormal y, shift y to normal */
+	    n = -1022-iy;
+	    if(n<=31) {
+	        hy = (hy<<n)|(ly>>(32-n));
+	        ly <<= n;
+	    } else {
+		hy = ly<<(n-32);
+		ly = 0;
+	    }
+	}
+
+    /* fix point fmod */
+	n = ix - iy;
+	q = 0;
+	while(n--) {
+	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	    if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
+	    else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;}
+	    q <<= 1;
+	}
+	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	if(hz>=0) {hx=hz;lx=lz;q++;}
+
+    /* convert back to floating value and restore the sign */
+	if((hx|lx)==0) {			/* return sign(x)*0 */
+	    *quo = (sxy ? -q : q);
+	    return Zero[(u_int32_t)sx>>31];
+	}
+	while(hx<0x00100000) {		/* normalize x */
+	    hx = hx+hx+(lx>>31); lx = lx+lx;
+	    iy -= 1;
+	}
+	if(iy>= -1022) {	/* normalize output */
+	    hx = ((hx-0x00100000)|((iy+1023)<<20));
+	} else {		/* subnormal output */
+	    n = -1022 - iy;
+	    if(n<=20) {
+		lx = (lx>>n)|((u_int32_t)hx<<(32-n));
+		hx >>= n;
+	    } else if (n<=31) {
+		lx = (hx<<(32-n))|(lx>>n); hx = sx;
+	    } else {
+		lx = hx>>(n-32); hx = sx;
+	    }
+	}
+fixup:
+	INSERT_WORDS(x,hx,lx);
+	y = fabs(y);
+	if (y < 0x1p-1021) {
+	    if (x+x>y || (x+x==y && (q & 1))) {
+		q++;
+		x-=y;
+	    }
+	} else if (x>0.5*y || (x==0.5*y && (q & 1))) {
+	    q++;
+	    x-=y;
+	}
+	GET_HIGH_WORD(hx,x);
+	SET_HIGH_WORD(x,hx^sx);
+	q &= 0x7fffffff;
+	*quo = (sxy ? -q : q);
+	return x;
+}
diff --git a/lib/msun/src/s_remquof.c b/lib/msun/src/s_remquof.c
new file mode 100644
index 0000000..3141642
--- /dev/null
+++ b/lib/msun/src/s_remquof.c
@@ -0,0 +1,121 @@
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include "math.h"
+#include "math_private.h"
+
+static const float Zero[] = {0.0, -0.0,};
+
+/*
+ * Return the IEEE remainder and set *quo to the last n bits of the
+ * quotient, rounded to the nearest integer.  We choose n=31 because
+ * we wind up computing all the integer bits of the quotient anyway as
+ * a side-effect of computing the remainder by the shift and subtract
+ * method.  In practice, this is far more bits than are needed to use
+ * remquo in reduction algorithms.
+ */
+float
+remquof(float x, float y, int *quo)
+{
+	int32_t n,hx,hy,hz,ix,iy,sx,i;
+	u_int32_t q,sxy;
+
+	GET_FLOAT_WORD(hx,x);
+	GET_FLOAT_WORD(hy,y);
+	sxy = (hx ^ hy) & 0x80000000;
+	sx = hx&0x80000000;		/* sign of x */
+	hx ^=sx;		/* |x| */
+	hy &= 0x7fffffff;	/* |y| */
+
+    /* purge off exception values */
+	if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */
+	    return (x*y)/(x*y);
+	if(hx<hy) {
+	    q = 0;
+	    goto fixup;	/* |x|<|y| return x or x-y */
+	} else if(hx==hy) {
+	    *quo = 1;
+	    return Zero[(u_int32_t)sx>>31];	/* |x|=|y| return x*0*/
+	}
+
+    /* determine ix = ilogb(x) */
+	if(hx<0x00800000) {	/* subnormal x */
+	    for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
+	} else ix = (hx>>23)-127;
+
+    /* determine iy = ilogb(y) */
+	if(hy<0x00800000) {	/* subnormal y */
+	    for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1;
+	} else iy = (hy>>23)-127;
+
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+	if(ix >= -126)
+	    hx = 0x00800000|(0x007fffff&hx);
+	else {		/* subnormal x, shift x to normal */
+	    n = -126-ix;
+	    hx <<= n;
+	}
+	if(iy >= -126)
+	    hy = 0x00800000|(0x007fffff&hy);
+	else {		/* subnormal y, shift y to normal */
+	    n = -126-iy;
+	    hy <<= n;
+	}
+
+    /* fix point fmod */
+	n = ix - iy;
+	q = 0;
+	while(n--) {
+	    hz=hx-hy;
+	    if(hz<0) hx = hx << 1;
+	    else {hx = hz << 1; q++;}
+	    q <<= 1;
+	}
+	hz=hx-hy;
+	if(hz>=0) {hx=hz;q++;}
+
+    /* convert back to floating value and restore the sign */
+	if(hx==0) {				/* return sign(x)*0 */
+	    *quo = (sxy ? -q : q);
+	    return Zero[(u_int32_t)sx>>31];
+	}
+	while(hx<0x00800000) {		/* normalize x */
+	    hx <<= 1;
+	    iy -= 1;
+	}
+	if(iy>= -126) {		/* normalize output */
+	    hx = ((hx-0x00800000)|((iy+127)<<23));
+	} else {		/* subnormal output */
+	    n = -126 - iy;
+	    hx >>= n;
+	}
+fixup:
+	SET_FLOAT_WORD(x,hx);
+	y = fabsf(y);
+	if (y < 0x1p-125f) {
+	    if (x+x>y || (x+x==y && (q & 1))) {
+		q++;
+		x-=y;
+	    }
+	} else if (x>0.5f*y || (x==0.5f*y && (q & 1))) {
+	    q++;
+	    x-=y;
+	}
+	GET_FLOAT_WORD(hx,x);
+	SET_FLOAT_WORD(x,hx^sx);
+	q &= 0x7fffffff;
+	*quo = (sxy ? -q : q);
+	return x;
+}