GPU Prime Sieve

[tapion64]

I've never programmed in CUDA before, so trying to read through how mfaktc does its prime sieving is a little difficult for me.

Essentially my goal is to duplicate the functionality of mfaktc, but instead of trial factoring it against a particular exponent, I want to do a multiplicative order of 2 calculation for that prime, and output the results. On the CPU, I've been using an implementation in Go that can be found here:

http://rosettacode.org/wiki/Multiplicative_order#Go

Unfortunately, I know nothing about how to do arbitrary precision calculations in CUDA, so I can't simply port that algorithm (and I'm sure it would need to be tweaked for performance anyways, judging on how the rest of the mfaktc code looks).

One question I do have, because I haven't fully followed how the prime sieving works: when we do the trial factoring, are we actually finding /every single prime/ between 2^n and 2^(n+1) and then taking the remainder? Or is there a sub-algorithm that's sieving out unnecessary primes for that exponent to check? I'm really hoping it's the former, since if it's the latter I wouldn't achieve what I want with this.

[tapion64]

So, the main component of this algorithm that slows it down on the CPU is factoring p-1. It should be possible to use msieve or something like it within CUDA to facilitate that, since we're working with 20-30 digit numbers.

So if the trial factoring prime sieve can hand us each and every prime between 2^64 and 2^65, and so forth, then we can run msieve to factor p-1, use the rest of Bach Shallit's algorithm to calculate the order, msieve on the order to determine primality, if the order is prime then we know that the prime divides Morder.

I just have to motivate myself to learn CUDA well enough to modify the mfaktc code. Would someone else like to do it? :p

[TheMawn]

Well I won't be able to help for very long, but it is known that any factors of 2^P-1 are of the form 2*k*P+1 (for prime P, I believe), so trial factoring only considers that form of number, i.e., k from 1 to whatever it takes to hit your upper bit limit.

A fun bit of trivia here is that for this exact reason, TF is actually faster for larger exponents because it takes less k's to hit a particular bit level.

[tapion64]

Originally Posted by TheMawn

Well I won't be able to help for very long, but it is known that any factors of 2^P-1 are of the form 2*k*P+1 (for prime P, I believe), so trial factoring only considers that form of number, i.e., k from 1 to whatever it takes to hit your upper bit limit.

A fun bit of trivia here is that for this exact reason, TF is actually faster for larger exponents because it takes less k's to hit a particular bit level.

Well, I thought that might be the case. Here I was hoping that it went through all of them. Although, now that I think about it, I might as well just make it all msieve. Come up with a reasonable Atkin sieve to pre-screen for primes, pass on to msieve to verify primality, use msieve on p-1, Bach Shallit algorithm to calculate the order, msieve on the order to verify primality. Sounds pretty reasonable to me. I still need to lrn2CUDA.

[TheMawn]

I'm struggling with the math a bit. I haven't done any number theory, only calculus, and I've repressed a lot of it. Let me see if I follow:

Multiplicative Order (MO) of n (mod m) is the smallest a such that n^a = 1 (mod m).

_{[SUB]My first problem is I don't really see the purpose of having a MO in the first place. I guess it's a nifty way to look at numbers, but that's a struggle for me.}[/SUB]

For our interests, we're taking n = 2 and a = P, so 2^P-1 = 0 (mod q) where q is a factor of M_P.

Going backward, the MO of 2 (mod q) is the smallest P such that 2^P = 1 (mod q)


Now, what you are proposing is a program that calculates the MO of 2 (mod q) for whatever q you punch in? Suppose that spits out a number C.

If the MO of 2 (mod q) is C, then you've found the smallest C such that 2^C - 1 = 0 (mod q).


Ooooooookayy...... If I am understanding this, the program you're proposing is going to find the smallest Mersenne Number for which a selected number is a factor?

...

I've got to imagine that the MO algorithm is either a lot slower than you think or that something has gone wrong here. (We need Doctor Silverman).

Do all numbers have a multiplicative order?

Wouldn't this immediately prove the existence of infinitely many Mersenne Primes? EDIT: Derp derp derp derp derp nevermind that. Holy crap je suis un Derp.

[TheMawn]

However, as for your initial question, I think in practical terms, the answer is all primes are candidates.

If you're tackling the problem as "Here's a factor. Let's find a Mersenne Number" instead of "Here's a Mersenne number. Let's find a factor", technically speaking all prime factors are candidate factors.

A number q can only be a factor for 2^P-1 if it is of the form q = 2kP+1, but q may be invalid for one Mersenne Number, but valid for the next.

For example, 89 is a factor of M₁₁, and 89 = 2*4*11+1. (the other factor is 23). 89 could also have been a factor for 22 (2*2*22+1=89) or 44 (2*1*44+1=89) but 22 and 44 are composite anyway.

I couldn't say for sure but my intuition tells me that, similar to 89, a number can only be a Mersenne Number factor once...

[tapion64]

Originally Posted by TheMawn

However, as for your initial question, I think in practical terms, the answer is all primes are candidates.

If you're tackling the problem as "Here's a factor. Let's find a Mersenne Number" instead of "Here's a Mersenne number. Let's find a factor", technically speaking all prime factors are candidate factors.

A number q can only be a factor for 2^P-1 if it is of the form q = 2kP+1, but q may be invalid for one Mersenne Number, but valid for the next.

For example, 89 is a factor of M₁₁, and 89 = 2*4*11+1. (the other factor is 23). 89 could also have been a factor for 22 (2*2*22+1=89) or 44 (2*1*44+1=89) but 22 and 44 are composite anyway.

I couldn't say for sure but my intuition tells me that, similar to 89, a number can only be a Mersenne Number factor once...

The point of the Multiplicative order is that once you know the order, you know the first (and only) Mersenne number it can be a factor of (if any). there's the possibility that its multiplicative order is composite (very highly likely, actually), in which case it will never divide a number of the form 2^p-1.

Once you know the multiplicative order k, you know that prime will only divide 2^k-1, 2^(2*k)-1, 2^(3*k)-1, and so forth. In other words, it will never divide another Mersenne number, if it divides one at all. But you're right, the feasibility of doing this depends on the both the prime sieve's speed and the speed of msieve to factor it. In any case, with the advent of msieve, I know we can find the multiplicative order in milliseconds. And from that we can immediately tell what Mersenne number it's a factor of, if any. The problem is that there are so many primes to eliminate. It would probably have to be a distributed project, not feasible to run on one GPU. But still, if it can eliminate even 2^64 through 2^65, that's something.

[LaurV]

Originally Posted by tapion64

The point of the Multiplicative order is that

is that you assume you have an efficient method to factor q-1, for any prime q. Which you don't have. See the other topic you opened. Behind 2^48 or so, this method is much MUCH slower than the method GIMPS (and mfaktX) uses. At 2^64 it gets like 60 THOUSAND times slower.

Maybe a supermod can join the two topics.