Optimizing Modulus (was Re: Would someone mind elaborating the explanation in the manual?)

[Sisyphus]

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.