1
/* mpfr_sub1 -- internal function to perform a "real" subtraction
3
Copyright 2001 Free Software Foundation.
4
Contributed by the Spaces project, INRIA Lorraine.
6
This file is part of the MPFR Library.
8
The MPFR Library is free software; you can redistribute it and/or modify
9
it under the terms of the GNU Lesser General Public License as published by
10
the Free Software Foundation; either version 2.1 of the License, or (at your
11
option) any later version.
13
The MPFR Library is distributed in the hope that it will be useful, but
14
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16
License for more details.
18
You should have received a copy of the GNU Lesser General Public License
19
along with the MPFR Library; see the file COPYING.LIB. If not, write to
20
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21
MA 02111-1307, USA. */
26
#include "mpfr-impl.h"
28
/* compute sign(b) * (|b| - |c|), with |b| > |c|, diff_exp = EXP(b) - EXP(c)
29
Returns 0 iff result is exact,
30
a negative value when the result is less than the exact value,
31
a positive value otherwise.
35
mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode,
39
mp_exp_unsigned_t diff_exp;
40
mp_prec_t cancel, cancel1;
41
mp_size_t cancel2, an, bn, cn, cn0;
42
mp_limb_t *ap, *bp, *cp;
43
mp_limb_t carry, bb, cc, borrow = 0;
44
int inexact = 0, shift_b, shift_c, is_exact = 1, down = 0, add_exp = 0;
50
an = 1 + (MPFR_PREC(a) - 1) / BITS_PER_MP_LIMB;
52
sign = mpfr_cmp2 (b, c, &cancel);
55
if (rnd_mode == GMP_RNDD)
63
/* If subtraction: sign(a) = sign * sign(b) */
64
if (sub && MPFR_SIGN(a) != sign * MPFR_SIGN(b))
67
if (sign < 0) /* swap b and c so that |b| > |c| */
73
/* If addition: sign(a) = sign of the larger argument in absolute value */
75
MPFR_SET_SAME_SIGN(a, b);
77
diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c);
79
/* reserve a space to store b aligned with the result, i.e. shifted by
80
(-cancel) % BITS_PER_MP_LIMB to the right */
81
bn = 1 + (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
82
shift_b = cancel % BITS_PER_MP_LIMB;
84
shift_b = BITS_PER_MP_LIMB - shift_b;
85
cancel1 = (cancel + shift_b) / BITS_PER_MP_LIMB;
86
/* the high cancel1 limbs from b should not be taken into account */
88
bp = MPFR_MANT(b); /* no need of an extra space */
91
bp = TMP_ALLOC ((bn + 1) * BYTES_PER_MP_LIMB);
92
bp[0] = mpn_rshift (bp + 1, MPFR_MANT(b), bn++, shift_b);
95
/* reserve a space to store c aligned with the result, i.e. shifted by
96
(diff_exp-cancel) % BITS_PER_MP_LIMB to the right */
97
cn = 1 + (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
98
shift_c = diff_exp - (cancel % BITS_PER_MP_LIMB);
99
shift_c = (shift_c + BITS_PER_MP_LIMB) % BITS_PER_MP_LIMB;
104
cp = TMP_ALLOC ((cn + 1) * BYTES_PER_MP_LIMB);
105
cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c);
109
printf("shift_b=%u shift_c=%u\n", shift_b, shift_c);
112
/* ensure ap != bp and ap != cp */
115
bp = (mp_ptr) TMP_ALLOC(bn * BYTES_PER_MP_LIMB);
116
MPN_COPY (bp, ap, bn);
117
/* ap == cp cannot occur since we would have b=c, which is detected
118
in mpfr_add or mpfr_sub */
122
cp = (mp_ptr) TMP_ALLOC (cn * BYTES_PER_MP_LIMB);
123
MPN_COPY(cp, ap, cn);
126
/* here we have shift_c = (diff_exp - cancel) % BITS_PER_MP_LIMB,
127
thus we want cancel2 = ceil((cancel - diff_exp) / BITS_PER_MP_LIMB) */
129
cancel2 = (long int) (cancel - (diff_exp - shift_c)) / BITS_PER_MP_LIMB;
130
/* the high cancel2 limbs from b should not be taken into account */
132
printf("cancel=%u cancel1=%u cancel2=%d\n", cancel, cancel1, cancel2);
136
<----------------+-----------|---->
137
<----------PREC(a)----------><-sh->
139
limbs bp[bn-cancel1-1]
140
<--...-----><----------------+-----------+----------->
142
limbs cp[cn-cancel2-1] cancel2 >= 0
143
<--...--><----------------+----------------+---------------->
144
(-cancel2) cancel2 < 0
145
limbs <----------------+---------------->
148
/* first part: put in ap[0..an-1] the value of high(b) - high(c),
149
where high(b) consists of the high an+cancel1 limbs of b,
150
and high(c) consists of the high an+cancel2 limbs of c.
153
/* copy high(b) into a */
154
if (an + cancel1 <= bn) /* a: <----------------+-----------|---->
155
b: <-----------------------------------------> */
156
MPN_COPY (ap, bp + bn - (an + cancel1), an);
157
else /* a: <----------------+-----------|---->
158
b: <-------------------------> */
159
if (cancel1 < bn) /* otherwise b does not overlap with a */
161
MPN_ZERO (ap, an + cancel1 - bn);
162
MPN_COPY (ap + an + cancel1 - bn, bp, bn - cancel1);
168
printf("after copying high(b), a="); mpfr_print_binary(a); putchar('\n');
171
/* subtract high(c) */
172
if (an + cancel2 > 0) /* otherwise c does not overlap with a */
178
if (an + cancel2 <= cn) /* a: <----------------------------->
179
c: <-----------------------------------------> */
180
mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an);
181
else /* a: <---------------------------->
182
c: <-------------------------> */
184
ap2 = ap + an + cancel2 - cn;
186
mpn_sub_n (ap2, ap2, cp, cn - cancel2);
189
else /* cancel2 < 0 */
191
if (an + cancel2 <= cn) /* a: <----------------------------->
192
c: <-----------------------------> */
193
borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an + cancel2);
194
else /* a: <---------------------------->
195
c: <----------------> */
197
ap2 = ap + an + cancel2 - cn;
198
borrow = mpn_sub_n (ap2, ap2, cp, cn);
200
ap2 = ap + an + cancel2;
201
mpn_sub_1 (ap2, ap2, -cancel2, borrow);
206
printf("after subtracting high(c), a="); mpfr_print_binary(a); putchar('\n');
209
/* now perform rounding */
210
sh = an * BITS_PER_MP_LIMB - MPFR_PREC(a); /* last unused bits from a */
211
carry = ap[0] & ((MP_LIMB_T_ONE << sh) - MP_LIMB_T_ONE);
214
if (rnd_mode == GMP_RNDN)
218
is_exact = (carry == 0);
219
/* can decide except when carry = 2^(sh-1) [middle]
220
or carry = 0 [truncate, but cannot decide inexact flag] */
221
down = (carry < (MP_LIMB_T_ONE << (sh - 1)));
222
if (carry > (MP_LIMB_T_ONE << (sh - 1)))
224
else if ((0 < carry) && down)
226
inexact = -1; /* result if smaller than exact value */
231
else /* directed rounding: set rnd_mode to RNDZ iff towards zero */
233
if (((rnd_mode == GMP_RNDD) && (MPFR_SIGN(a) > 0)) ||
234
((rnd_mode == GMP_RNDU) && (MPFR_SIGN(a) < 0)))
239
if (rnd_mode == GMP_RNDZ)
244
else /* round away */
249
/* we have to consider the low (bn - (an+cancel1)) limbs from b,
250
and the (cn - (an+cancel2)) limbs from c. */
253
cn -= (long int) an + cancel2;
255
printf("last %d bits from a are %lu, bn=%ld, cn=%ld\n", sh, carry, bn, cn);
258
for (k = 0; (bn > 0) || (cn > 0); k = 1)
260
bb = (bn > 0) ? bp[--bn] : 0;
261
if ((cn > 0) && (cn-- <= cn0))
269
if ((rnd_mode == GMP_RNDN) && !k && sh == 0)
271
mp_limb_t half = GMP_LIMB_HIGHBIT;
273
is_exact = (bb == cc);
275
/* add one ulp if bb > cc + half
276
truncate if cc - half < bb < cc + half
277
sub one ulp if bb < cc - half
297
printf(" bb=%lu cc=%lu down=%d is_exact=%d\n", bb, cc, down, is_exact);
301
if (rnd_mode == GMP_RNDZ)
303
else if (rnd_mode != GMP_RNDN) /* round away */
308
else /* round to nearest */
310
if (is_exact && sh == 0)
315
else if (down && sh == 0)
319
inexact = (is_exact) ? 1 : -1;
326
if (rnd_mode == GMP_RNDZ)
331
else if (rnd_mode != GMP_RNDN) /* round away */
333
else /* round to nearest */
351
if ((rnd_mode == GMP_RNDN) && !is_exact)
353
/* even rounding rule */
354
if ((ap[0] >> sh) & 1)
356
if (down) goto sub_one_ulp;
357
else goto add_one_ulp;
360
inexact = (down) ? 1 : -1;
366
sub_one_ulp: /* add one unit in last place to a */
367
mpn_sub_1 (ap, ap, an, MP_LIMB_T_ONE << sh);
371
add_one_ulp: /* add one unit in last place to a */
372
if (mpn_add_1 (ap, ap, an, MP_LIMB_T_ONE << sh)) /* result is a power of 2 */
374
ap[an-1] = GMP_LIMB_HIGHBIT;
377
inexact = 1; /* result larger than exact value */
380
if ((ap[an-1] >> (BITS_PER_MP_LIMB - 1)) == 0) /* case 1 - epsilon */
382
ap[an-1] = GMP_LIMB_HIGHBIT;
387
/* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking
388
care of underflows/overflows in that computation, and of the allowed
394
cancel -= add_exp; /* still valid as unsigned long */
395
exp_b = MPFR_EXP(b); /* save it in case a equals b */
396
MPFR_EXP(a) = MPFR_EXP(b) - cancel;
397
if ((MPFR_EXP(a) > exp_b) /* underflow in type mp_exp_t */
398
|| (MPFR_EXP(a) < __mpfr_emin))
401
return mpfr_set_underflow (a, rnd_mode, MPFR_SIGN(a));
404
else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
406
/* in case cancel = 0, add_exp can still be 1, in case b is just
407
below a power of two, c is very small, prec(a) < prec(b),
408
and rnd=away or nearest */
409
if (add_exp && MPFR_EXP(b) == __mpfr_emax)
412
return mpfr_set_overflow (a, rnd_mode, MPFR_SIGN(a));
414
MPFR_EXP(a) = MPFR_EXP(b) + add_exp;
418
printf ("result is a="); mpfr_print_binary(a); putchar('\n');
420
/* check that result is msb-normalized */
421
MPFR_ASSERTN(ap[an-1] > ~ap[an-1]);
422
return inexact * MPFR_SIGN(a);
added
removed
src/gmp/mpfr/sub1.c
