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   long double q = - truncl (x / y);
  36   long double r = fmal (q, y, x); /* = x + q * y, computed in one step */
  37   /* Correct possible rounding errors in the quotient x / y.  */
  38   if (y >= 0)
  39     {
  40       if (x >= 0)
  41         {
  42           /* Expect 0 <= r < y.  */
  43           if (r < 0)
  44             q += 1, r = fmal (q, y, x);
  45           else if (r >= y)
  46             q -= 1, r = fmal (q, y, x);
  47         }
  48       else
  49         {
  50           /* Expect - y < r <= 0.  */
  51           if (r > 0)
  52             q -= 1, r = fmal (q, y, x);
  53           else if (r <= - y)
  54             q += 1, r = fmal (q, y, x);
  55         }
  56     }
  57   else
  58     {
  59       if (x >= 0)
  60         {
  61           /* Expect 0 <= r < - y.  */
  62           if (r < 0)
  63             q -= 1, r = fmal (q, y, x);
  64           else if (r >= - y)
  65             q += 1, r = fmal (q, y, x);
  66         }
  67       else
  68         {
  69           /* Expect y < r <= 0.  */
  70           if (r > 0)
  71             q += 1, r = fmal (q, y, x);
  72           else if (r <= y)
  73             q -= 1, r = fmal (q, y, x);
  74         }
  75     }
  76   return r;
  77 }
  78
  79 #endif