contrib/mpfr/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 Free Software Foundation, Inc.
   6 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  22 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  23 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 (GMP_RNDZ) or a remainder/remquo (GMP_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, mp_rnd_t rnd_q,
  57            mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
  58 {
  59   mp_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 == GMP_RNDN || rnd_q == GMP_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_exp (mx, x);  /* x = mx*2^ex */
  91   ey = mpfr_get_z_exp (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_div_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 == GMP_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 == GMP_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)
 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       else
 137         /* Let X = mx*2^(ex-ey) and Y = my. Then both X and Y are integers.
 138            Assume X = R mod Y, then x = X*2^ey = R*2^ey mod (Y*2^ey=y).
 139            To be able to perform the rounding, we need the least significant
 140            bit of the quotient, i.e., one more bit in the remainder,
 141            which is obtained by dividing by 2Y. */
 142         mpz_mul_2exp (my, my, 1);       /* 2Y */
 143
 144       mpz_set_ui (r, 2);
 145       mpz_powm_ui (r, r, ex - ey, my);  /* 2^(ex-ey) mod my */
 146       mpz_mul (r, r, mx);
 147       mpz_mod (r, r, my);
 148
 149       if (quo)                  /* now 0 <= r < 2^WANTED_BITS*Y */
 150         {
 151           mpz_div_2exp (my, my, WANTED_BITS);   /* back to Y */
 152           mpz_tdiv_qr (mx, r, r, my);
 153           /* oldr = mx*my + newr */
 154           *quo = mpz_get_si (mx);
 155           q_is_odd = *quo & 1;
 156         }
 157       else                      /* now 0 <= r < 2Y */
 158         {
 159           mpz_div_2exp (my, my, 1);     /* back to Y */
 160           if (rnd_q == GMP_RNDN)
 161             {
 162               /* least significant bit of q */
 163               q_is_odd = mpz_cmpabs (r, my) >= 0;
 164               if (q_is_odd)
 165                 mpz_sub (r, r, my);
 166             }
 167         }
 168       /* now 0 <= |r| < |my|, and if needed,
 169          q_is_odd is the least significant bit of q */
 170     }
 171
 172   if (mpz_cmp_ui (r, 0) == 0)
 173     inex = mpfr_set_ui (rem, 0, GMP_RNDN);
 174   else
 175     {
 176       if (rnd_q == GMP_RNDN)
 177         {
 178           /* FIXME: the comparison 2*r < my could be done more efficiently
 179              at the mpn level */
 180           mpz_mul_2exp (r, r, 1);
 181           compare = mpz_cmpabs (r, my);
 182           mpz_div_2exp (r, r, 1);
 183           compare = ((compare > 0) ||
 184                      ((rnd_q == GMP_RNDN) && (compare == 0) && q_is_odd));
 185           /* if compare != 0, we need to subtract my to r, and add 1 to quo */
 186           if (compare)
 187             {
 188               mpz_sub (r, r, my);
 189               if (quo && (rnd_q == GMP_RNDN))
 190                 *quo += 1;
 191             }
 192         }
 193       inex = mpfr_set_z (rem, r, rnd);
 194       /* if ex > ey, rem should be multiplied by 2^ey, else by 2^ex */
 195       MPFR_EXP (rem) += (ex > ey) ? ey : ex;
 196     }
 197
 198   if (quo)
 199     *quo *= sign;
 200
 201   /* take into account sign of x */
 202   if (signx < 0)
 203     {
 204       mpfr_neg (rem, rem, GMP_RNDN);
 205       inex = -inex;
 206     }
 207
 208   mpz_clear (mx);
 209   mpz_clear (my);
 210   mpz_clear (r);
 211
 212   return inex;
 213 }
 214
 215 int
 216 mpfr_remainder (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
 217 {
 218   return mpfr_rem1 (rem, (long *) 0, GMP_RNDN, x, y, rnd);
 219 }
 220
 221 int
 222 mpfr_remquo (mpfr_ptr rem, long *quo,
 223              mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
 224 {
 225   return mpfr_rem1 (rem, quo, GMP_RNDN, x, y, rnd);
 226 }
 227
 228 int
 229 mpfr_fmod (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
 230 {
 231   return mpfr_rem1 (rem, (long *) 0, GMP_RNDZ, x, y, rnd);
 232 }