Better b^e mod m

Torbjorn Granlund tg at
Thu Oct 18 20:40:36 CEST 2012

nisse at (Niels Möller) writes:

  Torbjorn Granlund <tg at> writes:
  > nisse at (Niels Möller) writes:
  >   (I've also been considering if we could also use single-limb hensel
  >   reduction, doing b multiplies as
  >     (x' b) B^{-1}
  >   Then we get a factor of B error (which is going to be amplified while
  >   squaring).
  > I am afraid I don't follow you.  Specifically, what does B^{-1} denote?
  > Do you compute this inverse modulo something, do you intend to form
  > either (x' b)/B or floor((x' b)/B)?
  Everything mod m, division by B (mod m) is one iteration redc_1. So here
  I'm suggesting that multiplying the accumulator x by the single-limb
  base b is done as
    tp[n] = mpn_mul_1 (tp, xp, n, b);
    q = -tp[0] * minv; /* Inverse of m mod B */
    tp[n] += mpn_addmul_1 (tp, mp, n, q); /* FIXME: Handle any carry */ 
    MPN_COPY (xp, tp + 1, n);
  This actually computes x <-- x * (b / B) (mod m). So we need to deal
  with the gratuitous factor of B somewhere else.
  Not sure if this is useful, but in the cases where the extra power of B
  can be handled for free, I'd expect it to be slightly faster than using
  Euclidean reduction.
Should be, yes.

Another alternative might be to simply not reduce after the squaring:

      mpn_sqr (...)
      if (next-exp_bit-set)
        mpn_mul_1 (...)
      mpn_bdiv_r (...)

The intermediate result should be in redc form, the base argument should
be in plain form.


More information about the gmp-devel mailing list