Optimizing Modulus (was Re: Would someone mind elaborating the
explanation in the manual?)
Sisyphus
kalinabears at iinet.net.au
Tue Oct 28 16:54:35 CET 2003
David McKen wrote:
>
> --Begin
> remainder = mpz_fdiv_ui (up_counter, prime_product_3_29);
> if(!((remainder % 3) &&
> (remainder % 5) &&
> (remainder % 7) &&
> (remainder % 11) &&
> (remainder % 13) &&
> (remainder % 17) &&
> (remainder % 19) &&
> (remainder % 23) &&
> (remainder % 29)) )
> {
> update_func (up_counter, up_counter, 2);
> continue;
> }
> --End
>
I'm assuming that up_counter is odd.
Looks to me that what you really should be doing is sieving a range of
numbers (let's say the odd numbers in the range up_counter to
up_counter+10000).
You then find the first number in the range that is divisible by 3, and
cross it off, and cross off every third number after it. Then find the
first number in the range that is divisible by 5 and cross it off, and
cross off every fifth number after it. And so on for 7, 11, ... 29.
By then, the only numbers that haven't been crossed off are the numbers
that are not divisible by [3..29].
You keep track of which numbers have been "crossed off" using a bit
vector (5000 bits long) initialised with all bits set to (say) 1.
"Crossing off" simply involves changing the appropriate bit to 0.
You can then determine which numbers have been crossed off by examining
the bits in the vector.
eg bits 5, 7 and 8 set to 1 would tell you that up_counter+10,
up_counter+14 and up_counter+16 are not divisible by any of [3..29].
(In case I'm not being clear, 'up_counter' is now just the initial value
given - it doesn't get incremented.)
One inefficiency of this method is that many values are getting "crossed
off" more than once - which is wasteful. This can be addressed (by some
sort of 'wheel', I think), but I can't find a link at the moment, and
it's not something I've personally attempted. I'm not sure of the
savings to be made.
I find I can take (eg) a ten thousand decimal digit odd number, x, and
find all of the odd numbers in the range [x .. x+100000] not divisible
by any prime less than 2^20 - in about 2.1 seconds on my 1MHz pc.
That includes the time it takes to generate the primes in the first
place with a Sieve of Eratosthenes. You can vary the size of the range
without greatly affecting the timing. A range of [x .. x+200000] adds
about 0.01 seconds to the time taken. The big efect is, of course, how
deeply you sieve. Sieving (in this example) for all primes less than
2^25 takes 56 seconds. (I'm not making use of the time-saving technique
we started off discussing, btw. Something for me to consider.)
Anyway ... that was just by way of example ... I'm just a wannabe
programmer, so if someone asks "how come it's taking me so long ?",
don't be surprised.
And it's probably OT for this group. (I've found sci.crypt to be a very
useful newsgroup for matters pertaining to primes and factoring.) Sorry,
I can't really help with the on-topic aspects you raised.
Any thoughts, Kevin, on including a 'mpz_sieve_eratosthenes(max_prime)'
function with GMP ? I envisage a function that returns a mpz such that
if (and only if) bit x is set, then 2 * x is prime. Obviously such a mpz
does not encode the prime number '2' - but if you want to encode '2',
you double the size of the output mpz.
My sieve of eratosthenes is just a straight C impelementation, where the
output vector is, in fact, an array of ints. I don't know if there would
be significant gains in using the low level GMP functions to create an
mpz that encodes the same information.
Another issue might be that it is slower to iterate through an mpz bit
by bit than it is to iterate through an array of ints bit by bit.
Cheers,
Rob
--
Any emails containing attachments will be deleted from my ISP's mail
server before I even get to see them. If you wish to email me an
attachment, please provide advance warning so that I can make the
necessary arrangements.
More information about the gmp-discuss
mailing list