From 3b9141ee9171245e22c7b4becddf194699e25845 Mon Sep 17 00:00:00 2001 From: David Schultz Date: Fri, 25 Mar 2005 04:40:44 +0000 Subject: [PATCH] Implement and document remquo() and remquof(). --- lib/msun/Makefile | 4 +- lib/msun/amd64/Makefile.inc | 2 +- lib/msun/amd64/s_remquo.S | 65 +++++++++++++++ lib/msun/amd64/s_remquof.S | 65 +++++++++++++++ lib/msun/i387/Makefile.inc | 4 +- lib/msun/i387/s_remquo.S | 62 +++++++++++++++ lib/msun/i387/s_remquof.S | 62 +++++++++++++++ lib/msun/man/math.3 | 7 +- lib/msun/man/remainder.3 | 57 +++++++++++--- lib/msun/src/math.h | 2 + lib/msun/src/s_remquo.c | 152 ++++++++++++++++++++++++++++++++++++ lib/msun/src/s_remquof.c | 121 ++++++++++++++++++++++++++++ 12 files changed, 583 insertions(+), 20 deletions(-) create mode 100644 lib/msun/amd64/s_remquo.S create mode 100644 lib/msun/amd64/s_remquof.S create mode 100644 lib/msun/i387/s_remquo.S create mode 100644 lib/msun/i387/s_remquof.S create mode 100644 lib/msun/src/s_remquo.c create mode 100644 lib/msun/src/s_remquof.c diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 7c0f56baaeca..63b279268bab 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 8f7573889bc2..a747b678424d 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 000000000000..eb113d718e86 --- /dev/null +++ b/lib/msun/amd64/s_remquo.S @@ -0,0 +1,65 @@ +/*- + * Copyright (c) 2005 David Schultz + * 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 . + */ + +#include +__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 000000000000..0833f5b984b5 --- /dev/null +++ b/lib/msun/amd64/s_remquof.S @@ -0,0 +1,65 @@ +/*- + * Copyright (c) 2005 David Schultz + * 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 . + */ + +#include +__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 3be260ec2bb1..38cb0fc2e2d0 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 000000000000..fd8af5ef8b3c --- /dev/null +++ b/lib/msun/i387/s_remquo.S @@ -0,0 +1,62 @@ +/*- + * Copyright (c) 2005 David Schultz + * 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 . + */ + +#include +__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 000000000000..70a43f4b4b39 --- /dev/null +++ b/lib/msun/i387/s_remquof.S @@ -0,0 +1,62 @@ +/*- + * Copyright (c) 2005 David Schultz + * 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 . + */ + +#include +__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 0a5ad4112fd1..142f75e99a4a 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 4b5b877838bd..10ac4dedcf9e 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 9b8f73cb9787..b29ace27b646 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 000000000000..8756c2a3d0bb --- /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 +__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>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<>(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<>(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>31); lx = lx+lx;} + else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} + q <<= 1; + } + hz=hx-hy;lz=lx-ly; if(lx=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 000000000000..3141642116a2 --- /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 +__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>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; +}