sysdeps/ieee754/dbl-64/e_remainder.c

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