1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
|
/* mpfr_sub1 -- internal function to perform a "real" subtraction
Copyright 2001-2015 Free Software Foundation, Inc.
Contributed by the AriC and Caramel projects, INRIA.
This file is part of the GNU MPFR Library.
The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
The GNU MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
#include "mpfr-impl.h"
/* compute sign(b) * (|b| - |c|), with |b| > |c|, diff_exp = EXP(b) - EXP(c)
Returns 0 iff result is exact,
a negative value when the result is less than the exact value,
a positive value otherwise.
*/
int
mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
int sign;
mpfr_uexp_t diff_exp;
mpfr_prec_t cancel, cancel1;
mp_size_t cancel2, an, bn, cn, cn0;
mp_limb_t *ap, *bp, *cp;
mp_limb_t carry, bb, cc;
int inexact, shift_b, shift_c, add_exp = 0;
int cmp_low = 0; /* used for rounding to nearest: 0 if low(b) = low(c),
negative if low(b) < low(c), positive if low(b)>low(c) */
int sh, k;
MPFR_TMP_DECL(marker);
MPFR_TMP_MARK(marker);
ap = MPFR_MANT(a);
an = MPFR_LIMB_SIZE(a);
sign = mpfr_cmp2 (b, c, &cancel);
if (MPFR_UNLIKELY(sign == 0))
{
if (rnd_mode == MPFR_RNDD)
MPFR_SET_NEG (a);
else
MPFR_SET_POS (a);
MPFR_SET_ZERO (a);
MPFR_RET (0);
}
/*
* If subtraction: sign(a) = sign * sign(b)
* If addition: sign(a) = sign of the larger argument in absolute value.
*
* Both cases can be simplidied in:
* if (sign>0)
* if addition: sign(a) = sign * sign(b) = sign(b)
* if subtraction, b is greater, so sign(a) = sign(b)
* else
* if subtraction, sign(a) = - sign(b)
* if addition, sign(a) = sign(c) (since c is greater)
* But if it is an addition, sign(b) and sign(c) are opposed!
* So sign(a) = - sign(b)
*/
if (sign < 0) /* swap b and c so that |b| > |c| */
{
mpfr_srcptr t;
MPFR_SET_OPPOSITE_SIGN (a,b);
t = b; b = c; c = t;
}
else
MPFR_SET_SAME_SIGN (a,b);
/* Check if c is too small.
A more precise test is to replace 2 by
(rnd == MPFR_RNDN) + mpfr_power2_raw (b)
but it is more expensive and not very useful */
if (MPFR_UNLIKELY (MPFR_GET_EXP (c) <= MPFR_GET_EXP (b)
- (mpfr_exp_t) MAX (MPFR_PREC (a), MPFR_PREC (b)) - 2))
{
/* Remember, we can't have an exact result! */
/* A.AAAAAAAAAAAAAAAAA
= B.BBBBBBBBBBBBBBB
- C.CCCCCCCCCCCCC */
/* A = S*ABS(B) +/- ulp(a) */
MPFR_SET_EXP (a, MPFR_GET_EXP (b));
MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), MPFR_PREC (b),
rnd_mode, MPFR_SIGN (a),
if (MPFR_UNLIKELY ( ++MPFR_EXP (a) > __gmpfr_emax))
inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)));
/* inexact = mpfr_set4 (a, b, rnd_mode, MPFR_SIGN (a)); */
if (inexact == 0)
{
/* a = b (Exact)
But we know it isn't (Since we have to remove `c')
So if we round to Zero, we have to remove one ulp.
Otherwise the result is correctly rounded. */
if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
{
mpfr_nexttozero (a);
MPFR_RET (- MPFR_INT_SIGN (a));
}
MPFR_RET (MPFR_INT_SIGN (a));
}
else
{
/* A.AAAAAAAAAAAAAA
= B.BBBBBBBBBBBBBBB
- C.CCCCCCCCCCCCC */
/* It isn't exact so Prec(b) > Prec(a) and the last
Prec(b)-Prec(a) bits of `b' are not zeros.
Which means that removing c from b can't generate a carry
execpt in case of even rounding.
In all other case the result and the inexact flag should be
correct (We can't have an exact result).
In case of EVEN rounding:
1.BBBBBBBBBBBBBx10
- 1.CCCCCCCCCCCC
= 1.BBBBBBBBBBBBBx01 Rounded to Prec(b)
= 1.BBBBBBBBBBBBBx Nearest / Rounded to Prec(a)
Set gives:
1.BBBBBBBBBBBBB0 if inexact == EVEN_INEX (x == 0)
1.BBBBBBBBBBBBB1+1 if inexact == -EVEN_INEX (x == 1)
which means we get a wrong rounded result if x==1,
i.e. inexact= MPFR_EVEN_INEX */
if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX*MPFR_INT_SIGN (a)))
{
mpfr_nexttozero (a);
inexact = -MPFR_INT_SIGN (a);
}
MPFR_RET (inexact);
}
}
diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
/* reserve a space to store b aligned with the result, i.e. shifted by
(-cancel) % GMP_NUMB_BITS to the right */
bn = MPFR_LIMB_SIZE (b);
MPFR_UNSIGNED_MINUS_MODULO (shift_b, cancel);
cancel1 = (cancel + shift_b) / GMP_NUMB_BITS;
/* the high cancel1 limbs from b should not be taken into account */
if (MPFR_UNLIKELY (shift_b == 0))
{
bp = MPFR_MANT(b); /* no need of an extra space */
/* Ensure ap != bp */
if (MPFR_UNLIKELY (ap == bp))
{
bp = MPFR_TMP_LIMBS_ALLOC (bn);
MPN_COPY (bp, ap, bn);
}
}
else
{
bp = MPFR_TMP_LIMBS_ALLOC (bn + 1);
bp[0] = mpn_rshift (bp + 1, MPFR_MANT(b), bn++, shift_b);
}
/* reserve a space to store c aligned with the result, i.e. shifted by
(diff_exp-cancel) % GMP_NUMB_BITS to the right */
cn = MPFR_LIMB_SIZE(c);
if ((UINT_MAX % GMP_NUMB_BITS) == (GMP_NUMB_BITS-1)
&& ((-(unsigned) 1)%GMP_NUMB_BITS > 0))
shift_c = ((mpfr_uexp_t) diff_exp - cancel) % GMP_NUMB_BITS;
else
{
shift_c = diff_exp - (cancel % GMP_NUMB_BITS);
shift_c = (shift_c + GMP_NUMB_BITS) % GMP_NUMB_BITS;
}
MPFR_ASSERTD( shift_c >= 0 && shift_c < GMP_NUMB_BITS);
if (MPFR_UNLIKELY(shift_c == 0))
{
cp = MPFR_MANT(c);
/* Ensure ap != cp */
if (ap == cp)
{
cp = MPFR_TMP_LIMBS_ALLOC (cn);
MPN_COPY(cp, ap, cn);
}
}
else
{
cp = MPFR_TMP_LIMBS_ALLOC (cn + 1);
cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c);
}
#ifdef DEBUG
printf ("rnd=%s shift_b=%d shift_c=%d diffexp=%lu\n",
mpfr_print_rnd_mode (rnd_mode), shift_b, shift_c,
(unsigned long) diff_exp);
#endif
MPFR_ASSERTD (ap != cp);
MPFR_ASSERTD (bp != cp);
/* here we have shift_c = (diff_exp - cancel) % GMP_NUMB_BITS,
0 <= shift_c < GMP_NUMB_BITS
thus we want cancel2 = ceil((cancel - diff_exp) / GMP_NUMB_BITS) */
/* Possible optimization with a C99 compiler (i.e. well-defined
integer division): if MPFR_PREC_MAX is reduced to
((mpfr_prec_t)((mpfr_uprec_t)(~(mpfr_uprec_t)0)>>1) - GMP_NUMB_BITS + 1)
and diff_exp is of type mpfr_exp_t (no need for mpfr_uexp_t, since
the sum or difference of 2 exponents must be representable, as used
by the multiplication code), then the computation of cancel2 could
be simplified to
cancel2 = (cancel - (diff_exp - shift_c)) / GMP_NUMB_BITS;
because cancel, diff_exp and shift_c are all non-negative and
these variables are signed. */
MPFR_ASSERTD (cancel >= 0);
if (cancel >= diff_exp)
/* Note that cancel is signed and will be converted to mpfr_uexp_t
(type of diff_exp) in the expression below, so that this will
work even if cancel is very large and diff_exp = 0. */
cancel2 = (cancel - diff_exp + (GMP_NUMB_BITS - 1)) / GMP_NUMB_BITS;
else
cancel2 = - (mp_size_t) ((diff_exp - cancel) / GMP_NUMB_BITS);
/* the high cancel2 limbs from b should not be taken into account */
#ifdef DEBUG
printf ("cancel=%lu cancel1=%lu cancel2=%ld\n",
(unsigned long) cancel, (unsigned long) cancel1, (long) cancel2);
#endif
/* ap[an-1] ap[0]
<----------------+-----------|---->
<----------PREC(a)----------><-sh->
cancel1
limbs bp[bn-cancel1-1]
<--...-----><----------------+-----------+----------->
cancel2
limbs cp[cn-cancel2-1] cancel2 >= 0
<--...--><----------------+----------------+---------------->
(-cancel2) cancel2 < 0
limbs <----------------+---------------->
*/
/* first part: put in ap[0..an-1] the value of high(b) - high(c),
where high(b) consists of the high an+cancel1 limbs of b,
and high(c) consists of the high an+cancel2 limbs of c.
*/
/* copy high(b) into a */
if (MPFR_LIKELY(an + (mp_size_t) cancel1 <= bn))
/* a: <----------------+-----------|---->
b: <-----------------------------------------> */
MPN_COPY (ap, bp + bn - (an + cancel1), an);
else
/* a: <----------------+-----------|---->
b: <-------------------------> */
if ((mp_size_t) cancel1 < bn) /* otherwise b does not overlap with a */
{
MPN_ZERO (ap, an + cancel1 - bn);
MPN_COPY (ap + (an + cancel1 - bn), bp, bn - cancel1);
}
else
MPN_ZERO (ap, an);
#ifdef DEBUG
printf("after copying high(b), a="); mpfr_print_binary(a); putchar('\n');
#endif
/* subtract high(c) */
if (MPFR_LIKELY(an + cancel2 > 0)) /* otherwise c does not overlap with a */
{
mp_limb_t *ap2;
if (cancel2 >= 0)
{
if (an + cancel2 <= cn)
/* a: <----------------------------->
c: <-----------------------------------------> */
mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an);
else
/* a: <---------------------------->
c: <-------------------------> */
{
ap2 = ap + an + (cancel2 - cn);
if (cn > cancel2)
mpn_sub_n (ap2, ap2, cp, cn - cancel2);
}
}
else /* cancel2 < 0 */
{
mp_limb_t borrow;
if (an + cancel2 <= cn)
/* a: <----------------------------->
c: <-----------------------------> */
borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2),
an + cancel2);
else
/* a: <---------------------------->
c: <----------------> */
{
ap2 = ap + an + cancel2 - cn;
borrow = mpn_sub_n (ap2, ap2, cp, cn);
}
ap2 = ap + an + cancel2;
mpn_sub_1 (ap2, ap2, -cancel2, borrow);
}
}
#ifdef DEBUG
printf("after subtracting high(c), a=");
mpfr_print_binary(a);
putchar('\n');
#endif
/* now perform rounding */
sh = (mpfr_prec_t) an * GMP_NUMB_BITS - MPFR_PREC(a);
/* last unused bits from a */
carry = ap[0] & MPFR_LIMB_MASK (sh);
ap[0] -= carry;
if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
{
if (MPFR_LIKELY(sh))
{
/* can decide except when carry = 2^(sh-1) [middle]
or carry = 0 [truncate, but cannot decide inexact flag] */
if (carry > (MPFR_LIMB_ONE << (sh - 1)))
goto add_one_ulp;
else if ((0 < carry) && (carry < (MPFR_LIMB_ONE << (sh - 1))))
{
inexact = -1; /* result if smaller than exact value */
goto truncate;
}
/* now carry = 2^(sh-1), in which case cmp_low=2,
or carry = 0, in which case cmp_low=0 */
cmp_low = (carry == 0) ? 0 : 2;
}
}
else /* directed rounding: set rnd_mode to RNDZ iff toward zero */
{
if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, MPFR_IS_NEG(a)))
rnd_mode = MPFR_RNDZ;
if (carry)
{
if (rnd_mode == MPFR_RNDZ)
{
inexact = -1;
goto truncate;
}
else /* round away */
goto add_one_ulp;
}
}
/* we have to consider the low (bn - (an+cancel1)) limbs from b,
and the (cn - (an+cancel2)) limbs from c. */
bn -= an + cancel1;
cn0 = cn;
cn -= an + cancel2;
#ifdef DEBUG
printf ("last sh=%d bits from a are %lu, bn=%ld, cn=%ld\n",
sh, (unsigned long) carry, (long) bn, (long) cn);
#endif
/* for rounding to nearest, we couldn't conclude up to here in the following
cases:
1. sh = 0, then cmp_low=0: we can either truncate, subtract one ulp
or add one ulp: -1 ulp < low(b)-low(c) < 1 ulp
2. sh > 0 but the low sh bits from high(b)-high(c) equal 2^(sh-1):
-0.5 ulp <= -1/2^sh < low(b)-low(c)-0.5 < 1/2^sh <= 0.5 ulp
we can't decide the rounding, in that case cmp_low=2:
either we truncate and flag=-1, or we add one ulp and flag=1
3. the low sh>0 bits from high(b)-high(c) equal 0: we know we have to
truncate but we can't decide the ternary value, here cmp_low=0:
-0.5 ulp <= -1/2^sh < low(b)-low(c) < 1/2^sh <= 0.5 ulp
we always truncate and inexact can be any of -1,0,1
*/
/* note: here cn might exceed cn0, in which case we consider a zero limb */
for (k = 0; (bn > 0) || (cn > 0); k = 1)
{
/* if cmp_low < 0, we know low(b) - low(c) < 0
if cmp_low > 0, we know low(b) - low(c) > 0
(more precisely if cmp_low = 2, low(b) - low(c) = 0.5 ulp so far)
if cmp_low = 0, so far low(b) - low(c) = 0 */
/* get next limbs */
bb = (bn > 0) ? bp[--bn] : 0;
if ((cn > 0) && (cn-- <= cn0))
cc = cp[cn];
else
cc = 0;
/* cmp_low compares low(b) and low(c) */
if (cmp_low == 0) /* case 1 or 3 */
cmp_low = (bb < cc) ? -2+k : (bb > cc) ? 1 : 0;
/* Case 1 for k=0 splits into 7 subcases:
1a: bb > cc + half
1b: bb = cc + half
1c: 0 < bb - cc < half
1d: bb = cc
1e: -half < bb - cc < 0
1f: bb - cc = -half
1g: bb - cc < -half
Case 2 splits into 3 subcases:
2a: bb > cc
2b: bb = cc
2c: bb < cc
Case 3 splits into 3 subcases:
3a: bb > cc
3b: bb = cc
3c: bb < cc
*/
/* the case rounding to nearest with sh=0 is special since one couldn't
subtract above 1/2 ulp in the trailing limb of the result */
if (rnd_mode == MPFR_RNDN && sh == 0 && k == 0) /* case 1 for k=0 */
{
mp_limb_t half = MPFR_LIMB_HIGHBIT;
/* add one ulp if bb > cc + half
truncate if cc - half < bb < cc + half
sub one ulp if bb < cc - half
*/
if (cmp_low < 0) /* bb < cc: -1 ulp < low(b) - low(c) < 0,
cases 1e, 1f and 1g */
{
if (cc >= half)
cc -= half;
else /* since bb < cc < half, bb+half < 2*half */
bb += half;
/* now we have bb < cc + half:
we have to subtract one ulp if bb < cc,
and truncate if bb > cc */
}
else if (cmp_low >= 0) /* bb >= cc, cases 1a to 1d */
{
if (cc < half)
cc += half;
else /* since bb >= cc >= half, bb - half >= 0 */
bb -= half;
/* now we have bb > cc - half: we have to add one ulp if bb > cc,
and truncate if bb < cc */
if (cmp_low > 0)
cmp_low = 2;
}
}
#ifdef DEBUG
printf ("k=%u bb=%lu cc=%lu cmp_low=%d\n", k,
(unsigned long) bb, (unsigned long) cc, cmp_low);
#endif
if (cmp_low < 0) /* low(b) - low(c) < 0: either truncate or subtract
one ulp */
{
if (rnd_mode == MPFR_RNDZ)
goto sub_one_ulp; /* set inexact=-1 */
else if (rnd_mode != MPFR_RNDN) /* round away */
{
inexact = 1;
goto truncate;
}
else /* round to nearest */
{
/* If cmp_low < 0 and bb > cc, then -0.5 ulp < low(b)-low(c) < 0,
whatever the value of sh.
If sh>0, then cmp_low < 0 implies that the initial neglected
sh bits were 0 (otherwise cmp_low=2 initially), thus the
weight of the new bits is less than 0.5 ulp too.
If k > 0 (and sh=0) this means that either the first neglected
limbs bb and cc were equal (thus cmp_low was 0 for k=0),
or we had bb - cc = -0.5 ulp or 0.5 ulp.
The last case is not possible here since we would have
cmp_low > 0 which is sticky.
In the first case (where we have cmp_low = -1), we truncate,
whereas in the 2nd case we have cmp_low = -2 and we subtract
one ulp.
*/
if (bb > cc || sh > 0 || cmp_low == -1)
{ /* -0.5 ulp < low(b)-low(c) < 0,
bb > cc corresponds to cases 1e and 1f1
sh > 0 corresponds to cases 3c and 3b3
cmp_low = -1 corresponds to case 1d3 (also 3b3) */
inexact = 1;
goto truncate;
}
else if (bb < cc) /* here sh = 0 and low(b)-low(c) < -0.5 ulp,
this corresponds to cases 1g and 1f3 */
goto sub_one_ulp;
/* the only case where we can't conclude is sh=0 and bb=cc,
i.e., we have low(b) - low(c) = -0.5 ulp (up to now), thus
we don't know if we must truncate or subtract one ulp.
Note: for sh=0 we can't have low(b) - low(c) = -0.5 ulp up to
now, since low(b) - low(c) > 1/2^sh */
}
}
else if (cmp_low > 0) /* 0 < low(b) - low(c): either truncate or
add one ulp */
{
if (rnd_mode == MPFR_RNDZ)
{
inexact = -1;
goto truncate;
}
else if (rnd_mode != MPFR_RNDN) /* round away */
goto add_one_ulp;
else /* round to nearest */
{
if (bb > cc)
{
/* if sh=0, then bb>cc means that low(b)-low(c) > 0.5 ulp,
and similarly when cmp_low=2 */
if (cmp_low == 2) /* cases 1a, 1b1, 2a and 2b1 */
goto add_one_ulp;
/* sh > 0 and cmp_low > 0: this implies that the sh initial
neglected bits were 0, and the remaining low(b)-low(c)>0,
but its weight is less than 0.5 ulp */
else /* 0 < low(b) - low(c) < 0.5 ulp, this corresponds to
cases 3a, 1d1 and 3b1 */
{
inexact = -1;
goto truncate;
}
}
else if (bb < cc) /* 0 < low(b) - low(c) < 0.5 ulp, cases 1c,
1b3, 2b3 and 2c */
{
inexact = -1;
goto truncate;
}
/* the only case where we can't conclude is bb=cc, i.e.,
low(b) - low(c) = 0.5 ulp (up to now), thus we don't know
if we must truncate or add one ulp. */
}
}
/* after k=0, we cannot conclude in the following cases, we split them
according to the values of bb and cc for k=1:
1b. sh=0 and cmp_low = 1 and bb-cc = half [around 0.5 ulp]
1b1. bb > cc: add one ulp, inex = 1
1b2: bb = cc: cannot conclude
1b3: bb < cc: truncate, inex = -1
1d. sh=0 and cmp_low = 0 and bb-cc = 0 [around 0]
1d1: bb > cc: truncate, inex = -1
1d2: bb = cc: cannot conclude
1d3: bb < cc: truncate, inex = +1
1f. sh=0 and cmp_low = -1 and bb-cc = -half [around -0.5 ulp]
1f1: bb > cc: truncate, inex = +1
1f2: bb = cc: cannot conclude
1f3: bb < cc: sub one ulp, inex = -1
2b. sh > 0 and cmp_low = 2 and bb=cc [around 0.5 ulp]
2b1. bb > cc: add one ulp, inex = 1
2b2: bb = cc: cannot conclude
2b3: bb < cc: truncate, inex = -1
3b. sh > 0 and cmp_low = 0 [around 0]
3b1. bb > cc: truncate, inex = -1
3b2: bb = cc: cannot conclude
3b3: bb < cc: truncate, inex = +1
*/
}
if ((rnd_mode == MPFR_RNDN) && cmp_low != 0)
{
/* even rounding rule */
if ((ap[0] >> sh) & 1)
{
if (cmp_low < 0)
goto sub_one_ulp;
else
goto add_one_ulp;
}
else
inexact = (cmp_low > 0) ? -1 : 1;
}
else
inexact = 0;
goto truncate;
sub_one_ulp: /* sub one unit in last place to a */
mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh);
inexact = -1;
goto end_of_sub;
add_one_ulp: /* add one unit in last place to a */
if (MPFR_UNLIKELY(mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh)))
/* result is a power of 2: 11111111111111 + 1 = 1000000000000000 */
{
ap[an-1] = MPFR_LIMB_HIGHBIT;
add_exp = 1;
}
inexact = 1; /* result larger than exact value */
truncate:
if (MPFR_UNLIKELY((ap[an-1] >> (GMP_NUMB_BITS - 1)) == 0))
/* case 1 - epsilon */
{
ap[an-1] = MPFR_LIMB_HIGHBIT;
add_exp = 1;
}
end_of_sub:
/* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking
care of underflows/overflows in that computation, and of the allowed
exponent range */
if (MPFR_LIKELY(cancel))
{
mpfr_exp_t exp_a;
cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */
exp_a = MPFR_GET_EXP (b) - cancel;
if (MPFR_UNLIKELY(exp_a < __gmpfr_emin))
{
MPFR_TMP_FREE(marker);
if (rnd_mode == MPFR_RNDN &&
(exp_a < __gmpfr_emin - 1 ||
(inexact >= 0 && mpfr_powerof2_raw (a))))
rnd_mode = MPFR_RNDZ;
return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
}
MPFR_SET_EXP (a, exp_a);
}
else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
{
/* in case cancel = 0, add_exp can still be 1, in case b is just
below a power of two, c is very small, prec(a) < prec(b),
and rnd=away or nearest */
mpfr_exp_t exp_b;
exp_b = MPFR_GET_EXP (b);
if (MPFR_UNLIKELY(add_exp && exp_b == __gmpfr_emax))
{
MPFR_TMP_FREE(marker);
return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
}
MPFR_SET_EXP (a, exp_b + add_exp);
}
MPFR_TMP_FREE(marker);
#ifdef DEBUG
printf ("result is a="); mpfr_print_binary(a); putchar('\n');
#endif
/* check that result is msb-normalized */
MPFR_ASSERTD(ap[an-1] > ~ap[an-1]);
MPFR_RET (inexact * MPFR_INT_SIGN (a));
}
|