lib/fmod.c

   1 /* Remainder.
   2    Copyright (C) 2011-2014 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
  16
  17 #if ! defined USE_LONG_DOUBLE
  18 # include <config.h>
  19 #endif
  20
  21 /* Specification.  */
  22 #include <math.h>
  23
  24 #include <float.h>
  25 #include <stdlib.h>
  26
  27 #ifdef USE_LONG_DOUBLE
  28 # define FMOD fmodl
  29 # define DOUBLE long double
  30 # define MANT_DIG LDBL_MANT_DIG
  31 # define L_(literal) literal##L
  32 # define FREXP frexpl
  33 # define LDEXP ldexpl
  34 # define FABS fabsl
  35 # define TRUNC truncl
  36 # define ISNAN isnanl
  37 #else
  38 # define FMOD fmod
  39 # define DOUBLE double
  40 # define MANT_DIG DBL_MANT_DIG
  41 # define L_(literal) literal
  42 # define FREXP frexp
  43 # define LDEXP ldexp
  44 # define FABS fabs
  45 # define TRUNC trunc
  46 # define ISNAN isnand
  47 #endif
  48
  49 /* MSVC with option -fp:strict refuses to compile constant initializers that
  50    contain floating-point operations.  Pacify this compiler.  */
  51 #ifdef _MSC_VER
  52 # pragma fenv_access (off)
  53 #endif
  54
  55 #undef NAN
  56 #if defined _MSC_VER
  57 static DOUBLE zero;
  58 # define NAN (zero / zero)
  59 #else
  60 # define NAN (L_(0.0) / L_(0.0))
  61 #endif
  62
  63 /* To avoid rounding errors during the computation of x - q * y,
  64    there are three possibilities:
  65      - Use fma (- q, y, x).
  66      - Have q be a single bit at a time, and compute x - q * y
  67        through a subtraction.
  68      - Have q be at most MANT_DIG/2 bits at a time, and compute
  69        x - q * y by splitting y into two halves:
  70        y = y1 * 2^(MANT_DIG/2) + y0
  71        x - q * y = (x - q * y1 * 2^MANT_DIG/2) - q * y0.
  72    The latter is probably the most efficient.  */
  73
  74 /* Number of bits in a limb.  */
  75 #define LIMB_BITS (MANT_DIG / 2)
  76
  77 /* 2^LIMB_BITS.  */
  78 static const DOUBLE TWO_LIMB_BITS =
  79   /* Assume LIMB_BITS <= 3 * 31.
  80      Use the identity
  81        n = floor(n/3) + floor((n+1)/3) + floor((n+2)/3).  */
  82   (DOUBLE) (1U << (LIMB_BITS / 3))
  83   * (DOUBLE) (1U << ((LIMB_BITS + 1) / 3))
  84   * (DOUBLE) (1U << ((LIMB_BITS + 2) / 3));
  85
  86 /* 2^-LIMB_BITS.  */
  87 static const DOUBLE TWO_LIMB_BITS_INVERSE =
  88   /* Assume LIMB_BITS <= 3 * 31.
  89      Use the identity
  90        n = floor(n/3) + floor((n+1)/3) + floor((n+2)/3).  */
  91   L_(1.0) / ((DOUBLE) (1U << (LIMB_BITS / 3))
  92              * (DOUBLE) (1U << ((LIMB_BITS + 1) / 3))
  93              * (DOUBLE) (1U << ((LIMB_BITS + 2) / 3)));
  94
  95 DOUBLE
  96 FMOD (DOUBLE x, DOUBLE y)
  97 {
  98   if (isfinite (x) && isfinite (y) && y != L_(0.0))
  99     {
 100       if (x == L_(0.0))
 101         /* Return x, regardless of the sign of y.  */
 102         return x;
 103
 104       {
 105         int negate = ((!signbit (x)) ^ (!signbit (y)));
 106
 107         /* Take the absolute value of x and y.  */
 108         x = FABS (x);
 109         y = FABS (y);
 110
 111         /* Trivial case that requires no computation.  */
 112         if (x < y)
 113           return (negate ? - x : x);
 114
 115         {
 116           int yexp;
 117           DOUBLE ym;
 118           DOUBLE y1;
 119           DOUBLE y0;
 120           int k;
 121           DOUBLE x2;
 122           DOUBLE x1;
 123           DOUBLE x0;
 124
 125           /* Write y = 2^yexp * (y1 * 2^-LIMB_BITS + y0 * 2^-(2*LIMB_BITS))
 126              where y1 is an integer, 2^(LIMB_BITS-1) <= y1 < 2^LIMB_BITS,
 127              y1 has at most LIMB_BITS bits,
 128              0 <= y0 < 2^LIMB_BITS,
 129              y0 has at most (MANT_DIG + 1) / 2 bits.  */
 130           ym = FREXP (y, &yexp);
 131           ym = ym * TWO_LIMB_BITS;
 132           y1 = TRUNC (ym);
 133           y0 = (ym - y1) * TWO_LIMB_BITS;
 134
 135           /* Write
 136                x = 2^(yexp+(k-3)*LIMB_BITS)
 137                    * (x2 * 2^(2*LIMB_BITS) + x1 * 2^LIMB_BITS + x0)
 138              where x2, x1, x0 are each integers >= 0, < 2^LIMB_BITS.  */
 139           {
 140             int xexp;
 141             DOUBLE xm = FREXP (x, &xexp);
 142             /* Since we know x >= y, we know xexp >= yexp.  */
 143             xexp -= yexp;
 144             /* Compute k = ceil(xexp / LIMB_BITS).  */
 145             k = (xexp + LIMB_BITS - 1) / LIMB_BITS;
 146             /* Note: (k - 1) * LIMB_BITS + 1 <= xexp <= k * LIMB_BITS.  */
 147             /* Note: 0.5 <= xm < 1.0.  */
 148             xm = LDEXP (xm, xexp - (k - 1) * LIMB_BITS);
 149             /* Note: Now xm < 2^(xexp - (k - 1) * LIMB_BITS) <= 2^LIMB_BITS
 150                and xm >= 0.5 * 2^(xexp - (k - 1) * LIMB_BITS) >= 1.0
 151                and xm has at most MANT_DIG <= 2*LIMB_BITS+1 bits.  */
 152             x2 = TRUNC (xm);
 153             x1 = (xm - x2) * TWO_LIMB_BITS;
 154             /* Split off x0 from x1 later.  */
 155           }
 156
 157           /* Test whether [x2,x1,0] >= 2^LIMB_BITS * [y1,y0].  */
 158           if (x2 > y1 || (x2 == y1 && x1 >= y0))
 159             {
 160               /* Subtract 2^LIMB_BITS * [y1,y0] from [x2,x1,0].  */
 161               x2 -= y1;
 162               x1 -= y0;
 163               if (x1 < L_(0.0))
 164                 {
 165                   if (!(x2 >= L_(1.0)))
 166                     abort ();
 167                   x2 -= L_(1.0);
 168                   x1 += TWO_LIMB_BITS;
 169                 }
 170             }
 171
 172           /* Split off x0 from x1.  */
 173           {
 174             DOUBLE x1int = TRUNC (x1);
 175             x0 = TRUNC ((x1 - x1int) * TWO_LIMB_BITS);
 176             x1 = x1int;
 177           }
 178
 179           for (; k > 0; k--)
 180             {
 181               /* Multiprecision division of the limb sequence [x2,x1,x0]
 182                  by [y1,y0].  */
 183               /* Here [x2,x1,x0] < 2^LIMB_BITS * [y1,y0].  */
 184               /* The first guess takes into account only [x2,x1] and [y1].
 185
 186                  By Knuth's theorem, we know that
 187                    q* = min (floor ([x2,x1] / [y1]), 2^LIMB_BITS - 1)
 188                  and
 189                    q = floor ([x2,x1,x0] / [y1,y0])
 190                  are not far away:
 191                    q* - 2 <= q <= q* + 1.
 192
 193                  Proof:
 194                  a) q* * y1 <= floor ([x2,x1] / [y1]) * y1 <= [x2,x1].
 195                     Hence
 196                     [x2,x1,x0] - q* * [y1,y0]
 197                       = 2^LIMB_BITS * ([x2,x1] - q* * [y1]) + x0 - q* * y0
 198                       >= x0 - q* * y0
 199                       >= - q* * y0
 200                       > - 2^(2*LIMB_BITS)
 201                       >= - 2 * [y1,y0]
 202                     So
 203                       [x2,x1,x0] > (q* - 2) * [y1,y0].
 204                  b) If q* = floor ([x2,x1] / [y1]), then
 205                       [x2,x1] < (q* + 1) * y1
 206                     Hence
 207                     [x2,x1,x0] - q* * [y1,y0]
 208                       = 2^LIMB_BITS * ([x2,x1] - q* * [y1]) + x0 - q* * y0
 209                       <= 2^LIMB_BITS * (y1 - 1) + x0 - q* * y0
 210                       <= 2^LIMB_BITS * (2^LIMB_BITS-2) + (2^LIMB_BITS-1) - 0
 211                       < 2^(2*LIMB_BITS)
 212                       <= 2 * [y1,y0]
 213                     So
 214                       [x2,x1,x0] < (q* + 2) * [y1,y0].
 215                     and so
 216                       q < q* + 2
 217                     which implies
 218                       q <= q* + 1.
 219                     In the other case, q* = 2^LIMB_BITS - 1.  Then trivially
 220                       q < 2^LIMB_BITS = q* + 1.
 221
 222                  We know that floor ([x2,x1] / [y1]) >= 2^LIMB_BITS if and
 223                  only if x2 >= y1.  */
 224               DOUBLE q =
 225                 (x2 >= y1
 226                  ? TWO_LIMB_BITS - L_(1.0)
 227                  : TRUNC ((x2 * TWO_LIMB_BITS + x1) / y1));
 228               if (q > L_(0.0))
 229                 {
 230                   /* Compute
 231                      [x2,x1,x0] - q* * [y1,y0]
 232                        = 2^LIMB_BITS * ([x2,x1] - q* * [y1]) + x0 - q* * y0.  */
 233                   DOUBLE q_y1 = q * y1; /* exact, at most 2*LIMB_BITS bits */
 234                   DOUBLE q_y1_1 = TRUNC (q_y1 * TWO_LIMB_BITS_INVERSE);
 235                   DOUBLE q_y1_0 = q_y1 - q_y1_1 * TWO_LIMB_BITS;
 236                   DOUBLE q_y0 = q * y0; /* exact, at most MANT_DIG bits */
 237                   DOUBLE q_y0_1 = TRUNC (q_y0 * TWO_LIMB_BITS_INVERSE);
 238                   DOUBLE q_y0_0 = q_y0 - q_y0_1 * TWO_LIMB_BITS;
 239                   x2 -= q_y1_1;
 240                   x1 -= q_y1_0;
 241                   x1 -= q_y0_1;
 242                   x0 -= q_y0_0;
 243                   /* Move negative carry from x0 to x1 and from x1 to x2.  */
 244                   if (x0 < L_(0.0))
 245                     {
 246                       x0 += TWO_LIMB_BITS;
 247                       x1 -= L_(1.0);
 248                     }
 249                   if (x1 < L_(0.0))
 250                     {
 251                       x1 += TWO_LIMB_BITS;
 252                       x2 -= L_(1.0);
 253                       if (x1 < L_(0.0)) /* not sure this can happen */
 254                         {
 255                           x1 += TWO_LIMB_BITS;
 256                           x2 -= L_(1.0);
 257                         }
 258                     }
 259                   if (x2 < L_(0.0))
 260                     {
 261                       /* Reduce q by 1.  */
 262                       x1 += y1;
 263                       x0 += y0;
 264                       /* Move overflow from x0 to x1 and from x1 to x0.  */
 265                       if (x0 >= TWO_LIMB_BITS)
 266                         {
 267                           x0 -= TWO_LIMB_BITS;
 268                           x1 += L_(1.0);
 269                         }
 270                       if (x1 >= TWO_LIMB_BITS)
 271                         {
 272                           x1 -= TWO_LIMB_BITS;
 273                           x2 += L_(1.0);
 274                         }
 275                       if (x2 < L_(0.0))
 276                         {
 277                           /* Reduce q by 1 again.  */
 278                           x1 += y1;
 279                           x0 += y0;
 280                           /* Move overflow from x0 to x1 and from x1 to x0.  */
 281                           if (x0 >= TWO_LIMB_BITS)
 282                             {
 283                               x0 -= TWO_LIMB_BITS;
 284                               x1 += L_(1.0);
 285                             }
 286                           if (x1 >= TWO_LIMB_BITS)
 287                             {
 288                               x1 -= TWO_LIMB_BITS;
 289                               x2 += L_(1.0);
 290                             }
 291                           if (x2 < L_(0.0))
 292                             /* Shouldn't happen, because we proved that
 293                                q >= q* - 2.  */
 294                             abort ();
 295                         }
 296                     }
 297                 }
 298               if (x2 > L_(0.0)
 299                   || x1 > y1
 300                   || (x1 == y1 && x0 >= y0))
 301                 {
 302                   /* Increase q by 1.  */
 303                   x1 -= y1;
 304                   x0 -= y0;
 305                   /* Move negative carry from x0 to x1 and from x1 to x2.  */
 306                   if (x0 < L_(0.0))
 307                     {
 308                       x0 += TWO_LIMB_BITS;
 309                       x1 -= L_(1.0);
 310                     }
 311                   if (x1 < L_(0.0))
 312                     {
 313                       x1 += TWO_LIMB_BITS;
 314                       x2 -= L_(1.0);
 315                     }
 316                   if (x2 < L_(0.0))
 317                     abort ();
 318                   if (x2 > L_(0.0)
 319                       || x1 > y1
 320                       || (x1 == y1 && x0 >= y0))
 321                     /* Shouldn't happen, because we proved that
 322                        q <= q* + 1.  */
 323                     abort ();
 324                 }
 325               /* Here [x2,x1,x0] < [y1,y0].  */
 326               /* Next round.  */
 327               x2 = x1;
 328 #if (MANT_DIG + 1) / 2 > LIMB_BITS /* y0 can have a fractional bit */
 329               x1 = TRUNC (x0);
 330               x0 = (x0 - x1) * TWO_LIMB_BITS;
 331 #else
 332               x1 = x0;
 333               x0 = L_(0.0);
 334 #endif
 335               /* Here [x2,x1,x0] < 2^LIMB_BITS * [y1,y0].  */
 336             }
 337           /* Here k = 0.
 338              The result is
 339                2^(yexp-3*LIMB_BITS)
 340                * (x2 * 2^(2*LIMB_BITS) + x1 * 2^LIMB_BITS + x0).  */
 341           {
 342             DOUBLE r =
 343               LDEXP ((x2 * TWO_LIMB_BITS + x1) * TWO_LIMB_BITS + x0,
 344                      yexp - 3 * LIMB_BITS);
 345             return (negate ? - r : r);
 346           }
 347         }
 348       }
 349     }
 350   else
 351     {
 352       if (ISNAN (x) || ISNAN (y))
 353         return x + y; /* NaN */
 354       else if (isinf (y))
 355         return x;
 356       else
 357         /* x infinite or y zero */
 358         return NAN;
 359     }
 360 }