Talk:Lucas-Lehmer test: Difference between revisions
Content added Content deleted
(→Math in introduction: deleted (done, no need to discuss)) |
(→Speeding things up: new section) |
||
Line 327: | Line 327: | ||
The Java version is still limited by types. Integer.parseInt(args[0]) limits p to 2147483647. Also the fact that isMersennePrime takes an int limits it there too. For full arbitrary precision every int needs to be a BigInteger or BigDecimal and a square root method will need to be created for them. The limitation is OK I think (I don't think we'll be getting up to 2<sup>2147483647</sup> - 1 anytime soon), but the claim "any arbitrary prime" is false because of the use of ints. --[[User:Mwn3d|Mwn3d]] 07:45, 21 February 2008 (MST) |
The Java version is still limited by types. Integer.parseInt(args[0]) limits p to 2147483647. Also the fact that isMersennePrime takes an int limits it there too. For full arbitrary precision every int needs to be a BigInteger or BigDecimal and a square root method will need to be created for them. The limitation is OK I think (I don't think we'll be getting up to 2<sup>2147483647</sup> - 1 anytime soon), but the claim "any arbitrary prime" is false because of the use of ints. --[[User:Mwn3d|Mwn3d]] 07:45, 21 February 2008 (MST) |
||
== Speeding things up == |
|||
The main loop in Lucas-Lehmer is doing (n*n) mod M where M=2^p-1, and p > 1. '''But we can do it without division'''. |
|||
We compute the remainder of a division by M. Now, intuitively, dividing by 2^p-1 is almost like dividing by 2^p, except the latter is much faster since it's a shift. |
|||
Let's compute how much the two divisions differ. |
|||
We will call S = n*n. Notice that since the remainder mod M is computed again and again, the value of n must be < M at the beginning of a loop, that is at most 2^p-2, |
|||
thus S = n*n <= 2^(2*p) - 4*2^p + 4 = 2^p * (2^p - 2) + 4 - 2*2^p |
|||
When dividing S by M, you get quotient q1 and remainder r1 with S = q1*M + r1 and 0 <= r1 < M |
|||
When dividing S by M+1, you get likewise S = q2*(M+1) + r2 and 0 <= r2 <= M |
|||
In the latter, we divide by a larger number, so the quotient must be less, or maybe equal, that is, q2 <= q1. |
|||
Subtract the two equalities, giving |
|||
0 = (q2 - q1)*M + q2 + r2 - r1 |
|||
(q1 - q2)*M = q2 + r2 - r1 |
|||
Since S = 2^p * (2^p - 2) + 4 - 2*2^p |
|||
<= 2^p * (2^p - 2), |
|||
then the quotient q2 is less than 2^p - 2 (remember, when computing q2, we divide by M+1 = 2^p). |
|||
Now, 0 <= q2 <= 2^p - 2 |
|||
0 <= r2 <= 2^p - 1 |
|||
0 <= r1 <= 2^p - 2 |
|||
Thus the right hand side is >= 0, and <= 2*2^p - 3. |
|||
The left hand side is a multiple of M = 2^p - 1. |
|||
Therefore, this multiple must be 0*M or 1*M, certainly not 2*M = 2*2^p - 2 which would be > 2*2^p - 3, and not any other higher multiple would do. |
|||
So we have proved that q1 - q2 = 0 or 1. |
|||
This means that division by 2^p is almost equivalent (regarding the quotient) to dividing by 2^p-1: it's the same quotient, or maybe too short by 1. |
|||
Now, the remainder S mod M. |
|||
We start with a quotient q = S div 2^p, or simply q = S >> p (right shift). |
|||
The remainder is S - q*M = S - q*(2^p - 1) = S - q*2^p + q, and the multiplication by 2^p is a left shift. |
|||
And this remainder may be a bit too large, if our quotient is a bit too small (by one): in this case we must subtract M. |
|||
So, in pseudo-code, we are done if we do: |
|||
S = n*n |
|||
q = S >> p |
|||
r = S - (q << p) + q |
|||
if r >= M then r = r - M |
|||
We can go a bit further: taking S >> p then q << p is simply keeping the higher bits of S. |
|||
But then we subtract these higher bits from S, so we only keep the lower bits, that is we do (S & mask), and this mask is simply M ! (remember, M = 2^p - 1, a bit mask of p bits equal to "1") |
|||
The pseudo-code can thus be written |
|||
S = n*n |
|||
r = (S & M) + (S >> p) |
|||
if r >= M then r = r - M |
|||
And we have computed a remainder mod M without any division, only a few addition/subtraction/shift/bitwise-and, which will be much faster (each has a linear time complexity). |
|||
How much faster ? For exponents between 1 and 2000, in Python, the job is done 2.87 times as fast. For exponents between 1 and 5000, it's 3.42 times as fast. And it gets better and better, since the comlexity is lower. |
|||
[[User:Arbautjc|Arbautjc]] ([[User talk:Arbautjc|talk]]) 22:04, 15 November 2013 (UTC) |