1 /* mpfr_sub1sp -- internal function to perform a "real" substraction
2 All the op must have the same precision
3
4 Copyright 2003, 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_sub1sp with mpfr_sub1 */
28 #ifdef WANT_ASSERT
29 # if WANT_ASSERT >= 2
30
31 int mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode);
mpfr_sub1sp(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)32 int mpfr_sub1sp (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_sub1 (tmpa, tmpb, tmpc, rnd_mode);
48 inexact = mpfr_sub1sp2(a, b, c, rnd_mode);
49
50 if (mpfr_cmp (tmpa, a) || inexact != inexact2)
51 {
52 fprintf (stderr, "sub1 & sub1sp return different values for %s\n"
53 "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
54 mpfr_print_rnd_mode (rnd_mode), (unsigned long) MPFR_PREC (a),
55 (unsigned long) MPFR_PREC (b), (unsigned long) MPFR_PREC (c));
56 mpfr_fprint_binary (stderr, tmpb);
57 fprintf (stderr, "\nC = ");
58 mpfr_fprint_binary (stderr, tmpc);
59 fprintf (stderr, "\nSub1 : ");
60 mpfr_fprint_binary (stderr, tmpa);
61 fprintf (stderr, "\nSub1sp: ");
62 mpfr_fprint_binary (stderr, a);
63 fprintf (stderr, "\nInexact sp = %d | Inexact = %d\n",
64 inexact, inexact2);
65 MPFR_ASSERTN (0);
66 }
67 mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
68 return inexact;
69 }
70 # define mpfr_sub1sp mpfr_sub1sp2
71 # endif
72 #endif
73
74 /* Debugging support */
75 #ifdef DEBUG
76 # undef DEBUG
77 # define DEBUG(x) (x)
78 #else
79 # define DEBUG(x) /**/
80 #endif
81
82 /* Rounding Sub */
83
84 /*
85 compute sgn(b)*(|b| - |c|) if |b|>|c| else -sgn(b)*(|c| -|b|)
86 Returns 0 iff result is exact,
87 a negative value when the result is less than the exact value,
88 a positive value otherwise.
89 */
90
91 /* A0...Ap-1
92 * Cp Cp+1 ....
93 * <- C'p+1 ->
94 * Cp = -1 if calculated from c mantissa
95 * Cp = 0 if 0 from a or c
96 * Cp = 1 if calculated from a.
97 * C'p+1 = First bit not null or 0 if there isn't one
98 *
99 * Can't have Cp=-1 and C'p+1=1*/
100
101 /* RND = MPFR_RNDZ:
102 * + if Cp=0 and C'p+1=0,1, Truncate.
103 * + if Cp=0 and C'p+1=-1, SubOneUlp
104 * + if Cp=-1, SubOneUlp
105 * + if Cp=1, AddOneUlp
106 * RND = MPFR_RNDA (Away)
107 * + if Cp=0 and C'p+1=0,-1, Truncate
108 * + if Cp=0 and C'p+1=1, AddOneUlp
109 * + if Cp=1, AddOneUlp
110 * + if Cp=-1, Truncate
111 * RND = MPFR_RNDN
112 * + if Cp=0, Truncate
113 * + if Cp=1 and C'p+1=1, AddOneUlp
114 * + if Cp=1 and C'p+1=-1, Truncate
115 * + if Cp=1 and C'p+1=0, Truncate if Ap-1=0, AddOneUlp else
116 * + if Cp=-1 and C'p+1=-1, SubOneUlp
117 * + if Cp=-1 and C'p+1=0, Truncate if Ap-1=0, SubOneUlp else
118 *
119 * If AddOneUlp:
120 * If carry, then it is 11111111111 + 1 = 10000000000000
121 * ap[n-1]=MPFR_HIGHT_BIT
122 * If SubOneUlp:
123 * If we lose one bit, it is 1000000000 - 1 = 0111111111111
124 * Then shift, and put as last bit x which is calculated
125 * according Cp, Cp-1 and rnd_mode.
126 * If Truncate,
127 * If it is a power of 2,
128 * we may have to suboneulp in some special cases.
129 *
130 * To simplify, we don't use Cp = 1.
131 *
132 */
133
134 int
mpfr_sub1sp(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)135 mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
136 {
137 mpfr_exp_t bx,cx;
138 mpfr_uexp_t d;
139 mpfr_prec_t p, sh, cnt;
140 mp_size_t n;
141 mp_limb_t *ap, *bp, *cp;
142 mp_limb_t limb;
143 int inexact;
144 mp_limb_t bcp,bcp1; /* Cp and C'p+1 */
145 mp_limb_t bbcp = (mp_limb_t) -1, bbcp1 = (mp_limb_t) -1; /* Cp+1 and C'p+2,
146 gcc claims that they might be used uninitialized. We fill them with invalid
147 values, which should produce a failure if so. See README.dev file. */
148
149 MPFR_TMP_DECL(marker);
150
151 MPFR_TMP_MARK(marker);
152
153 MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
154 MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
155 MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
156
157 /* Read prec and num of limbs */
158 p = MPFR_PREC (b);
159 n = MPFR_PREC2LIMBS (p);
160
161 /* Fast cmp of |b| and |c|*/
162 bx = MPFR_GET_EXP (b);
163 cx = MPFR_GET_EXP (c);
164 if (MPFR_UNLIKELY(bx == cx))
165 {
166 mp_size_t k = n - 1;
167 /* Check mantissa since exponent are equals */
168 bp = MPFR_MANT(b);
169 cp = MPFR_MANT(c);
170 while (k>=0 && MPFR_UNLIKELY(bp[k] == cp[k]))
171 k--;
172 if (MPFR_UNLIKELY(k < 0))
173 /* b == c ! */
174 {
175 /* Return exact number 0 */
176 if (rnd_mode == MPFR_RNDD)
177 MPFR_SET_NEG(a);
178 else
179 MPFR_SET_POS(a);
180 MPFR_SET_ZERO(a);
181 MPFR_RET(0);
182 }
183 else if (bp[k] > cp[k])
184 goto BGreater;
185 else
186 {
187 MPFR_ASSERTD(bp[k]<cp[k]);
188 goto CGreater;
189 }
190 }
191 else if (MPFR_UNLIKELY(bx < cx))
192 {
193 /* Swap b and c and set sign */
194 mpfr_srcptr t;
195 mpfr_exp_t tx;
196 CGreater:
197 MPFR_SET_OPPOSITE_SIGN(a,b);
198 t = b; b = c; c = t;
199 tx = bx; bx = cx; cx = tx;
200 }
201 else
202 {
203 /* b > c */
204 BGreater:
205 MPFR_SET_SAME_SIGN(a,b);
206 }
207
208 /* Now b > c */
209 MPFR_ASSERTD(bx >= cx);
210 d = (mpfr_uexp_t) bx - cx;
211 DEBUG (printf ("New with diff=%lu\n", (unsigned long) d));
212
213 if (MPFR_UNLIKELY(d <= 1))
214 {
215 if (MPFR_LIKELY(d < 1))
216 {
217 /* <-- b -->
218 <-- c --> : exact sub */
219 ap = MPFR_MANT(a);
220 mpn_sub_n (ap, MPFR_MANT(b), MPFR_MANT(c), n);
221 /* Normalize */
222 ExactNormalize:
223 limb = ap[n-1];
224 if (MPFR_LIKELY(limb))
225 {
226 /* First limb is not zero. */
227 count_leading_zeros(cnt, limb);
228 /* cnt could be == 0 <= SubD1Lose */
229 if (MPFR_LIKELY(cnt))
230 {
231 mpn_lshift(ap, ap, n, cnt); /* Normalize number */
232 bx -= cnt; /* Update final expo */
233 }
234 /* Last limb should be ok */
235 MPFR_ASSERTD(!(ap[0] & MPFR_LIMB_MASK((unsigned int) (-p)
236 % GMP_NUMB_BITS)));
237 }
238 else
239 {
240 /* First limb is zero */
241 mp_size_t k = n-1, len;
242 /* Find the first limb not equal to zero.
243 FIXME:It is assume it exists (since |b| > |c| and same prec)*/
244 do
245 {
246 MPFR_ASSERTD( k > 0 );
247 limb = ap[--k];
248 }
249 while (limb == 0);
250 MPFR_ASSERTD(limb != 0);
251 count_leading_zeros(cnt, limb);
252 k++;
253 len = n - k; /* Number of last limb */
254 MPFR_ASSERTD(k >= 0);
255 if (MPFR_LIKELY(cnt))
256 mpn_lshift(ap+len, ap, k, cnt); /* Normalize the High Limb*/
257 else
258 {
259 /* Must use DECR since src and dest may overlap & dest>=src*/
260 MPN_COPY_DECR(ap+len, ap, k);
261 }
262 MPN_ZERO(ap, len); /* Zeroing the last limbs */
263 bx -= cnt + len*GMP_NUMB_BITS; /* Update Expo */
264 /* Last limb should be ok */
265 MPFR_ASSERTD(!(ap[len]&MPFR_LIMB_MASK((unsigned int) (-p)
266 % GMP_NUMB_BITS)));
267 }
268 /* Check expo underflow */
269 if (MPFR_UNLIKELY(bx < __gmpfr_emin))
270 {
271 MPFR_TMP_FREE(marker);
272 /* inexact=0 */
273 DEBUG( printf("(D==0 Underflow)\n") );
274 if (rnd_mode == MPFR_RNDN &&
275 (bx < __gmpfr_emin - 1 ||
276 (/*inexact >= 0 &&*/ mpfr_powerof2_raw (a))))
277 rnd_mode = MPFR_RNDZ;
278 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
279 }
280 MPFR_SET_EXP (a, bx);
281 /* No rounding is necessary since the result is exact */
282 MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
283 MPFR_TMP_FREE(marker);
284 return 0;
285 }
286 else /* if (d == 1) */
287 {
288 /* | <-- b -->
289 | <-- c --> */
290 mp_limb_t c0, mask;
291 mp_size_t k;
292 MPFR_UNSIGNED_MINUS_MODULO(sh, p);
293 /* If we lose at least one bit, compute 2*b-c (Exact)
294 * else compute b-c/2 */
295 bp = MPFR_MANT(b);
296 cp = MPFR_MANT(c);
297 k = n-1;
298 limb = bp[k] - cp[k]/2;
299 if (limb > MPFR_LIMB_HIGHBIT)
300 {
301 /* We can't lose precision: compute b-c/2 */
302 /* Shift c in the allocated temporary block */
303 SubD1NoLose:
304 c0 = cp[0] & (MPFR_LIMB_ONE<<sh);
305 cp = MPFR_TMP_LIMBS_ALLOC (n);
306 mpn_rshift(cp, MPFR_MANT(c), n, 1);
307 if (MPFR_LIKELY(c0 == 0))
308 {
309 /* Result is exact: no need of rounding! */
310 ap = MPFR_MANT(a);
311 mpn_sub_n (ap, bp, cp, n);
312 MPFR_SET_EXP(a, bx); /* No expo overflow! */
313 /* No truncate or normalize is needed */
314 MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
315 /* No rounding is necessary since the result is exact */
316 MPFR_TMP_FREE(marker);
317 return 0;
318 }
319 ap = MPFR_MANT(a);
320 mask = ~MPFR_LIMB_MASK(sh);
321 cp[0] &= mask; /* Delete last bit of c */
322 mpn_sub_n (ap, bp, cp, n);
323 MPFR_SET_EXP(a, bx); /* No expo overflow! */
324 MPFR_ASSERTD( !(ap[0] & ~mask) ); /* Check last bits */
325 /* No normalize is needed */
326 MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
327 /* Rounding is necessary since c0 = 1*/
328 /* Cp =-1 and C'p+1=0 */
329 bcp = 1; bcp1 = 0;
330 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
331 {
332 /* Even Rule apply: Check Ap-1 */
333 if (MPFR_LIKELY( (ap[0] & (MPFR_LIMB_ONE<<sh)) == 0) )
334 goto truncate;
335 else
336 goto sub_one_ulp;
337 }
338 MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
339 if (rnd_mode == MPFR_RNDZ)
340 goto sub_one_ulp;
341 else
342 goto truncate;
343 }
344 else if (MPFR_LIKELY(limb < MPFR_LIMB_HIGHBIT))
345 {
346 /* We lose at least one bit of prec */
347 /* Calcul of 2*b-c (Exact) */
348 /* Shift b in the allocated temporary block */
349 SubD1Lose:
350 bp = MPFR_TMP_LIMBS_ALLOC (n);
351 mpn_lshift (bp, MPFR_MANT(b), n, 1);
352 ap = MPFR_MANT(a);
353 mpn_sub_n (ap, bp, cp, n);
354 bx--;
355 goto ExactNormalize;
356 }
357 else
358 {
359 /* Case: limb = 100000000000 */
360 /* Check while b[k] == c'[k] (C' is C shifted by 1) */
361 /* If b[k]<c'[k] => We lose at least one bit*/
362 /* If b[k]>c'[k] => We don't lose any bit */
363 /* If k==-1 => We don't lose any bit
364 AND the result is 100000000000 0000000000 00000000000 */
365 mp_limb_t carry;
366 do {
367 carry = cp[k]&MPFR_LIMB_ONE;
368 k--;
369 } while (k>=0 &&
370 bp[k]==(carry=cp[k]/2+(carry<<(GMP_NUMB_BITS-1))));
371 if (MPFR_UNLIKELY(k<0))
372 {
373 /*If carry then (sh==0 and Virtual c'[-1] > Virtual b[-1]) */
374 if (MPFR_UNLIKELY(carry)) /* carry = cp[0]&MPFR_LIMB_ONE */
375 {
376 /* FIXME: Can be faster? */
377 MPFR_ASSERTD(sh == 0);
378 goto SubD1Lose;
379 }
380 /* Result is a power of 2 */
381 ap = MPFR_MANT (a);
382 MPN_ZERO (ap, n);
383 ap[n-1] = MPFR_LIMB_HIGHBIT;
384 MPFR_SET_EXP (a, bx); /* No expo overflow! */
385 /* No Normalize is needed*/
386 /* No Rounding is needed */
387 MPFR_TMP_FREE (marker);
388 return 0;
389 }
390 /* carry = cp[k]/2+(cp[k-1]&1)<<(GMP_NUMB_BITS-1) = c'[k]*/
391 else if (bp[k] > carry)
392 goto SubD1NoLose;
393 else
394 {
395 MPFR_ASSERTD(bp[k]<carry);
396 goto SubD1Lose;
397 }
398 }
399 }
400 }
401 else if (MPFR_UNLIKELY(d >= p))
402 {
403 ap = MPFR_MANT(a);
404 MPFR_UNSIGNED_MINUS_MODULO(sh, p);
405 /* We can't set A before since we use cp for rounding... */
406 /* Perform rounding: check if a=b or a=b-ulp(b) */
407 if (MPFR_UNLIKELY(d == p))
408 {
409 /* cp == -1 and c'p+1 = ? */
410 bcp = 1;
411 /* We need Cp+1 later for a very improbable case. */
412 bbcp = (MPFR_MANT(c)[n-1] & (MPFR_LIMB_ONE<<(GMP_NUMB_BITS-2)));
413 /* We need also C'p+1 for an even more unprobable case... */
414 if (MPFR_LIKELY( bbcp ))
415 bcp1 = 1;
416 else
417 {
418 cp = MPFR_MANT(c);
419 if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
420 {
421 mp_size_t k = n-1;
422 do {
423 k--;
424 } while (k>=0 && cp[k]==0);
425 bcp1 = (k>=0);
426 }
427 else
428 bcp1 = 1;
429 }
430 DEBUG( printf("(D=P) Cp=-1 Cp+1=%d C'p+1=%d \n", bbcp!=0, bcp1!=0) );
431 bp = MPFR_MANT (b);
432
433 /* Even if src and dest overlap, it is ok using MPN_COPY */
434 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
435 {
436 if (MPFR_UNLIKELY( bcp && bcp1==0 ))
437 /* Cp=-1 and C'p+1=0: Even rule Apply! */
438 /* Check Ap-1 = Bp-1 */
439 if ((bp[0] & (MPFR_LIMB_ONE<<sh)) == 0)
440 {
441 MPN_COPY(ap, bp, n);
442 goto truncate;
443 }
444 MPN_COPY(ap, bp, n);
445 goto sub_one_ulp;
446 }
447 MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
448 if (rnd_mode == MPFR_RNDZ)
449 {
450 MPN_COPY(ap, bp, n);
451 goto sub_one_ulp;
452 }
453 else
454 {
455 MPN_COPY(ap, bp, n);
456 goto truncate;
457 }
458 }
459 else
460 {
461 /* Cp=0, Cp+1=-1 if d==p+1, C'p+1=-1 */
462 bcp = 0; bbcp = (d==p+1); bcp1 = 1;
463 DEBUG( printf("(D>P) Cp=%d Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0) );
464 /* Need to compute C'p+2 if d==p+1 and if rnd_mode=NEAREST
465 (Because of a very improbable case) */
466 if (MPFR_UNLIKELY(d==p+1 && rnd_mode==MPFR_RNDN))
467 {
468 cp = MPFR_MANT(c);
469 if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
470 {
471 mp_size_t k = n-1;
472 do {
473 k--;
474 } while (k>=0 && cp[k]==0);
475 bbcp1 = (k>=0);
476 }
477 else
478 bbcp1 = 1;
479 DEBUG( printf("(D>P) C'p+2=%d\n", bbcp1!=0) );
480 }
481 /* Copy mantissa B in A */
482 MPN_COPY(ap, MPFR_MANT(b), n);
483 /* Round */
484 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
485 goto truncate;
486 MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
487 if (rnd_mode == MPFR_RNDZ)
488 goto sub_one_ulp;
489 else /* rnd_mode = AWAY */
490 goto truncate;
491 }
492 }
493 else
494 {
495 mpfr_uexp_t dm;
496 mp_size_t m;
497 mp_limb_t mask;
498
499 /* General case: 2 <= d < p */
500 MPFR_UNSIGNED_MINUS_MODULO(sh, p);
501 cp = MPFR_TMP_LIMBS_ALLOC (n);
502
503 /* Shift c in temporary allocated place */
504 dm = d % GMP_NUMB_BITS;
505 m = d / GMP_NUMB_BITS;
506 if (MPFR_UNLIKELY(dm == 0))
507 {
508 /* dm = 0 and m > 0: Just copy */
509 MPFR_ASSERTD(m!=0);
510 MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
511 MPN_ZERO(cp+n-m, m);
512 }
513 else if (MPFR_LIKELY(m == 0))
514 {
515 /* dm >=2 and m == 0: just shift */
516 MPFR_ASSERTD(dm >= 2);
517 mpn_rshift(cp, MPFR_MANT(c), n, dm);
518 }
519 else
520 {
521 /* dm > 0 and m > 0: shift and zero */
522 mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
523 MPN_ZERO(cp+n-m, m);
524 }
525
526 DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
527 DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
528 DEBUG( mpfr_print_mant_binary("After ", cp, p) );
529
530 /* Compute bcp=Cp and bcp1=C'p+1 */
531 if (MPFR_LIKELY(sh))
532 {
533 /* Try to compute them from C' rather than C (FIXME: Faster?) */
534 bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
535 if (MPFR_LIKELY( cp[0] & MPFR_LIMB_MASK(sh-1) ))
536 bcp1 = 1;
537 else
538 {
539 /* We can't compute C'p+1 from C'. Compute it from C */
540 /* Start from bit x=p-d+sh in mantissa C
541 (+sh since we have already looked sh bits in C'!) */
542 mpfr_prec_t x = p-d+sh-1;
543 if (MPFR_LIKELY(x>p))
544 /* We are already looked at all the bits of c, so C'p+1 = 0*/
545 bcp1 = 0;
546 else
547 {
548 mp_limb_t *tp = MPFR_MANT(c);
549 mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
550 mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
551 DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
552 (unsigned long) x, (long) kx,
553 (unsigned long) sx));
554 /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
555 if (tp[kx] & MPFR_LIMB_MASK(sx))
556 bcp1 = 1;
557 else
558 {
559 /*kx += (sx==0);*/
560 /*If sx==0, tp[kx] hasn't been checked*/
561 do {
562 kx--;
563 } while (kx>=0 && tp[kx]==0);
564 bcp1 = (kx >= 0);
565 }
566 }
567 }
568 }
569 else
570 {
571 /* Compute Cp and C'p+1 from C with sh=0 */
572 mp_limb_t *tp = MPFR_MANT(c);
573 /* Start from bit x=p-d in mantissa C */
574 mpfr_prec_t x = p-d;
575 mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
576 mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
577 MPFR_ASSERTD(p >= d);
578 bcp = (tp[kx] & (MPFR_LIMB_ONE<<sx));
579 /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
580 if (tp[kx] & MPFR_LIMB_MASK(sx))
581 bcp1 = 1;
582 else
583 {
584 /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
585 do {
586 kx--;
587 } while (kx>=0 && tp[kx]==0);
588 bcp1 = (kx>=0);
589 }
590 }
591 DEBUG( printf("sh=%lu Cp=%d C'p+1=%d\n", sh, bcp!=0, bcp1!=0) );
592
593 /* Check if we can lose a bit, and if so compute Cp+1 and C'p+2 */
594 bp = MPFR_MANT(b);
595 if (MPFR_UNLIKELY((bp[n-1]-cp[n-1]) <= MPFR_LIMB_HIGHBIT))
596 {
597 /* We can lose a bit so we precompute Cp+1 and C'p+2 */
598 /* Test for trivial case: since C'p+1=0, Cp+1=0 and C'p+2 =0 */
599 if (MPFR_LIKELY(bcp1 == 0))
600 {
601 bbcp = 0;
602 bbcp1 = 0;
603 }
604 else /* bcp1 != 0 */
605 {
606 /* We can lose a bit:
607 compute Cp+1 and C'p+2 from mantissa C */
608 mp_limb_t *tp = MPFR_MANT(c);
609 /* Start from bit x=(p+1)-d in mantissa C */
610 mpfr_prec_t x = p+1-d;
611 mp_size_t kx = n-1 - (x/GMP_NUMB_BITS);
612 mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
613 MPFR_ASSERTD(p > d);
614 DEBUG (printf ("(pre) x=%lu Kx=%ld Sx=%lu\n",
615 (unsigned long) x, (long) kx,
616 (unsigned long) sx));
617 bbcp = (tp[kx] & (MPFR_LIMB_ONE<<sx)) ;
618 /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
619 /* If Cp+1=0, since C'p+1!=0, C'p+2=1 ! */
620 if (MPFR_LIKELY(bbcp==0 || (tp[kx]&MPFR_LIMB_MASK(sx))))
621 bbcp1 = 1;
622 else
623 {
624 /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
625 do {
626 kx--;
627 } while (kx>=0 && tp[kx]==0);
628 bbcp1 = (kx>=0);
629 DEBUG (printf ("(Pre) Scan done for %ld\n", (long) kx));
630 }
631 } /*End of Bcp1 != 0*/
632 DEBUG( printf("(Pre) Cp+1=%d C'p+2=%d\n", bbcp!=0, bbcp1!=0) );
633 } /* End of "can lose a bit" */
634
635 /* Clean shifted C' */
636 mask = ~MPFR_LIMB_MASK (sh);
637 cp[0] &= mask;
638
639 /* Substract the mantissa c from b in a */
640 ap = MPFR_MANT(a);
641 mpn_sub_n (ap, bp, cp, n);
642 DEBUG( mpfr_print_mant_binary("Sub= ", ap, p) );
643
644 /* Normalize: we lose at max one bit*/
645 if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
646 {
647 /* High bit is not set and we have to fix it! */
648 /* Ap >= 010000xxx001 */
649 mpn_lshift(ap, ap, n, 1);
650 /* Ap >= 100000xxx010 */
651 if (MPFR_UNLIKELY(bcp!=0)) /* Check if Cp = -1 */
652 /* Since Cp == -1, we have to substract one more */
653 {
654 mpn_sub_1(ap, ap, n, MPFR_LIMB_ONE<<sh);
655 MPFR_ASSERTD(MPFR_LIMB_MSB(ap[n-1]) != 0);
656 }
657 /* Ap >= 10000xxx001 */
658 /* Final exponent -1 since we have shifted the mantissa */
659 bx--;
660 /* Update bcp and bcp1 */
661 MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
662 MPFR_ASSERTN(bbcp1 != (mp_limb_t) -1);
663 bcp = bbcp;
664 bcp1 = bbcp1;
665 /* We dont't have anymore a valid Cp+1!
666 But since Ap >= 100000xxx001, the final sub can't unnormalize!*/
667 }
668 MPFR_ASSERTD( !(ap[0] & ~mask) );
669
670 /* Rounding */
671 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
672 {
673 if (MPFR_LIKELY(bcp==0))
674 goto truncate;
675 else if ((bcp1) || ((ap[0] & (MPFR_LIMB_ONE<<sh)) != 0))
676 goto sub_one_ulp;
677 else
678 goto truncate;
679 }
680
681 /* Update rounding mode */
682 MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
683 if (rnd_mode == MPFR_RNDZ && (MPFR_LIKELY(bcp || bcp1)))
684 goto sub_one_ulp;
685 goto truncate;
686 }
687 MPFR_RET_NEVER_GO_HERE ();
688
689 /* Sub one ulp to the result */
690 sub_one_ulp:
691 mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
692 /* Result should be smaller than exact value: inexact=-1 */
693 inexact = -1;
694 /* Check normalisation */
695 if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
696 {
697 /* ap was a power of 2, and we lose a bit */
698 /* Now it is 0111111111111111111[00000 */
699 mpn_lshift(ap, ap, n, 1);
700 bx--;
701 /* And the lost bit x depends on Cp+1, and Cp */
702 /* Compute Cp+1 if it isn't already compute (ie d==1) */
703 /* FIXME: Is this case possible? */
704 if (MPFR_UNLIKELY(d == 1))
705 bbcp = 0;
706 DEBUG( printf("(SubOneUlp)Cp=%d, Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0));
707 /* Compute the last bit (Since we have shifted the mantissa)
708 we need one more bit!*/
709 MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
710 if ( (rnd_mode == MPFR_RNDZ && bcp==0)
711 || (rnd_mode==MPFR_RNDN && bbcp==0)
712 || (bcp && bcp1==0) ) /*Exact result*/
713 {
714 ap[0] |= MPFR_LIMB_ONE<<sh;
715 if (rnd_mode == MPFR_RNDN)
716 inexact = 1;
717 DEBUG( printf("(SubOneUlp) Last bit set\n") );
718 }
719 /* Result could be exact if C'p+1 = 0 and rnd == Zero
720 since we have had one more bit to the result */
721 /* Fixme: rnd_mode == MPFR_RNDZ needed ? */
722 if (bcp1==0 && rnd_mode==MPFR_RNDZ)
723 {
724 DEBUG( printf("(SubOneUlp) Exact result\n") );
725 inexact = 0;
726 }
727 }
728
729 goto end_of_sub;
730
731 truncate:
732 /* Check if the result is an exact power of 2: 100000000000
733 in which cases, we could have to do sub_one_ulp due to some nasty reasons:
734 If Result is a Power of 2:
735 + If rnd = AWAY,
736 | If Cp=-1 and C'p+1 = 0, SubOneUlp and the result is EXACT.
737 If Cp=-1 and C'p+1 =-1, SubOneUlp and the result is above.
738 Otherwise truncate
739 + If rnd = NEAREST,
740 If Cp= 0 and Cp+1 =-1 and C'p+2=-1, SubOneUlp and the result is above
741 If cp=-1 and C'p+1 = 0, SubOneUlp and the result is exact.
742 Otherwise truncate.
743 X bit should always be set if SubOneUlp*/
744 if (MPFR_UNLIKELY(ap[n-1] == MPFR_LIMB_HIGHBIT))
745 {
746 mp_size_t k = n-1;
747 do {
748 k--;
749 } while (k>=0 && ap[k]==0);
750 if (MPFR_UNLIKELY(k<0))
751 {
752 /* It is a power of 2! */
753 /* Compute Cp+1 if it isn't already compute (ie d==1) */
754 /* FIXME: Is this case possible? */
755 if (d == 1)
756 bbcp=0;
757 DEBUG( printf("(Truncate) Cp=%d, Cp+1=%d C'p+1=%d C'p+2=%d\n", \
758 bcp!=0, bbcp!=0, bcp1!=0, bbcp1!=0) );
759 MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
760 MPFR_ASSERTN((rnd_mode != MPFR_RNDN) || (bcp != 0) || (bbcp == 0) || (bbcp1 != (mp_limb_t) -1));
761 if (((rnd_mode != MPFR_RNDZ) && bcp)
762 ||
763 ((rnd_mode == MPFR_RNDN) && (bcp == 0) && (bbcp) && (bbcp1)))
764 {
765 DEBUG( printf("(Truncate) Do sub\n") );
766 mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
767 mpn_lshift(ap, ap, n, 1);
768 ap[0] |= MPFR_LIMB_ONE<<sh;
769 bx--;
770 /* FIXME: Explain why it works (or why not)... */
771 inexact = (bcp1 == 0) ? 0 : (rnd_mode==MPFR_RNDN) ? -1 : 1;
772 goto end_of_sub;
773 }
774 }
775 }
776
777 /* Calcul of Inexact flag.*/
778 inexact = MPFR_LIKELY(bcp || bcp1) ? 1 : 0;
779
780 end_of_sub:
781 /* Update Expo */
782 /* FIXME: Is this test really useful?
783 If d==0 : Exact case. This is never called.
784 if 1 < d < p : bx=MPFR_EXP(b) or MPFR_EXP(b)-1 > MPFR_EXP(c) > emin
785 if d == 1 : bx=MPFR_EXP(b). If we could lose any bits, the exact
786 normalisation is called.
787 if d >= p : bx=MPFR_EXP(b) >= MPFR_EXP(c) + p > emin
788 After SubOneUlp, we could have one bit less.
789 if 1 < d < p : bx >= MPFR_EXP(b)-2 >= MPFR_EXP(c) > emin
790 if d == 1 : bx >= MPFR_EXP(b)-1 = MPFR_EXP(c) > emin.
791 if d >= p : bx >= MPFR_EXP(b)-1 > emin since p>=2.
792 */
793 MPFR_ASSERTD( bx >= __gmpfr_emin);
794 /*
795 if (MPFR_UNLIKELY(bx < __gmpfr_emin))
796 {
797 DEBUG( printf("(Final Underflow)\n") );
798 if (rnd_mode == MPFR_RNDN &&
799 (bx < __gmpfr_emin - 1 ||
800 (inexact >= 0 && mpfr_powerof2_raw (a))))
801 rnd_mode = MPFR_RNDZ;
802 MPFR_TMP_FREE(marker);
803 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
804 }
805 */
806 MPFR_SET_EXP (a, bx);
807
808 MPFR_TMP_FREE(marker);
809 MPFR_RET (inexact * MPFR_INT_SIGN (a));
810 }
811