1/* mpfr_add1sp -- internal function to perform a "real" addition 2 All the op must have the same precision 3 4Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 5Contributed by the AriC and Caramel projects, INRIA. 6 7This file is part of the GNU MPFR Library. 8 9The GNU MPFR Library is free software; you can redistribute it and/or modify 10it under the terms of the GNU Lesser General Public License as published by 11the Free Software Foundation; either version 3 of the License, or (at your 12option) any later version. 13 14The GNU MPFR Library is distributed in the hope that it will be useful, but 15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17License for more details. 18 19You should have received a copy of the GNU Lesser General Public License 20along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2251 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 31int mpfr_add1sp2 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t); 32int 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. */ 88int 89mpfr_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