Talk:Lucas-Lehmer test: Difference between revisions
Content added Content deleted
m (→Speeding things up: to keep indentation as is) |
|||
Line 332: | Line 332: | ||
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'''. |
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'''. |
||
<lang>We compute the remainder of a division by M. Now, intuitively, dividing by 2^p-1 is almost |
<lang>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. |
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 |
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 |
||
thus S = n*n <= 2^(2*p) - 4*2^p + 4 = 2^p * (2^p - 2) + 4 - 2*2^p |
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, you get quotient q1 and remainder r1 with S = q1*M + r1 and 0 <= r1 < M |
||
Line 356: | Line 357: | ||
The left hand side is a multiple of M = 2^p - 1. |
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 |
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. |
So we have proved that q1 - q2 = 0 or 1. |
||
This means that division by 2^p is almost equivalent (regarding the quotient) |
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. |
Now, the remainder S mod M. |
||
Line 374: | Line 377: | ||
We can go a bit further: taking S >> p then q << p is simply keeping the higher bits of S. |
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, |
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 |
The pseudo-code can thus be written |
||
Line 382: | Line 386: | ||
if r >= M then r = r - M |
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, |
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 |
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. |
|||
</lang> |
</lang> |
||
[[User:Arbautjc|Arbautjc]] ([[User talk:Arbautjc|talk]]) 22:04, 15 November 2013 (UTC) |
[[User:Arbautjc|Arbautjc]] ([[User talk:Arbautjc|talk]]) 22:04, 15 November 2013 (UTC) |