contrib/mpfr/src/rem1.c

   1 /* mpfr_rem1 -- internal function
   2    mpfr_fmod -- compute the floating-point remainder of x/y
   3    mpfr_remquo and mpfr_remainder -- argument reduction functions
   4
   5 Copyright 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   6 Contributed by the AriC and Caramel projects, INRIA.
   7
   8 This file is part of the GNU MPFR Library.
   9
  10 The GNU MPFR Library is free software; you can redistribute it and/or modify
  11 it under the terms of the GNU Lesser General Public License as published by
  12 the Free Software Foundation; either version 3 of the License, or (at your
  13 option) any later version.
  14
  15 The GNU MPFR Library is distributed in the hope that it will be useful, but
  16 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  17 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  18 License for more details.
  19
  20 You should have received a copy of the GNU Lesser General Public License
  21 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  22 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  23 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  24
  25 # include "mpfr-impl.h"
  26
  27 /* we return as many bits as we can, keeping just one bit for the sign */
  28 # define WANTED_BITS (sizeof(long) * CHAR_BIT - 1)
  29
  30 /*
  31   rem1 works as follows:
  32   The first rounding mode rnd_q indicate if we are actually computing
  33   a fmod (MPFR_RNDZ) or a remainder/remquo (MPFR_RNDN).
  34
  35   Let q = x/y rounded to an integer in the direction rnd_q.
  36   Put x - q*y in rem, rounded according to rnd.
  37   If quo is not null, the value stored in *quo has the sign of q,
  38   and agrees with q with the 2^n low order bits.
  39   In other words, *quo = q (mod 2^n) and *quo q >= 0.
  40   If rem is zero, then it has the sign of x.
  41   The returned 'int' is the inexact flag giving the place of rem wrt x - q*y.
  42
  43   If x or y is NaN: *quo is undefined, rem is NaN.
  44   If x is Inf, whatever y: *quo is undefined, rem is NaN.
  45   If y is Inf, x not NaN nor Inf: *quo is 0, rem is x.
  46   If y is 0, whatever x: *quo is undefined, rem is NaN.
  47   If x is 0, whatever y (not NaN nor 0): *quo is 0, rem is x.
  48
  49   Otherwise if x and y are neither NaN, Inf nor 0, q is always defined,
  50   thus *quo is.
  51   Since |x - q*y| <= y/2, no overflow is possible.
  52   Only an underflow is possible when y is very small.
  53  */
  54
  55 static int
  56 mpfr_rem1 (mpfr_ptr rem, long *quo, mpfr_rnd_t rnd_q,
  57            mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
  58 {
  59   mpfr_exp_t ex, ey;
  60   int compare, inex, q_is_odd, sign, signx = MPFR_SIGN (x);
  61   mpz_t mx, my, r;
  62
  63   MPFR_ASSERTD (rnd_q == MPFR_RNDN || rnd_q == MPFR_RNDZ);
  64
  65   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y)))
  66     {
  67       if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x)
  68           || MPFR_IS_ZERO (y))
  69         {
  70           /* for remquo, quo is undefined */
  71           MPFR_SET_NAN (rem);
  72           MPFR_RET_NAN;
  73         }
  74       else                      /* either y is Inf and x is 0 or non-special,
  75                                    or x is 0 and y is non-special,
  76                                    in both cases the quotient is zero. */
  77         {
  78           if (quo)
  79             *quo = 0;
  80           return mpfr_set (rem, x, rnd);
  81         }
  82     }
  83
  84   /* now neither x nor y is NaN, Inf or zero */
  85
  86   mpz_init (mx);
  87   mpz_init (my);
  88   mpz_init (r);
  89
  90   ex = mpfr_get_z_2exp (mx, x);  /* x = mx*2^ex */
  91   ey = mpfr_get_z_2exp (my, y);  /* y = my*2^ey */
  92
  93   /* to get rid of sign problems, we compute it separately:
  94      quo(-x,-y) = quo(x,y), rem(-x,-y) = -rem(x,y)
  95      quo(-x,y) = -quo(x,y), rem(-x,y)  = -rem(x,y)
  96      thus quo = sign(x/y)*quo(|x|,|y|), rem = sign(x)*rem(|x|,|y|) */
  97   sign = (signx == MPFR_SIGN (y)) ? 1 : -1;
  98   mpz_abs (mx, mx);
  99   mpz_abs (my, my);
 100   q_is_odd = 0;
 101
 102   /* divide my by 2^k if possible to make operations mod my easier */
 103   {
 104     unsigned long k = mpz_scan1 (my, 0);
 105     ey += k;
 106     mpz_fdiv_q_2exp (my, my, k);
 107   }
 108
 109   if (ex <= ey)
 110     {
 111       /* q = x/y = mx/(my*2^(ey-ex)) */
 112       mpz_mul_2exp (my, my, ey - ex);   /* divide mx by my*2^(ey-ex) */
 113       if (rnd_q == MPFR_RNDZ)
 114         /* 0 <= |r| <= |my|, r has the same sign as mx */
 115         mpz_tdiv_qr (mx, r, mx, my);
 116       else
 117         /* 0 <= |r| <= |my|, r has the same sign as my */
 118         mpz_fdiv_qr (mx, r, mx, my);
 119
 120       if (rnd_q == MPFR_RNDN)
 121         q_is_odd = mpz_tstbit (mx, 0);
 122       if (quo)                  /* mx is the quotient */
 123         {
 124           mpz_tdiv_r_2exp (mx, mx, WANTED_BITS);
 125           *quo = mpz_get_si (mx);
 126         }
 127     }
 128   else                          /* ex > ey */
 129     {
 130       if (quo) /* remquo case */
 131         /* for remquo, to get the low WANTED_BITS more bits of the quotient,
 132            we first compute R =  X mod Y*2^WANTED_BITS, where X and Y are
 133            defined below. Then the low WANTED_BITS of the quotient are
 134            floor(R/Y). */
 135         mpz_mul_2exp (my, my, WANTED_BITS);     /* 2^WANTED_BITS*Y */
 136
 137       else if (rnd_q == MPFR_RNDN) /* remainder case */
 138         /* Let X = mx*2^(ex-ey) and Y = my. Then both X and Y are integers.
 139            Assume X = R mod Y, then x = X*2^ey = R*2^ey mod (Y*2^ey=y).
 140            To be able to perform the rounding, we need the least significant
 141            bit of the quotient, i.e., one more bit in the remainder,
 142            which is obtained by dividing by 2Y. */
 143         mpz_mul_2exp (my, my, 1);       /* 2Y */
 144
 145       mpz_set_ui (r, 2);
 146       mpz_powm_ui (r, r, ex - ey, my);  /* 2^(ex-ey) mod my */
 147       mpz_mul (r, r, mx);
 148       mpz_mod (r, r, my);
 149
 150       if (quo)                  /* now 0 <= r < 2^WANTED_BITS*Y */
 151         {
 152           mpz_fdiv_q_2exp (my, my, WANTED_BITS);   /* back to Y */
 153           mpz_tdiv_qr (mx, r, r, my);
 154           /* oldr = mx*my + newr */
 155           *quo = mpz_get_si (mx);
 156           q_is_odd = *quo & 1;
 157         }
 158       else if (rnd_q == MPFR_RNDN) /* now 0 <= r < 2Y in the remainder case */
 159         {
 160           mpz_fdiv_q_2exp (my, my, 1);     /* back to Y */
 161           /* least significant bit of q */
 162           q_is_odd = mpz_cmpabs (r, my) >= 0;
 163           if (q_is_odd)
 164             mpz_sub (r, r, my);
 165         }
 166       /* now 0 <= |r| < |my|, and if needed,
 167          q_is_odd is the least significant bit of q */
 168     }
 169
 170   if (mpz_cmp_ui (r, 0) == 0)
 171     {
 172       inex = mpfr_set_ui (rem, 0, MPFR_RNDN);
 173       /* take into account sign of x */
 174       if (signx < 0)
 175         mpfr_neg (rem, rem, MPFR_RNDN);
 176     }
 177   else
 178     {
 179       if (rnd_q == MPFR_RNDN)
 180         {
 181           /* FIXME: the comparison 2*r < my could be done more efficiently
 182              at the mpn level */
 183           mpz_mul_2exp (r, r, 1);
 184           compare = mpz_cmpabs (r, my);
 185           mpz_fdiv_q_2exp (r, r, 1);
 186           compare = ((compare > 0) ||
 187                      ((rnd_q == MPFR_RNDN) && (compare == 0) && q_is_odd));
 188           /* if compare != 0, we need to subtract my to r, and add 1 to quo */
 189           if (compare)
 190             {
 191               mpz_sub (r, r, my);
 192               if (quo && (rnd_q == MPFR_RNDN))
 193                 *quo += 1;
 194             }
 195         }
 196       /* take into account sign of x */
 197       if (signx < 0)
 198         mpz_neg (r, r);
 199       inex = mpfr_set_z_2exp (rem, r, ex > ey ? ey : ex, rnd);
 200     }
 201
 202   if (quo)
 203     *quo *= sign;
 204
 205   mpz_clear (mx);
 206   mpz_clear (my);
 207   mpz_clear (r);
 208
 209   return inex;
 210 }
 211
 212 int
 213 mpfr_remainder (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
 214 {
 215   return mpfr_rem1 (rem, (long *) 0, MPFR_RNDN, x, y, rnd);
 216 }
 217
 218 int
 219 mpfr_remquo (mpfr_ptr rem, long *quo,
 220              mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
 221 {
 222   return mpfr_rem1 (rem, quo, MPFR_RNDN, x, y, rnd);
 223 }
 224
 225 int
 226 mpfr_fmod (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
 227 {
 228   return mpfr_rem1 (rem, (long *) 0, MPFR_RNDZ, x, y, rnd);
 229 }