diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/e_remainder.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_remainder.c | 175 |
1 files changed, 103 insertions, 72 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c index 6418081182..631c2fcf01 100644 --- a/sysdeps/ieee754/dbl-64/e_remainder.c +++ b/sysdeps/ieee754/dbl-64/e_remainder.c @@ -1,80 +1,111 @@ -/* @(#)e_remainder.c 5.1 93/09/24 */ /* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * IBM Accurate Mathematical Library + * Copyright (c) International Business Machines Corp., 2001 * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_remainder.c,v 1.8 1995/05/10 20:46:05 jtc Exp $"; -#endif - -/* __ieee754_remainder(x,p) - * Return : - * returns x REM p = x - [x/p]*p as if in infinite - * precise arithmetic, where [x/p] is the (infinite bit) - * integer nearest x/p (in half way case choose the even one). - * Method : - * Based on fmod() return x-[x/p]chopped*p exactlp. + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +/**************************************************************************/ +/* MODULE_NAME urem.c */ +/* */ +/* FUNCTION: uremainder */ +/* */ +/* An ultimate remainder routine. Given two IEEE double machine numbers x */ +/* ,y it computes the correctly rounded (to nearest) value of remainder */ +/* of dividing x by y. */ +/* Assumption: Machine arithmetic operations are performed in */ +/* round to nearest mode of IEEE 754 standard. */ +/* */ +/* ************************************************************************/ -#include "math.h" -#include "math_private.h" +#include "endian.h" +#include "mydefs.h" +#include "urem.h" +#include "MathLib.h" -#ifdef __STDC__ -static const double zero = 0.0; -#else -static double zero = 0.0; -#endif - - -#ifdef __STDC__ - double __ieee754_remainder(double x, double p) -#else - double __ieee754_remainder(x,p) - double x,p; -#endif +/**************************************************************************/ +/* An ultimate remainder routine. Given two IEEE double machine numbers x */ +/* ,y it computes the correctly rounded (to nearest) value of remainder */ +/**************************************************************************/ +double __ieee754_remainder(double x, double y) { - int32_t hx,hp; - u_int32_t sx,lx,lp; - double p_half; - - EXTRACT_WORDS(hx,lx,x); - EXTRACT_WORDS(hp,lp,p); - sx = hx&0x80000000; - hp &= 0x7fffffff; - hx &= 0x7fffffff; - - /* purge off exception values */ - if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */ - if((hx>=0x7ff00000)|| /* x not finite */ - ((hp>=0x7ff00000)&& /* p is NaN */ - (((hp-0x7ff00000)|lp)!=0))) - return (x*p)/(x*p); - + double z,d,xx,yy; + int4 kx,ky,m,n,nn,n1,m1,l; + mynumber u,t,w={0,0},v={0,0},ww={0,0},r; + u.x=x; + t.x=y; + kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign for x*/ + t.i[HIGH_HALF]&=0x7fffffff; /*no sign for y */ + ky=t.i[HIGH_HALF]; + /*------ |x| < 2^1024 and 2^-970 < |y| < 2^1024 ------------------*/ + if (kx<0x7ff00000 && ky<0x7ff00000 && ky>=0x03500000) { + if (kx+0x00100000<ky) return x; + if ((kx-0x01500000)<ky) { + z=x/t.x; + v.i[HIGH_HALF]=t.i[HIGH_HALF]; + d=(z+big.x)-big.x; + xx=(x-d*v.x)-d*(t.x-v.x); + if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x); + else { + if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x; + else return xx; + } + } /* (kx<(ky+0x01500000)) */ + else { + r.x=1.0/t.x; + n=t.i[HIGH_HALF]; + nn=(n&0x7ff00000)+0x01400000; + w.i[HIGH_HALF]=n; + ww.x=t.x-w.x; + l=(kx-nn)&0xfff00000; + n1=ww.i[HIGH_HALF]; + m1=r.i[HIGH_HALF]; + while (l>0) { + r.i[HIGH_HALF]=m1-l; + z=u.x*r.x; + w.i[HIGH_HALF]=n+l; + ww.i[HIGH_HALF]=(n1)?n1+l:n1; + d=(z+big.x)-big.x; + u.x=(u.x-d*w.x)-d*ww.x; + l=(u.i[HIGH_HALF]&0x7ff00000)-nn; + } + r.i[HIGH_HALF]=m1; + w.i[HIGH_HALF]=n; + ww.i[HIGH_HALF]=n1; + z=u.x*r.x; + d=(z+big.x)-big.x; + u.x=(u.x-d*w.x)-d*ww.x; + if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x); + else + if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x; + else + {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);} + } - if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */ - if (((hx-hp)|(lx-lp))==0) return zero*x; - x = fabs(x); - p = fabs(p); - if (hp<0x00200000) { - if(x+x>p) { - x-=p; - if(x+x>=p) x -= p; - } - } else { - p_half = 0.5*p; - if(x>p_half) { - x-=p; - if(x>=p_half) x -= p; - } - } - GET_HIGH_WORD(hx,x); - SET_HIGH_WORD(x,hx^sx); - return x; + } /* (kx<0x7ff00000&&ky<0x7ff00000&&ky>=0x03500000) */ + else { + if (kx<0x7ff00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) { + y=ABS(y)*t128.x; + z=uremainder(x,y)*t128.x; + z=uremainder(z,y)*tm128.x; + return z; + } + else { /* if x is too big */ + if (kx>=0x7ff00000||(ky==0&&t.i[LOW_HALF]==0)||ky>0x7ff00000|| + (ky==0x7ff00000&&t.i[LOW_HALF]!=0)) + return (u.i[HIGH_HALF]&0x80000000)?nNAN.x:NAN.x; + else return x; + } + } } |