sysdeps/ieee754/dbl-64/e_remainder.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2014 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
  43 __ieee754_remainder (double x, double y)
  44 {
  45   double z, d, xx;
  46   int4 kx, ky, n, nn, n1, m1, l;
  47   mynumber u, t, w = { { 0, 0 } }, v = { { 0, 0 } }, ww = { { 0, 0 } }, r;
  48   u.x = x;
  49   t.x = y;
  50   kx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign  for x*/
  51   t.i[HIGH_HALF] &= 0x7fffffff;   /*no sign for y */
  52   ky = t.i[HIGH_HALF];
  53   /*------ |x| < 2^1023  and   2^-970 < |y| < 2^1024 ------------------*/
  54   if (kx < 0x7fe00000 && ky < 0x7ff00000 && ky >= 0x03500000)
  55     {
  56       SET_RESTORE_ROUND_NOEX (FE_TONEAREST);
  57       if (kx + 0x00100000 < ky)
  58         return x;
  59       if ((kx - 0x01500000) < ky)
  60         {
  61           z = x / t.x;
  62           v.i[HIGH_HALF] = t.i[HIGH_HALF];
  63           d = (z + big.x) - big.x;
  64           xx = (x - d * v.x) - d * (t.x - v.x);
  65           if (d - z != 0.5 && d - z != -0.5)
  66             return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x);
  67           else
  68             {
  69               if (ABS (xx) > 0.5 * t.x)
  70                 return (z > d) ? xx - t.x : xx + t.x;
  71               else
  72                 return xx;
  73             }
  74         } /*    (kx<(ky+0x01500000))         */
  75       else
  76         {
  77           r.x = 1.0 / t.x;
  78           n = t.i[HIGH_HALF];
  79           nn = (n & 0x7ff00000) + 0x01400000;
  80           w.i[HIGH_HALF] = n;
  81           ww.x = t.x - w.x;
  82           l = (kx - nn) & 0xfff00000;
  83           n1 = ww.i[HIGH_HALF];
  84           m1 = r.i[HIGH_HALF];
  85           while (l > 0)
  86             {
  87               r.i[HIGH_HALF] = m1 - l;
  88               z = u.x * r.x;
  89               w.i[HIGH_HALF] = n + l;
  90               ww.i[HIGH_HALF] = (n1) ? n1 + l : n1;
  91               d = (z + big.x) - big.x;
  92               u.x = (u.x - d * w.x) - d * ww.x;
  93               l = (u.i[HIGH_HALF] & 0x7ff00000) - nn;
  94             }
  95           r.i[HIGH_HALF] = m1;
  96           w.i[HIGH_HALF] = n;
  97           ww.i[HIGH_HALF] = n1;
  98           z = u.x * r.x;
  99           d = (z + big.x) - big.x;
 100           u.x = (u.x - d * w.x) - d * ww.x;
 101           if (ABS (u.x) < 0.5 * t.x)
 102             return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x);
 103           else
 104           if (ABS (u.x) > 0.5 * t.x)
 105             return (d > z) ? u.x + t.x : u.x - t.x;
 106           else
 107             {
 108               z = u.x / t.x; d = (z + big.x) - big.x;
 109               return ((u.x - d * w.x) - d * ww.x);
 110             }
 111         }
 112     } /*   (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000)     */
 113   else
 114     {
 115       if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0))
 116         {
 117           y = ABS (y) * t128.x;
 118           z = __ieee754_remainder (x, y) * t128.x;
 119           z = __ieee754_remainder (z, y) * tm128.x;
 120           return z;
 121         }
 122       else
 123         {
 124           if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 &&
 125               (ky > 0 || t.i[LOW_HALF] != 0))
 126             {
 127               y = ABS (y);
 128               z = 2.0 * __ieee754_remainder (0.5 * x, y);
 129               d = ABS (z);
 130               if (d <= ABS (d - y))
 131                 return z;
 132               else
 133                 return (z > 0) ? z - y : z + y;
 134             }
 135           else /* if x is too big */
 136             {
 137               if (ky == 0 && t.i[LOW_HALF] == 0) /* y = 0 */
 138                 return (x * y) / (x * y);
 139               else if (kx >= 0x7ff00000         /* x not finite */
 140                        || (ky > 0x7ff00000      /* y is NaN */
 141                            || (ky == 0x7ff00000 && t.i[LOW_HALF] != 0)))
 142                 return (x * y) / (x * y);
 143               else
 144                 return x;
 145             }
 146         }
 147     }
 148 }
 149 strong_alias (__ieee754_remainder, __remainder_finite)