mpz_nextprime

[Niels Möller]

Seth Troisi <braintwo at gmail.com> writes:

> The idea is to try and avoid performing modulo by storing the distance to
> the
> next multiple of each prime. This avoids divisions, unfortunately it
> requires
> writing back to the table each time a multiple is passed.

Have you tried divisibility test based on 2-adic inverse? It works like
this:

Let p be an odd prime, and x any number that fits in a uint16_t (say),
and you want to know if p divides x.

Precompute 

  p_inv = p^-1 (mod 2^16) 

and 

  p_limit = floor (2^16 / p) = floor ( (2^16 - 1) / p)

Then p divides x if and only if x * p_inv mod 2^16 <= p_limit.

The reason this works is that (i) multiplication by p_inv mod 2^16 is a
one-to-one mapping, and (ii) whenever p divides x, it agrees with normal
division, p_inv * x mod 2^16 = x / p.

Therefore all the multiples of p are mapped to the interval 0 <= x *
p_inv mod 2^16 <= p_limit, and non-multiples are mapped to larger values.

So if you tabulate p_inv and p_limit together with the prime table, 

  r = (moduli[i] + incr) % prime;
  if (r == 0) goto next;

could be replaced with

  if ((uint16_t)((moduli[i] + incr) * p_inv) <= p_limit)
    goto next;

which a multiply instead of a division. I would expect it to be faster
than the loop

  while (incr > distance[i])
    {
      distance[i] += prime;
    }

even though it's a bit unclear to me how many times it typically runs
and how well predictable the branch is.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.