sysdeps/ieee754/dbl-64/e_remainder.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2013 Free Software Foundation, Inc.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  18  */
  19 /**************************************************************************/
  20 /*  MODULE_NAME urem.c                                                    */
  21 /*                                                                        */
  22 /*  FUNCTION: uremainder                                                  */
  23 /*                                                                        */
  24 /* An ultimate remainder routine. Given two IEEE double machine numbers x */
  25 /* ,y   it computes the correctly rounded (to nearest) value of remainder */
  26 /* of dividing x by y.                                                    */
  27 /* Assumption: Machine arithmetic operations are performed in             */
  28 /* round to nearest mode of IEEE 754 standard.                            */
  29 /*                                                                        */
  30 /* ************************************************************************/
  31
  32 #include "endian.h"
  33 #include "mydefs.h"
  34 #include "urem.h"
  35 #include "MathLib.h"
  36 #include <math_private.h>
  37
  38 /**************************************************************************/
  39 /* An ultimate remainder routine. Given two IEEE double machine numbers x */
  40 /* ,y   it computes the correctly rounded (to nearest) value of remainder */
  41 /**************************************************************************/
  42 double __ieee754_remainder(double x, double y)
  43 {
  44   double z,d,xx;
  45   int4 kx,ky,n,nn,n1,m1,l;
  46   mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r;
  47   u.x=x;
  48   t.x=y;
  49   kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign  for x*/
  50   t.i[HIGH_HALF]&=0x7fffffff;   /*no sign for y */
  51   ky=t.i[HIGH_HALF];
  52   /*------ |x| < 2^1023  and   2^-970 < |y| < 2^1024 ------------------*/
  53   if (kx<0x7fe00000 && ky<0x7ff00000 && ky>=0x03500000) {
  54     if (kx+0x00100000<ky) return x;
  55     if ((kx-0x01500000)<ky) {
  56       z=x/t.x;
  57       v.i[HIGH_HALF]=t.i[HIGH_HALF];
  58       d=(z+big.x)-big.x;
  59       xx=(x-d*v.x)-d*(t.x-v.x);
  60       if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x);
  61       else {
  62         if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x;
  63         else return xx;
  64       }
  65     }   /*    (kx<(ky+0x01500000))         */
  66     else  {
  67       r.x=1.0/t.x;
  68       n=t.i[HIGH_HALF];
  69       nn=(n&0x7ff00000)+0x01400000;
  70       w.i[HIGH_HALF]=n;
  71       ww.x=t.x-w.x;
  72       l=(kx-nn)&0xfff00000;
  73       n1=ww.i[HIGH_HALF];
  74       m1=r.i[HIGH_HALF];
  75       while (l>0) {
  76         r.i[HIGH_HALF]=m1-l;
  77         z=u.x*r.x;
  78         w.i[HIGH_HALF]=n+l;
  79         ww.i[HIGH_HALF]=(n1)?n1+l:n1;
  80         d=(z+big.x)-big.x;
  81         u.x=(u.x-d*w.x)-d*ww.x;
  82         l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
  83       }
  84       r.i[HIGH_HALF]=m1;
  85       w.i[HIGH_HALF]=n;
  86       ww.i[HIGH_HALF]=n1;
  87       z=u.x*r.x;
  88       d=(z+big.x)-big.x;
  89       u.x=(u.x-d*w.x)-d*ww.x;
  90       if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
  91       else
  92         if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
  93         else
  94         {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
  95     }
  96
  97   }   /*   (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000)     */
  98   else {
  99     if (kx<0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
 100       y=ABS(y)*t128.x;
 101       z=__ieee754_remainder(x,y)*t128.x;
 102       z=__ieee754_remainder(z,y)*tm128.x;
 103       return z;
 104     }
 105   else {
 106     if ((kx&0x7ff00000)==0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
 107       y=ABS(y);
 108       z=2.0*__ieee754_remainder(0.5*x,y);
 109       d = ABS(z);
 110       if (d <= ABS(d-y)) return z;
 111       else return (z>0)?z-y:z+y;
 112     }
 113     else { /* if x is too big */
 114       if (ky==0 && t.i[LOW_HALF] == 0)          /* y = 0 */
 115         return (x * y) / (x * y);
 116       else if (kx >= 0x7ff00000                 /* x not finite */
 117                || (ky>0x7ff00000                /* y is NaN */
 118                    || (ky == 0x7ff00000 && t.i[LOW_HALF] != 0)))
 119         return (x * y) / (x * y);
 120       else
 121         return x;
 122     }
 123    }
 124   }
 125 }
 126 strong_alias (__ieee754_remainder, __remainder_finite)