lib/remainderl.c

   1 /* Remainder.
   2    Copyright (C) 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 remainderl (long double x, long double y)
  26 {
  27   return remainder (x, y);
  28 }
  29
  30 #else
  31
  32 long double
  33 remainderl (long double x, long double y)
  34 {
  35   long double q = - roundl (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   long double half_y = 0.5L * y;
  39   if (y >= 0)
  40     {
  41       /* Expect -y/2 <= r <= y/2.  */
  42       if (r > half_y)
  43         q -= 1, r = fmal (q, y, x);
  44       else if (r < - half_y)
  45         q += 1, r = fmal (q, y, x);
  46     }
  47   else
  48     {
  49       /* Expect y/2 <= r <= -y/2.  */
  50       if (r > - half_y)
  51         q += 1, r = fmal (q, y, x);
  52       else if (r < half_y)
  53         q -= 1, r = fmal (q, y, x);
  54     }
  55   return r;
  56 }
  57
  58 #endif