fast inversion
Vincent Lefevre
vincent at vinc17.net
Wed Apr 29 08:07:22 UTC 2015
On 2015-04-29 07:01:39 +0200, bodrato at mail.dm.unipi.it wrote:
> diff -r 6e11cd70e19e gmp-h.in
> --- a/gmp-h.in Mon Apr 27 22:46:53 2015 +0200
> +++ b/gmp-h.in Tue Apr 28 22:41:52 2015 +0200
> @@ -2191,13 +2191,10 @@
> mp_limb_t
> mpn_neg (mp_ptr __gmp_rp, mp_srcptr __gmp_up, mp_size_t __gmp_n)
> {
> - mp_limb_t __gmp_ul, __gmp_cy;
> - __gmp_cy = 0;
> - do {
> - __gmp_ul = *__gmp_up++;
> - *__gmp_rp++ = -__gmp_ul - __gmp_cy;
> - __gmp_cy |= __gmp_ul != 0;
> - } while (--__gmp_n != 0);
> + mp_limb_t __gmp_cy;
> + mpn_com (__gmp_rp, __gmp_up, __gmp_n);
> + __gmp_cy = 1;
> + __gmp_cy -= mpn_add_1 (__gmp_rp, __gmp_rp, __gmp_n, __gmp_cy);
> return __gmp_cy;
> }
> #endif
>
> It would slow down mpn_neg for small values of n...
and possibly for large values of n since mpn_add_1 may do n operations
in the worst case. IMHO, you should do the "carry propagation" first
(basically, loop while the least significant limbs are 0), then call
mpn_com if there are remaining limbs. In my rewrite of mpfr_sum in MPFR,
I did something similar to implement com(x) - 1, which is not provided
by GMP:
mp_limb_t corr2, c;
mp_size_t i = 1;
/* We want to compute com(x) - 1, but GMP doesn't have an
operation for that. The fact is that a sequence of low
significant bits 1 is invariant. Starting at the first
low significant bit 0, we can do the complement with
mpn_com. */
corr2 = MPFR_LIMB_ONE << sd;
c = ~ (sump[0] | MPFR_LIMB_MASK (sd));
sump[0] = c - corr2;
if (c == 0)
{
while (MPFR_ASSERTD (i < sn), sump[i] == MPFR_LIMB_MAX)
i++;
sump[i] = (~ sump[i]) - 1;
i++;
}
if (i < sn)
mpn_com (sump + i, sump + i, sn - i);
--
Vincent Lefèvre <vincent at vinc17.net> - Web: <https://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <https://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)
More information about the gmp-devel
mailing list