aboutsummaryrefslogtreecommitdiffstats
path: root/libm
diff options
context:
space:
mode:
Diffstat (limited to 'libm')
-rw-r--r--libm/i387/fenv.c3
-rw-r--r--libm/i387/fenv.h240
-rw-r--r--libm/include/i387/fenv.h5
-rw-r--r--libm/src/s_logb.c4
-rw-r--r--libm/src/s_remquo.c217
-rw-r--r--libm/src/s_remquof.c153
6 files changed, 192 insertions, 430 deletions
diff --git a/libm/i387/fenv.c b/libm/i387/fenv.c
index 2794faf81..aabe270d5 100644
--- a/libm/i387/fenv.c
+++ b/libm/i387/fenv.c
@@ -153,7 +153,8 @@ feholdexcept(fenv_t *envp)
int
feupdateenv(const fenv_t *envp)
{
- int mxcsr, status;
+ int mxcsr;
+ short status;
__fnstsw(&status);
if (__HAS_SSE())
diff --git a/libm/i387/fenv.h b/libm/i387/fenv.h
deleted file mode 100644
index b124366ac..000000000
--- a/libm/i387/fenv.h
+++ /dev/null
@@ -1,240 +0,0 @@
-/*-
- * Copyright (c) 2004-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.
- *
- * $FreeBSD: src/lib/msun/i387/fenv.h,v 1.4 2005/03/17 22:21:46 das Exp $
- */
-
-#ifndef _FENV_H_
-#define _FENV_H_
-
-#include <sys/cdefs.h>
-#include <sys/_types.h>
-
-/*
- * To preserve binary compatibility with FreeBSD 5.3, we pack the
- * mxcsr into some reserved fields, rather than changing sizeof(fenv_t).
- */
-typedef struct {
- __uint16_t __control;
- __uint16_t __mxcsr_hi;
- __uint16_t __status;
- __uint16_t __mxcsr_lo;
- __uint32_t __tag;
- char __other[16];
-} fenv_t;
-
-#define __get_mxcsr(env) (((env).__mxcsr_hi << 16) | \
- ((env).__mxcsr_lo))
-#define __set_mxcsr(env, x) do { \
- (env).__mxcsr_hi = (__uint32_t)(x) >> 16; \
- (env).__mxcsr_lo = (__uint16_t)(x); \
-} while (0)
-
-typedef __uint16_t fexcept_t;
-
-/* Exception flags */
-#define FE_INVALID 0x01
-#define FE_DENORMAL 0x02
-#define FE_DIVBYZERO 0x04
-#define FE_OVERFLOW 0x08
-#define FE_UNDERFLOW 0x10
-#define FE_INEXACT 0x20
-#define FE_ALL_EXCEPT (FE_DIVBYZERO | FE_DENORMAL | FE_INEXACT | \
- FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW)
-
-/* Rounding modes */
-#define FE_TONEAREST 0x0000
-#define FE_DOWNWARD 0x0400
-#define FE_UPWARD 0x0800
-#define FE_TOWARDZERO 0x0c00
-#define _ROUND_MASK (FE_TONEAREST | FE_DOWNWARD | \
- FE_UPWARD | FE_TOWARDZERO)
-
-/*
- * As compared to the x87 control word, the SSE unit's control word
- * has the rounding control bits offset by 3 and the exception mask
- * bits offset by 7.
- */
-#define _SSE_ROUND_SHIFT 3
-#define _SSE_EMASK_SHIFT 7
-
-/* After testing for SSE support once, we cache the result in __has_sse. */
-enum __sse_support { __SSE_YES, __SSE_NO, __SSE_UNK };
-extern enum __sse_support __has_sse;
-int __test_sse(void);
-#ifdef __SSE__
-#define __HAS_SSE() 1
-#else
-#define __HAS_SSE() (__has_sse == __SSE_YES || \
- (__has_sse == __SSE_UNK && __test_sse()))
-#endif
-
-__BEGIN_DECLS
-
-/* Default floating-point environment */
-extern const fenv_t __fe_dfl_env;
-#define FE_DFL_ENV (&__fe_dfl_env)
-
-#define __fldcw(__cw) __asm __volatile("fldcw %0" : : "m" (__cw))
-#define __fldenv(__env) __asm __volatile("fldenv %0" : : "m" (__env))
-#define __fnclex() __asm __volatile("fnclex")
-#define __fnstenv(__env) __asm __volatile("fnstenv %0" : "=m" (*(__env)))
-#define __fnstcw(__cw) __asm __volatile("fnstcw %0" : "=m" (*(__cw)))
-#define __fnstsw(__sw) __asm __volatile("fnstsw %0" : "=am" (*(__sw)))
-#define __fwait() __asm __volatile("fwait")
-#define __ldmxcsr(__csr) __asm __volatile("ldmxcsr %0" : : "m" (__csr))
-#define __stmxcsr(__csr) __asm __volatile("stmxcsr %0" : "=m" (*(__csr)))
-
-static __inline int
-feclearexcept(int __excepts)
-{
- fenv_t __env;
- int __mxcsr;
-
- if (__excepts == FE_ALL_EXCEPT) {
- __fnclex();
- } else {
- __fnstenv(&__env);
- __env.__status &= ~__excepts;
- __fldenv(__env);
- }
- if (__HAS_SSE()) {
- __stmxcsr(&__mxcsr);
- __mxcsr &= ~__excepts;
- __ldmxcsr(__mxcsr);
- }
- return (0);
-}
-
-static __inline int
-fegetexceptflag(fexcept_t *__flagp, int __excepts)
-{
- int __mxcsr, __status;
-
- __fnstsw(&__status);
- if (__HAS_SSE())
- __stmxcsr(&__mxcsr);
- else
- __mxcsr = 0;
- *__flagp = (__mxcsr | __status) & __excepts;
- return (0);
-}
-
-int fesetexceptflag(const fexcept_t *__flagp, int __excepts);
-int feraiseexcept(int __excepts);
-
-static __inline int
-fetestexcept(int __excepts)
-{
- int __mxcsr, __status;
-
- __fnstsw(&__status);
- if (__HAS_SSE())
- __stmxcsr(&__mxcsr);
- else
- __mxcsr = 0;
- return ((__status | __mxcsr) & __excepts);
-}
-
-static __inline int
-fegetround(void)
-{
- int __control;
-
- /*
- * We assume that the x87 and the SSE unit agree on the
- * rounding mode. Reading the control word on the x87 turns
- * out to be about 5 times faster than reading it on the SSE
- * unit on an Opteron 244.
- */
- __fnstcw(&__control);
- return (__control & _ROUND_MASK);
-}
-
-static __inline int
-fesetround(int __round)
-{
- int __mxcsr, __control;
-
- if (__round & ~_ROUND_MASK)
- return (-1);
-
- __fnstcw(&__control);
- __control &= ~_ROUND_MASK;
- __control |= __round;
- __fldcw(__control);
-
- if (__HAS_SSE()) {
- __stmxcsr(&__mxcsr);
- __mxcsr &= ~(_ROUND_MASK << _SSE_ROUND_SHIFT);
- __mxcsr |= __round << _SSE_ROUND_SHIFT;
- __ldmxcsr(__mxcsr);
- }
-
- return (0);
-}
-
-int fegetenv(fenv_t *__envp);
-int feholdexcept(fenv_t *__envp);
-
-static __inline int
-fesetenv(const fenv_t *__envp)
-{
- fenv_t __env = *__envp;
- int __mxcsr;
-
- __mxcsr = __get_mxcsr(__env);
- __set_mxcsr(__env, 0xffffffff);
- __fldenv(__env);
- if (__HAS_SSE())
- __ldmxcsr(__mxcsr);
- return (0);
-}
-
-int feupdateenv(const fenv_t *__envp);
-
-#if __BSD_VISIBLE
-
-int feenableexcept(int __mask);
-int fedisableexcept(int __mask);
-
-static __inline int
-fegetexcept(void)
-{
- int __control;
-
- /*
- * We assume that the masks for the x87 and the SSE unit are
- * the same.
- */
- __fnstcw(&__control);
- return (~__control & FE_ALL_EXCEPT);
-}
-
-#endif /* __BSD_VISIBLE */
-
-__END_DECLS
-
-#endif /* !_FENV_H_ */
diff --git a/libm/include/i387/fenv.h b/libm/include/i387/fenv.h
index b124366ac..4281f10eb 100644
--- a/libm/include/i387/fenv.h
+++ b/libm/include/i387/fenv.h
@@ -102,7 +102,7 @@ extern const fenv_t __fe_dfl_env;
#define __fnclex() __asm __volatile("fnclex")
#define __fnstenv(__env) __asm __volatile("fnstenv %0" : "=m" (*(__env)))
#define __fnstcw(__cw) __asm __volatile("fnstcw %0" : "=m" (*(__cw)))
-#define __fnstsw(__sw) __asm __volatile("fnstsw %0" : "=am" (*(__sw)))
+#define __fnstsw(__sw) __asm __volatile("fnstsw %0" : "=a" (*(__sw)))
#define __fwait() __asm __volatile("fwait")
#define __ldmxcsr(__csr) __asm __volatile("ldmxcsr %0" : : "m" (__csr))
#define __stmxcsr(__csr) __asm __volatile("stmxcsr %0" : "=m" (*(__csr)))
@@ -148,7 +148,8 @@ int feraiseexcept(int __excepts);
static __inline int
fetestexcept(int __excepts)
{
- int __mxcsr, __status;
+ int __mxcsr;
+ short __status;
__fnstsw(&__status);
if (__HAS_SSE())
diff --git a/libm/src/s_logb.c b/libm/src/s_logb.c
index 30edb8749..57f91c44d 100644
--- a/libm/src/s_logb.c
+++ b/libm/src/s_logb.c
@@ -36,9 +36,9 @@ logb(double x)
if(ix>=0x7ff00000) return x*x;
if(ix<0x00100000) {
x *= two54; /* convert subnormal x to normal */
- GET_FLOAT_WORD(ix,x);
+ GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff;
- return (float) ((ix>>20)-1023-54);
+ return (double) ((ix>>20)-1023-54);
} else
return (double) ((ix>>20)-1023);
}
diff --git a/libm/src/s_remquo.c b/libm/src/s_remquo.c
index eee65df4f..37d55e8f1 100644
--- a/libm/src/s_remquo.c
+++ b/libm/src/s_remquo.c
@@ -5,7 +5,7 @@
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@@ -29,124 +29,123 @@ static const double Zero[] = {0.0, -0.0,};
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;
+ 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| */
+ 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*/
- }
- }
+ 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 = (sxy ? -1 : 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;
+ 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;
+ 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;
- }
- }
-
+ 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++;}
-
+ 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;
- }
- }
+ if((hx|lx)==0) { /* return sign(x)*0 */
+ q &= 0x7fffffff;
+ *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 = 0;
+ } else {
+ lx = hx>>(n-32); hx = 0;
+ }
+ }
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;
+ 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/libm/src/s_remquof.c b/libm/src/s_remquof.c
index 5d722ceaf..36269a622 100644
--- a/libm/src/s_remquof.c
+++ b/libm/src/s_remquof.c
@@ -5,7 +5,7 @@
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@@ -29,93 +29,94 @@ static const float Zero[] = {0.0, -0.0,};
float
remquof(float x, float y, int *quo)
{
- int32_t n,hx,hy,hz,ix,iy,sx,i;
- u_int32_t q,sxy;
+ 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| */
+ 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*/
- }
+ 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 = (sxy ? -1 : 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;
+ 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;
+ 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;
- }
+ 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++;}
+ 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;
- }
+ if(hx==0) { /* return sign(x)*0 */
+ q &= 0x7fffffff;
+ *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;
+ 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;
}
generated by cgit v1.2.3 (git 2.25.1) at 2025-06-17 03:57:24 +0000