d/freebsd-skq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Implement and document remquo() and remquof().

Browse Source
This commit is contained in:
David Schultz 2005-03-25 04:40:44 +00:00
parent 04f716a436
commit 3b9141ee91
12 changed files with 583 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
lib/msun/Makefile
View File

@ -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

2
lib/msun/amd64/Makefile.inc
View File

@ -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

65
lib/msun/amd64/s_remquo.S Normal file
View File

@ -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

65
lib/msun/amd64/s_remquof.S Normal file
View File

@ -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

4
lib/msun/i387/Makefile.inc
View File

@ -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

62
lib/msun/i387/s_remquo.S Normal file
View File

@ -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

62
lib/msun/i387/s_remquof.S Normal file
View File

@ -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

7
lib/msun/man/math.3
View File

@ -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.

57
lib/msun/man/remainder.3
View File

@ -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 .

2
lib/msun/src/math.h
View File

@ -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);

152
lib/msun/src/s_remquo.c Normal file
View File

@ -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;
}

121
lib/msun/src/s_remquof.c Normal file
View File

@ -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;
}