fast inversion

[Vincent Lefevre]

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)