lib/remainder.c

   1 /* Remainder.
   2    Copyright (C) 2012-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 || defined USE_FLOAT)
  18 # include <config.h>
  19 #endif
  20
  21 /* Specification.  */
  22 #include <math.h>
  23
  24 #ifdef USE_LONG_DOUBLE
  25 # define REMAINDER remainderl
  26 # define DOUBLE long double
  27 # define L_(literal) literal##L
  28 # define FABS fabsl
  29 # define FMOD fmodl
  30 # define ISNAN isnanl
  31 #elif ! defined USE_FLOAT
  32 # define REMAINDER remainder
  33 # define DOUBLE double
  34 # define L_(literal) literal
  35 # define FABS fabs
  36 # define FMOD fmod
  37 # define ISNAN isnand
  38 #else /* defined USE_FLOAT */
  39 # define REMAINDER remainderf
  40 # define DOUBLE float
  41 # define L_(literal) literal##f
  42 # define FABS fabsf
  43 # define FMOD fmodf
  44 # define ISNAN isnanf
  45 #endif
  46
  47 #undef NAN
  48 #if defined _MSC_VER
  49 static DOUBLE zero;
  50 # define NAN (zero / zero)
  51 #else
  52 # define NAN (L_(0.0) / L_(0.0))
  53 #endif
  54
  55 DOUBLE
  56 REMAINDER (DOUBLE x, DOUBLE y)
  57 {
  58   if (isfinite (x) && isfinite (y) && y != L_(0.0))
  59     {
  60       if (x == L_(0.0))
  61         /* Return x, regardless of the sign of y.  */
  62         return x;
  63
  64       {
  65         int negate = ((!signbit (x)) ^ (!signbit (y)));
  66         DOUBLE r;
  67
  68         /* Take the absolute value of x and y.  */
  69         x = FABS (x);
  70         y = FABS (y);
  71
  72         /* Trivial case that requires no computation.  */
  73         if (x <= L_(0.5) * y)
  74           return (negate ? - x : x);
  75
  76         /* With a fixed y, the function x -> remainder(x,y) has a period 2*y.
  77            Therefore we can reduce the argument x modulo 2*y.  And it's no
  78            problem if 2*y overflows, since fmod(x,Inf) = x.  */
  79         x = FMOD (x, L_(2.0) * y);
  80
  81         /* Consider the 3 cases:
  82              0 <= x <= 0.5 * y
  83              0.5 * y < x < 1.5 * y
  84              1.5 * y <= x <= 2.0 * y  */
  85         if (x <= L_(0.5) * y)
  86           r = x;
  87         else
  88           {
  89             r = x - y;
  90             if (r > L_(0.5) * y)
  91               r = x - L_(2.0) * y;
  92           }
  93         return (negate ? - r : r);
  94       }
  95     }
  96   else
  97     {
  98       if (ISNAN (x) || ISNAN (y))
  99         return x + y; /* NaN */
 100       else if (isinf (y))
 101         return x;
 102       else
 103         /* x infinite or y zero */
 104         return NAN;
 105     }
 106 }