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
5Copyright 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
6Contributed by the AriC and Caramel projects, INRIA.
7
8This file is part of the GNU MPFR Library.
9
10The GNU MPFR Library is free software; you can redistribute it and/or modify
11it under the terms of the GNU Lesser General Public License as published by
12the Free Software Foundation; either version 3 of the License, or (at your
13option) any later version.
14
15The GNU MPFR Library is distributed in the hope that it will be useful, but
16WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
18License for more details.
19
20You should have received a copy of the GNU Lesser General Public License
21along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
22http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2351 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
55static int
56mpfr_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
212int
213mpfr_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
218int
219mpfr_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
225int
226mpfr_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}
230