Back to index

glibc  2.9
e_remainder.c
Go to the documentation of this file.
00001 /*
00002  * IBM Accurate Mathematical Library
00003  * written by International Business Machines Corp.
00004  * Copyright (C) 2001 Free Software Foundation
00005  * 
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU Lesser General Public License as published by
00008  * the Free Software Foundation; either version 2.1 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00019  */
00020 /**************************************************************************/
00021 /*  MODULE_NAME urem.c                                                    */
00022 /*                                                                        */
00023 /*  FUNCTION: uremainder                                                  */
00024 /*                                                                        */
00025 /* An ultimate remainder routine. Given two IEEE double machine numbers x */
00026 /* ,y   it computes the correctly rounded (to nearest) value of remainder */
00027 /* of dividing x by y.                                                    */
00028 /* Assumption: Machine arithmetic operations are performed in             */
00029 /* round to nearest mode of IEEE 754 standard.                            */
00030 /*                                                                        */
00031 /* ************************************************************************/
00032 
00033 #include "endian.h"
00034 #include "mydefs.h"
00035 #include "urem.h"
00036 #include "MathLib.h"
00037 #include "math_private.h"
00038 
00039 /**************************************************************************/
00040 /* An ultimate remainder routine. Given two IEEE double machine numbers x */
00041 /* ,y   it computes the correctly rounded (to nearest) value of remainder */
00042 /**************************************************************************/
00043 double __ieee754_remainder(double x, double y)
00044 {
00045   double z,d,xx;
00046 #if 0
00047   double yy;
00048 #endif
00049   int4 kx,ky,n,nn,n1,m1,l;
00050 #if 0
00051   int4 m;
00052 #endif
00053   mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r;
00054   u.x=x;
00055   t.x=y;
00056   kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign  for x*/
00057   t.i[HIGH_HALF]&=0x7fffffff;   /*no sign for y */
00058   ky=t.i[HIGH_HALF];
00059   /*------ |x| < 2^1023  and   2^-970 < |y| < 2^1024 ------------------*/
00060   if (kx<0x7fe00000 && ky<0x7ff00000 && ky>=0x03500000) {
00061     if (kx+0x00100000<ky) return x;
00062     if ((kx-0x01500000)<ky) {
00063       z=x/t.x;
00064       v.i[HIGH_HALF]=t.i[HIGH_HALF];
00065       d=(z+big.x)-big.x;
00066       xx=(x-d*v.x)-d*(t.x-v.x);
00067       if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x);
00068       else {
00069        if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x;
00070        else return xx;
00071       }
00072     }   /*    (kx<(ky+0x01500000))         */
00073     else  {
00074       r.x=1.0/t.x;
00075       n=t.i[HIGH_HALF];
00076       nn=(n&0x7ff00000)+0x01400000;
00077       w.i[HIGH_HALF]=n;
00078       ww.x=t.x-w.x;
00079       l=(kx-nn)&0xfff00000;
00080       n1=ww.i[HIGH_HALF];
00081       m1=r.i[HIGH_HALF];
00082       while (l>0) {
00083        r.i[HIGH_HALF]=m1-l;
00084        z=u.x*r.x;
00085        w.i[HIGH_HALF]=n+l;
00086        ww.i[HIGH_HALF]=(n1)?n1+l:n1;
00087        d=(z+big.x)-big.x;
00088        u.x=(u.x-d*w.x)-d*ww.x;
00089        l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
00090       }
00091       r.i[HIGH_HALF]=m1;
00092       w.i[HIGH_HALF]=n;
00093       ww.i[HIGH_HALF]=n1;
00094       z=u.x*r.x;
00095       d=(z+big.x)-big.x;
00096       u.x=(u.x-d*w.x)-d*ww.x;
00097       if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
00098       else
00099         if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
00100         else
00101         {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
00102     }
00103 
00104   }   /*   (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000)     */
00105   else {
00106     if (kx<0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
00107       y=ABS(y)*t128.x;
00108       z=__ieee754_remainder(x,y)*t128.x;
00109       z=__ieee754_remainder(z,y)*tm128.x;
00110       return z;
00111     }
00112   else {
00113     if ((kx&0x7ff00000)==0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
00114       y=ABS(y);
00115       z=2.0*__ieee754_remainder(0.5*x,y);
00116       d = ABS(z);
00117       if (d <= ABS(d-y)) return z;
00118       else return (z>0)?z-y:z+y;
00119     }
00120     else { /* if x is too big */
00121       if (kx == 0x7ff00000 && u.i[LOW_HALF] == 0 && y == 1.0)
00122        return x / x;
00123       if (kx>=0x7ff00000||(ky==0&&t.i[LOW_HALF]==0)||ky>0x7ff00000||
00124          (ky==0x7ff00000&&t.i[LOW_HALF]!=0))
00125        return (u.i[HIGH_HALF]&0x80000000)?nNAN.x:NAN.x;
00126       else return x;
00127     }
00128    }
00129   }
00130 }