mpz_nextprime
Niels Möller
nisse at lysator.liu.se
Wed Sep 4 04:44:11 UTC 2019
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.
More information about the gmp-devel
mailing list