1/* mpfr_add1 -- internal function to perform a "real" addition
2
3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#include "mpfr-impl.h"
24
25/* compute sign(b) * (|b| + |c|), assuming b and c have same sign,
26   and are not NaN, Inf, nor zero. */
27int
28mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
29{
30  mp_limb_t *ap, *bp, *cp;
31  mpfr_prec_t aq, bq, cq, aq2;
32  mp_size_t an, bn, cn;
33  mpfr_exp_t difw, exp;
34  int sh, rb, fb, inex;
35  mpfr_uexp_t diff_exp;
36  MPFR_TMP_DECL(marker);
37
38  MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
39  MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
40
41  MPFR_TMP_MARK(marker);
42
43  aq = MPFR_PREC(a);
44  bq = MPFR_PREC(b);
45  cq = MPFR_PREC(c);
46
47  an = MPFR_PREC2LIMBS (aq); /* number of limbs of a */
48  aq2 = (mpfr_prec_t) an * GMP_NUMB_BITS;
49  sh = aq2 - aq;                  /* non-significant bits in low limb */
50
51  bn = MPFR_PREC2LIMBS (bq); /* number of limbs of b */
52  cn = MPFR_PREC2LIMBS (cq); /* number of limbs of c */
53
54  ap = MPFR_MANT(a);
55  bp = MPFR_MANT(b);
56  cp = MPFR_MANT(c);
57
58  if (MPFR_UNLIKELY(ap == bp))
59    {
60      bp = MPFR_TMP_LIMBS_ALLOC (bn);
61      MPN_COPY (bp, ap, bn);
62      if (ap == cp)
63        { cp = bp; }
64    }
65  else if (MPFR_UNLIKELY(ap == cp))
66    {
67      cp = MPFR_TMP_LIMBS_ALLOC (cn);
68      MPN_COPY(cp, ap, cn);
69    }
70
71  exp = MPFR_GET_EXP (b);
72  MPFR_SET_SAME_SIGN(a, b);
73  MPFR_UPDATE2_RND_MODE(rnd_mode, MPFR_SIGN(b));
74  /* now rnd_mode is either MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */
75  /* Note: exponents can be negative, but the unsigned subtraction is
76     a modular subtraction, so that one gets the correct result. */
77  diff_exp = (mpfr_uexp_t) exp - MPFR_GET_EXP(c);
78
79  /*
80   * 1. Compute the significant part A', the non-significant bits of A
81   * are taken into account.
82   *
83   * 2. Perform the rounding. At each iteration, we remember:
84   *     _ r = rounding bit
85   *     _ f = following bits (same value)
86   * where the result has the form: [number A]rfff...fff + a remaining
87   * value in the interval [0,2) ulp. We consider the most significant
88   * bits of the remaining value to update the result; a possible carry
89   * is immediately taken into account and A is updated accordingly. As
90   * soon as the bits f don't have the same value, A can be rounded.
91   * Variables:
92   *     _ rb = rounding bit (0 or 1).
93   *     _ fb = following bits (0 or 1), then sticky bit.
94   * If fb == 0, the only thing that can change is the sticky bit.
95   */
96
97  rb = fb = -1; /* means: not initialized */
98
99  if (MPFR_UNLIKELY (MPFR_UEXP (aq2) <= diff_exp))
100    { /* c does not overlap with a' */
101      if (MPFR_UNLIKELY(an > bn))
102        { /* a has more limbs than b */
103          /* copy b to the most significant limbs of a */
104          MPN_COPY(ap + (an - bn), bp, bn);
105          /* zero the least significant limbs of a */
106          MPN_ZERO(ap, an - bn);
107        }
108      else /* an <= bn */
109        {
110          /* copy the most significant limbs of b to a */
111          MPN_COPY(ap, bp + (bn - an), an);
112        }
113    }
114  else /* aq2 > diff_exp */
115    { /* c overlaps with a' */
116      mp_limb_t *a2p;
117      mp_limb_t cc;
118      mpfr_prec_t dif;
119      mp_size_t difn, k;
120      int shift;
121
122      /* copy c (shifted) into a */
123
124      dif = aq2 - diff_exp;
125      /* dif is the number of bits of c which overlap with a' */
126
127      difn = MPFR_PREC2LIMBS (dif);
128      /* only the highest difn limbs from c have to be considered */
129      if (MPFR_UNLIKELY(difn > cn))
130        {
131          /* c doesn't have enough limbs; take into account the virtual
132             zero limbs now by zeroing the least significant limbs of a' */
133          MPFR_ASSERTD(difn - cn <= an);
134          MPN_ZERO(ap, difn - cn);
135          difn = cn;
136        }
137      k = diff_exp / GMP_NUMB_BITS;
138
139      /* zero the most significant k limbs of a */
140      a2p = ap + (an - k);
141      MPN_ZERO(a2p, k);
142
143      shift = diff_exp % GMP_NUMB_BITS;
144
145      if (MPFR_LIKELY(shift))
146        {
147          MPFR_ASSERTD(a2p - difn >= ap);
148          cc = mpn_rshift(a2p - difn, cp + (cn - difn), difn, shift);
149          if (MPFR_UNLIKELY(a2p - difn > ap))
150            *(a2p - difn - 1) = cc;
151        }
152      else
153        MPN_COPY(a2p - difn, cp + (cn - difn), difn);
154
155      /* add b to a */
156      cc = MPFR_UNLIKELY(an > bn)
157        ? mpn_add_n(ap + (an - bn), ap + (an - bn), bp, bn)
158        : mpn_add_n(ap, ap, bp + (bn - an), an);
159
160      if (MPFR_UNLIKELY(cc)) /* carry */
161        {
162          if (MPFR_UNLIKELY(exp == __gmpfr_emax))
163            {
164              inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
165              goto end_of_add;
166            }
167          exp++;
168          rb = (ap[0] >> sh) & 1; /* LSB(a) --> rounding bit after the shift */
169          if (MPFR_LIKELY(sh))
170            {
171              mp_limb_t mask, bb;
172
173              mask = MPFR_LIMB_MASK (sh);
174              bb = ap[0] & mask;
175              ap[0] &= (~mask) << 1;
176              if (bb == 0)
177                fb = 0;
178              else if (bb == mask)
179                fb = 1;
180            }
181          mpn_rshift(ap, ap, an, 1);
182          ap[an-1] += MPFR_LIMB_HIGHBIT;
183          if (sh && fb < 0)
184            goto rounding;
185        } /* cc */
186    } /* aq2 > diff_exp */
187
188  /* non-significant bits of a */
189  if (MPFR_LIKELY(rb < 0 && sh))
190    {
191      mp_limb_t mask, bb;
192
193      mask = MPFR_LIMB_MASK (sh);
194      bb = ap[0] & mask;
195      ap[0] &= ~mask;
196      rb = bb >> (sh - 1);
197      if (MPFR_LIKELY(sh > 1))
198        {
199          mask >>= 1;
200          bb &= mask;
201          if (bb == 0)
202            fb = 0;
203          else if (bb == mask)
204            fb = 1;
205          else
206            goto rounding;
207        }
208    }
209
210  /* determine rounding and sticky bits (and possible carry) */
211
212  difw = (mpfr_exp_t) an - (mpfr_exp_t) (diff_exp / GMP_NUMB_BITS);
213  /* difw is the number of limbs from b (regarded as having an infinite
214     precision) that have already been combined with c; -n if the next
215     n limbs from b won't be combined with c. */
216
217  if (MPFR_UNLIKELY(bn > an))
218    { /* there are still limbs from b that haven't been taken into account */
219      mp_size_t bk;
220
221      if (fb == 0 && difw <= 0)
222        {
223          fb = 1; /* c hasn't been taken into account ==> sticky bit != 0 */
224          goto rounding;
225        }
226
227      bk = bn - an; /* index of lowest considered limb from b, > 0 */
228      while (difw < 0)
229        { /* ulp(next limb from b) > msb(c) */
230          mp_limb_t bb;
231
232          bb = bp[--bk];
233
234          MPFR_ASSERTD(fb != 0);
235          if (fb > 0)
236            {
237              if (bb != MP_LIMB_T_MAX)
238                {
239                  fb = 1; /* c hasn't been taken into account
240                             ==> sticky bit != 0 */
241                  goto rounding;
242                }
243            }
244          else /* fb not initialized yet */
245            {
246              if (rb < 0) /* rb not initialized yet */
247                {
248                  rb = bb >> (GMP_NUMB_BITS - 1);
249                  bb |= MPFR_LIMB_HIGHBIT;
250                }
251              fb = 1;
252              if (bb != MP_LIMB_T_MAX)
253                goto rounding;
254            }
255
256          if (bk == 0)
257            { /* b has entirely been read */
258              fb = 1; /* c hasn't been taken into account
259                         ==> sticky bit != 0 */
260              goto rounding;
261            }
262
263          difw++;
264        } /* while */
265      MPFR_ASSERTD(bk > 0 && difw >= 0);
266
267      if (difw <= cn)
268        {
269          mp_size_t ck;
270          mp_limb_t cprev;
271          int difs;
272
273          ck = cn - difw;
274          difs = diff_exp % GMP_NUMB_BITS;
275
276          if (difs == 0 && ck == 0)
277            goto c_read;
278
279          cprev = ck == cn ? 0 : cp[ck];
280
281          if (fb < 0)
282            {
283              mp_limb_t bb, cc;
284
285              if (difs)
286                {
287                  cc = cprev << (GMP_NUMB_BITS - difs);
288                  if (--ck >= 0)
289                    {
290                      cprev = cp[ck];
291                      cc += cprev >> difs;
292                    }
293                }
294              else
295                cc = cp[--ck];
296
297              bb = bp[--bk] + cc;
298
299              if (bb < cc /* carry */
300                  && (rb < 0 || (rb ^= 1) == 0)
301                  && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
302                {
303                  if (exp == __gmpfr_emax)
304                    {
305                      inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
306                      goto end_of_add;
307                    }
308                  exp++;
309                  ap[an-1] = MPFR_LIMB_HIGHBIT;
310                  rb = 0;
311                }
312
313              if (rb < 0) /* rb not initialized yet */
314                {
315                  rb = bb >> (GMP_NUMB_BITS - 1);
316                  bb <<= 1;
317                  bb |= bb >> (GMP_NUMB_BITS - 1);
318                }
319
320              fb = bb != 0;
321              if (fb && bb != MP_LIMB_T_MAX)
322                goto rounding;
323            } /* fb < 0 */
324
325          while (bk > 0)
326            {
327              mp_limb_t bb, cc;
328
329              if (difs)
330                {
331                  if (ck < 0)
332                    goto c_read;
333                  cc = cprev << (GMP_NUMB_BITS - difs);
334                  if (--ck >= 0)
335                    {
336                      cprev = cp[ck];
337                      cc += cprev >> difs;
338                    }
339                }
340              else
341                {
342                  if (ck == 0)
343                    goto c_read;
344                  cc = cp[--ck];
345                }
346
347              bb = bp[--bk] + cc;
348              if (bb < cc) /* carry */
349                {
350                  fb ^= 1;
351                  if (fb)
352                    goto rounding;
353                  rb ^= 1;
354                  if (rb == 0 && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
355                    {
356                      if (MPFR_UNLIKELY(exp == __gmpfr_emax))
357                        {
358                          inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
359                          goto end_of_add;
360                        }
361                      exp++;
362                      ap[an-1] = MPFR_LIMB_HIGHBIT;
363                    }
364                } /* bb < cc */
365
366              if (!fb && bb != 0)
367                {
368                  fb = 1;
369                  goto rounding;
370                }
371              if (fb && bb != MP_LIMB_T_MAX)
372                goto rounding;
373            } /* while */
374
375          /* b has entirely been read */
376
377          if (fb || ck < 0)
378            goto rounding;
379          if (difs && cprev << (GMP_NUMB_BITS - difs))
380            {
381              fb = 1;
382              goto rounding;
383            }
384          while (ck)
385            {
386              if (cp[--ck])
387                {
388                  fb = 1;
389                  goto rounding;
390                }
391            } /* while */
392        } /* difw <= cn */
393      else
394        { /* c has entirely been read */
395        c_read:
396          if (fb < 0) /* fb not initialized yet */
397            {
398              mp_limb_t bb;
399
400              MPFR_ASSERTD(bk > 0);
401              bb = bp[--bk];
402              if (rb < 0) /* rb not initialized yet */
403                {
404                  rb = bb >> (GMP_NUMB_BITS - 1);
405                  bb &= ~MPFR_LIMB_HIGHBIT;
406                }
407              fb = bb != 0;
408            } /* fb < 0 */
409          if (fb)
410            goto rounding;
411          while (bk)
412            {
413              if (bp[--bk])
414                {
415                  fb = 1;
416                  goto rounding;
417                }
418            } /* while */
419        } /* difw > cn */
420    } /* bn > an */
421  else if (fb != 1) /* if fb == 1, the sticky bit is 1 (no possible carry) */
422    { /* b has entirely been read */
423      if (difw > cn)
424        { /* c has entirely been read */
425          if (rb < 0)
426            rb = 0;
427          fb = 0;
428        }
429      else if (diff_exp > MPFR_UEXP (aq2))
430        { /* b is followed by at least a zero bit, then by c */
431          if (rb < 0)
432            rb = 0;
433          fb = 1;
434        }
435      else
436        {
437          mp_size_t ck;
438          int difs;
439
440          MPFR_ASSERTD(difw >= 0 && cn >= difw);
441          ck = cn - difw;
442          difs = diff_exp % GMP_NUMB_BITS;
443
444          if (difs == 0 && ck == 0)
445            { /* c has entirely been read */
446              if (rb < 0)
447                rb = 0;
448              fb = 0;
449            }
450          else
451            {
452              mp_limb_t cc;
453
454              cc = difs ? (MPFR_ASSERTD(ck < cn),
455                           cp[ck] << (GMP_NUMB_BITS - difs)) : cp[--ck];
456              if (rb < 0)
457                {
458                  rb = cc >> (GMP_NUMB_BITS - 1);
459                  cc &= ~MPFR_LIMB_HIGHBIT;
460                }
461              while (cc == 0)
462                {
463                  if (ck == 0)
464                    {
465                      fb = 0;
466                      goto rounding;
467                    }
468                  cc = cp[--ck];
469                } /* while */
470              fb = 1;
471            }
472        }
473    } /* fb != 1 */
474
475 rounding:
476  /* rnd_mode should be one of MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */
477  if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
478    {
479      if (fb == 0)
480        {
481          if (rb == 0)
482            {
483              inex = 0;
484              goto set_exponent;
485            }
486          /* round to even */
487          if (ap[0] & (MPFR_LIMB_ONE << sh))
488            goto rndn_away;
489          else
490            goto rndn_zero;
491        }
492      if (rb == 0)
493        {
494        rndn_zero:
495          inex = MPFR_IS_NEG(a) ? 1 : -1;
496          goto set_exponent;
497        }
498      else
499        {
500        rndn_away:
501          inex = MPFR_IS_POS(a) ? 1 : -1;
502          goto add_one_ulp;
503        }
504    }
505  else if (rnd_mode == MPFR_RNDZ)
506    {
507      inex = rb || fb ? (MPFR_IS_NEG(a) ? 1 : -1) : 0;
508      goto set_exponent;
509    }
510  else
511    {
512      MPFR_ASSERTN (rnd_mode == MPFR_RNDA);
513      inex = rb || fb ? (MPFR_IS_POS(a) ? 1 : -1) : 0;
514      if (inex)
515        goto add_one_ulp;
516      else
517        goto set_exponent;
518    }
519
520 add_one_ulp: /* add one unit in last place to a */
521  if (MPFR_UNLIKELY(mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh)))
522    {
523      if (MPFR_UNLIKELY(exp == __gmpfr_emax))
524        {
525          inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
526          goto end_of_add;
527        }
528      exp++;
529      ap[an-1] = MPFR_LIMB_HIGHBIT;
530    }
531
532 set_exponent:
533  MPFR_SET_EXP (a, exp);
534
535 end_of_add:
536  MPFR_TMP_FREE(marker);
537  MPFR_RET (inex);
538}
539