1 /* mpfr_add1sp -- internal function to perform a "real" addition
2 All the op must have the same precision
3
4 Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramel projects, INRIA.
6
7 This file is part of the GNU MPFR Library.
8
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23
24 #define MPFR_NEED_LONGLONG_H
25 #include "mpfr-impl.h"
26
27 /* Check if we have to check the result of mpfr_add1sp with mpfr_add1 */
28 #ifdef WANT_ASSERT
29 # if WANT_ASSERT >= 2
30
31 int mpfr_add1sp2 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
mpfr_add1sp(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)32 int mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
33 {
34 mpfr_t tmpa, tmpb, tmpc;
35 int inexb, inexc, inexact, inexact2;
36
37 mpfr_init2 (tmpa, MPFR_PREC (a));
38 mpfr_init2 (tmpb, MPFR_PREC (b));
39 mpfr_init2 (tmpc, MPFR_PREC (c));
40
41 inexb = mpfr_set (tmpb, b, MPFR_RNDN);
42 MPFR_ASSERTN (inexb == 0);
43
44 inexc = mpfr_set (tmpc, c, MPFR_RNDN);
45 MPFR_ASSERTN (inexc == 0);
46
47 inexact2 = mpfr_add1 (tmpa, tmpb, tmpc, rnd_mode);
48 inexact = mpfr_add1sp2 (a, b, c, rnd_mode);
49
50 if (mpfr_cmp (tmpa, a) || inexact != inexact2)
51 {
52 fprintf (stderr, "add1 & add1sp return different values for %s\n"
53 "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
54 mpfr_print_rnd_mode (rnd_mode),
55 (unsigned long) MPFR_PREC (a),
56 (unsigned long) MPFR_PREC (b),
57 (unsigned long) MPFR_PREC (c));
58 mpfr_fprint_binary (stderr, tmpb);
59 fprintf (stderr, "\nC = ");
60 mpfr_fprint_binary (stderr, tmpc);
61 fprintf (stderr, "\n\nadd1 : ");
62 mpfr_fprint_binary (stderr, tmpa);
63 fprintf (stderr, "\nadd1sp: ");
64 mpfr_fprint_binary (stderr, a);
65 fprintf (stderr, "\nInexact sp = %d | Inexact = %d\n",
66 inexact, inexact2);
67 MPFR_ASSERTN (0);
68 }
69 mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
70 return inexact;
71 }
72 # define mpfr_add1sp mpfr_add1sp2
73 # endif
74 #endif
75
76 /* Debugging support */
77 #ifdef DEBUG
78 # undef DEBUG
79 # define DEBUG(x) (x)
80 #else
81 # define DEBUG(x) /**/
82 #endif
83
84 /* compute sign(b) * (|b| + |c|)
85 Returns 0 iff result is exact,
86 a negative value when the result is less than the exact value,
87 a positive value otherwise. */
88 int
mpfr_add1sp(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)89 mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
90 {
91 mpfr_uexp_t d;
92 mpfr_prec_t p;
93 unsigned int sh;
94 mp_size_t n;
95 mp_limb_t *ap, *cp;
96 mpfr_exp_t bx;
97 mp_limb_t limb;
98 int inexact;
99 MPFR_TMP_DECL(marker);
100
101 MPFR_TMP_MARK(marker);
102
103 MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
104 MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
105 MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
106 MPFR_ASSERTD(MPFR_GET_EXP(b) >= MPFR_GET_EXP(c));
107
108 /* Read prec and num of limbs */
109 p = MPFR_PREC(b);
110 n = MPFR_PREC2LIMBS (p);
111 MPFR_UNSIGNED_MINUS_MODULO(sh, p);
112 bx = MPFR_GET_EXP(b);
113 d = (mpfr_uexp_t) (bx - MPFR_GET_EXP(c));
114
115 DEBUG (printf ("New add1sp with diff=%lu\n", (unsigned long) d));
116
117 if (MPFR_UNLIKELY(d == 0))
118 {
119 /* d==0 */
120 DEBUG( mpfr_print_mant_binary("C= ", MPFR_MANT(c), p) );
121 DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
122 bx++; /* exp + 1 */
123 ap = MPFR_MANT(a);
124 limb = mpn_add_n(ap, MPFR_MANT(b), MPFR_MANT(c), n);
125 DEBUG( mpfr_print_mant_binary("A= ", ap, p) );
126 MPFR_ASSERTD(limb != 0); /* There must be a carry */
127 limb = ap[0]; /* Get LSB (In fact, LSW) */
128 mpn_rshift(ap, ap, n, 1); /* Shift mantissa A */
129 ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */
130 ap[0] &= ~MPFR_LIMB_MASK(sh); /* Clear LSB bit */
131 if (MPFR_LIKELY((limb&(MPFR_LIMB_ONE<<sh)) == 0)) /* Check exact case */
132 { inexact = 0; goto set_exponent; }
133 /* Zero: Truncate
134 Nearest: Even Rule => truncate or add 1
135 Away: Add 1 */
136 if (MPFR_LIKELY(rnd_mode==MPFR_RNDN))
137 {
138 if (MPFR_LIKELY((ap[0]&(MPFR_LIMB_ONE<<sh))==0))
139 { inexact = -1; goto set_exponent; }
140 else
141 goto add_one_ulp;
142 }
143 MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
144 if (rnd_mode==MPFR_RNDZ)
145 { inexact = -1; goto set_exponent; }
146 else
147 goto add_one_ulp;
148 }
149 else if (MPFR_UNLIKELY (d >= p))
150 {
151 if (MPFR_LIKELY (d > p))
152 {
153 /* d > p : Copy B in A */
154 /* Away: Add 1
155 Nearest: Trunc
156 Zero: Trunc */
157 if (MPFR_LIKELY (rnd_mode==MPFR_RNDN
158 || MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b))))
159 {
160 copy_set_exponent:
161 ap = MPFR_MANT (a);
162 MPN_COPY (ap, MPFR_MANT(b), n);
163 inexact = -1;
164 goto set_exponent;
165 }
166 else
167 {
168 copy_add_one_ulp:
169 ap = MPFR_MANT(a);
170 MPN_COPY (ap, MPFR_MANT(b), n);
171 goto add_one_ulp;
172 }
173 }
174 else
175 {
176 /* d==p : Copy B in A */
177 /* Away: Add 1
178 Nearest: Even Rule if C is a power of 2, else Add 1
179 Zero: Trunc */
180 if (MPFR_LIKELY(rnd_mode==MPFR_RNDN))
181 {
182 /* Check if C was a power of 2 */
183 cp = MPFR_MANT(c);
184 if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
185 {
186 mp_size_t k = n-1;
187 do {
188 k--;
189 } while (k>=0 && cp[k]==0);
190 if (MPFR_UNLIKELY(k<0))
191 /* Power of 2: Even rule */
192 if ((MPFR_MANT (b)[0]&(MPFR_LIMB_ONE<<sh))==0)
193 goto copy_set_exponent;
194 }
195 /* Not a Power of 2 */
196 goto copy_add_one_ulp;
197 }
198 else if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b)))
199 goto copy_set_exponent;
200 else
201 goto copy_add_one_ulp;
202 }
203 }
204 else
205 {
206 mp_limb_t mask;
207 mp_limb_t bcp, bcp1; /* Cp and C'p+1 */
208
209 /* General case: 1 <= d < p */
210 cp = MPFR_TMP_LIMBS_ALLOC (n);
211
212 /* Shift c in temporary allocated place */
213 {
214 mpfr_uexp_t dm;
215 mp_size_t m;
216
217 dm = d % GMP_NUMB_BITS;
218 m = d / GMP_NUMB_BITS;
219 if (MPFR_UNLIKELY(dm == 0))
220 {
221 /* dm = 0 and m > 0: Just copy */
222 MPFR_ASSERTD(m!=0);
223 MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
224 MPN_ZERO(cp+n-m, m);
225 }
226 else if (MPFR_LIKELY(m == 0))
227 {
228 /* dm >=1 and m == 0: just shift */
229 MPFR_ASSERTD(dm >= 1);
230 mpn_rshift(cp, MPFR_MANT(c), n, dm);
231 }
232 else
233 {
234 /* dm > 0 and m > 0: shift and zero */
235 mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
236 MPN_ZERO(cp+n-m, m);
237 }
238 }
239
240 DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
241 DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
242 DEBUG( mpfr_print_mant_binary("After ", cp, p) );
243
244 /* Compute bcp=Cp and bcp1=C'p+1 */
245 if (MPFR_LIKELY (sh > 0))
246 {
247 /* Try to compute them from C' rather than C */
248 bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
249 if (MPFR_LIKELY(cp[0]&MPFR_LIMB_MASK(sh-1)))
250 bcp1 = 1;
251 else
252 {
253 /* We can't compute C'p+1 from C'. Compute it from C */
254 /* Start from bit x=p-d+sh in mantissa C
255 (+sh since we have already looked sh bits in C'!) */
256 mpfr_prec_t x = p-d+sh-1;
257 if (MPFR_LIKELY(x>p))
258 /* We are already looked at all the bits of c, so C'p+1 = 0*/
259 bcp1 = 0;
260 else
261 {
262 mp_limb_t *tp = MPFR_MANT(c);
263 mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
264 mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
265 DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
266 (unsigned long) x, (long) kx,
267 (unsigned long) sx));
268 /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
269 if (tp[kx] & MPFR_LIMB_MASK(sx))
270 bcp1 = 1;
271 else
272 {
273 /*kx += (sx==0);*/
274 /*If sx==0, tp[kx] hasn't been checked*/
275 do {
276 kx--;
277 } while (kx>=0 && tp[kx]==0);
278 bcp1 = (kx >= 0);
279 }
280 }
281 }
282 }
283 else /* sh == 0 */
284 {
285 /* Compute Cp and C'p+1 from C with sh=0 */
286 mp_limb_t *tp = MPFR_MANT(c);
287 /* Start from bit x=p-d in mantissa C */
288 mpfr_prec_t x = p-d;
289 mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
290 mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
291 MPFR_ASSERTD(p >= d);
292 bcp = tp[kx] & (MPFR_LIMB_ONE<<sx);
293 /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
294 if (tp[kx]&MPFR_LIMB_MASK(sx))
295 bcp1 = 1;
296 else
297 {
298 do {
299 kx--;
300 } while (kx>=0 && tp[kx]==0);
301 bcp1 = (kx>=0);
302 }
303 }
304 DEBUG (printf("sh=%u Cp=%lu C'p+1=%lu\n", sh,
305 (unsigned long) bcp, (unsigned long) bcp1));
306
307 /* Clean shifted C' */
308 mask = ~MPFR_LIMB_MASK(sh);
309 cp[0] &= mask;
310
311 /* Add the mantissa c from b in a */
312 ap = MPFR_MANT(a);
313 limb = mpn_add_n (ap, MPFR_MANT(b), cp, n);
314 DEBUG( mpfr_print_mant_binary("Add= ", ap, p) );
315
316 /* Check for overflow */
317 if (MPFR_UNLIKELY (limb))
318 {
319 limb = ap[0] & (MPFR_LIMB_ONE<<sh); /* Get LSB */
320 mpn_rshift (ap, ap, n, 1); /* Shift mantissa*/
321 bx++; /* Fix exponent */
322 ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */
323 ap[0] &= mask; /* Clear LSB bit */
324 bcp1 |= bcp; /* Recompute C'p+1 */
325 bcp = limb; /* Recompute Cp */
326 DEBUG (printf ("(Overflow) Cp=%lu C'p+1=%lu\n",
327 (unsigned long) bcp, (unsigned long) bcp1));
328 DEBUG (mpfr_print_mant_binary ("Add= ", ap, p));
329 }
330
331 /* Round:
332 Zero: Truncate but could be exact.
333 Away: Add 1 if Cp or C'p+1 !=0
334 Nearest: Truncate but could be exact if Cp==0
335 Add 1 if C'p+1 !=0,
336 Even rule else */
337 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
338 {
339 if (MPFR_LIKELY(bcp == 0))
340 { inexact = MPFR_LIKELY(bcp1) ? -1 : 0; goto set_exponent; }
341 else if (MPFR_UNLIKELY(bcp1==0) && (ap[0]&(MPFR_LIMB_ONE<<sh))==0)
342 { inexact = -1; goto set_exponent; }
343 else
344 goto add_one_ulp;
345 }
346 MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
347 if (rnd_mode == MPFR_RNDZ)
348 {
349 inexact = MPFR_LIKELY(bcp || bcp1) ? -1 : 0;
350 goto set_exponent;
351 }
352 else
353 {
354 if (MPFR_UNLIKELY(bcp==0 && bcp1==0))
355 { inexact = 0; goto set_exponent; }
356 else
357 goto add_one_ulp;
358 }
359 }
360 MPFR_ASSERTN(0);
361
362 add_one_ulp:
363 /* add one unit in last place to a */
364 DEBUG( printf("AddOneUlp\n") );
365 if (MPFR_UNLIKELY( mpn_add_1(ap, ap, n, MPFR_LIMB_ONE<<sh) ))
366 {
367 /* Case 100000x0 = 0x1111x1 + 1*/
368 DEBUG( printf("Pow of 2\n") );
369 bx++;
370 ap[n-1] = MPFR_LIMB_HIGHBIT;
371 }
372 inexact = 1;
373
374 set_exponent:
375 if (MPFR_UNLIKELY(bx > __gmpfr_emax)) /* Check for overflow */
376 {
377 DEBUG( printf("Overflow\n") );
378 MPFR_TMP_FREE(marker);
379 MPFR_SET_SAME_SIGN(a,b);
380 return mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
381 }
382 MPFR_SET_EXP (a, bx);
383 MPFR_SET_SAME_SIGN(a,b);
384
385 MPFR_TMP_FREE(marker);
386 MPFR_RET (inexact * MPFR_INT_SIGN (a));
387 }
388