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); 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 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