lib/fmodl.c

   1 /* Remainder.
   2    Copyright (C) 2011-2012 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 #include <config.h>
  18
  19 /* Specification.  */
  20 #include <math.h>
  21
  22 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
  23
  24 long double
  25 fmodl (long double x, long double y)
  26 {
  27   return fmod (x, y);
  28 }
  29
  30 #else
  31
  32 long double
  33 fmodl (long double x, long double y)
  34 {
  35   if (isinf (y))
  36     return x;
  37   else
  38     {
  39       long double q = - truncl (x / y);
  40       long double r = fmal (q, y, x); /* = x + q * y, computed in one step */
  41       /* Correct possible rounding errors in the quotient x / y.  */
  42       if (y >= 0)
  43         {
  44           if (x >= 0)
  45             {
  46               /* Expect 0 <= r < y.  */
  47               if (r < 0)
  48                 q += 1, r = fmal (q, y, x);
  49               else if (r >= y)
  50                 q -= 1, r = fmal (q, y, x);
  51             }
  52           else
  53             {
  54               /* Expect - y < r <= 0.  */
  55               if (r > 0)
  56                 q -= 1, r = fmal (q, y, x);
  57               else if (r <= - y)
  58                 q += 1, r = fmal (q, y, x);
  59             }
  60         }
  61       else
  62         {
  63           if (x >= 0)
  64             {
  65               /* Expect 0 <= r < - y.  */
  66               if (r < 0)
  67                 q -= 1, r = fmal (q, y, x);
  68               else if (r >= - y)
  69                 q += 1, r = fmal (q, y, x);
  70             }
  71           else
  72             {
  73               /* Expect y < r <= 0.  */
  74               if (r > 0)
  75                 q += 1, r = fmal (q, y, x);
  76               else if (r <= y)
  77                 q -= 1, r = fmal (q, y, x);
  78             }
  79         }
  80       return r;
  81     }
  82 }
  83
  84 #endif