E
Beispiel: http://www.math.rwth-aachen.de/~Gabriele.Nebe/Vorl/SEMCA/Mersenne.pdf
https://de.wikipedia.org/wiki/Lucas-Lehmer-Test
https://de.wikiversity.org/wiki/Benutzer:Tanik/Lucas-Lehmer-Test_für_Mersennezahlen/Beispiel1
Mit dem (verbesserten) Lucas-Lehmer-Test:
#define NOMINMAX // due to conflict with max()
#include <windows.h> // warning: no C++ standard!
#include "MyBigUnsigned.h"
#include <iostream>
#include <cstdint>
using namespace std;
void textcolor(unsigned short color = 15)
{
SetConsoleTextAttribute(::GetStdHandle(STD_OUTPUT_HANDLE), color);
}
BigUnsigned calculateMersenneNumber(BigUnsigned p)
{
BigUnsigned x = 2;
for (BigUnsigned i = 1; i < p; ++i)
{
x <<= 1;
}
return (x - 1);
}
bool LucasLehmerTest(BigUnsigned p)
{
BigUnsigned m = calculateMersenneNumber(p), s = 4;
for (BigUnsigned i = 2; i < p; ++i)
s = (s*s - 2) % m;
return (s == 0);
}
int main()
{
uint64_t startCount = 0, lastCount = 0, delta = 0, frequency;
QueryPerformanceFrequency((LARGE_INTEGER*)&frequency);
cout << "cpu frequency: " << frequency << " Hz" << endl << endl;
BigUnsigned p;
cout << "Lucas-Lehmer-Test (Mersenne-Zahlen)" << endl; // https://de.wikipedia.org/wiki/Lucas-Lehmer-Test
QueryPerformanceCounter((LARGE_INTEGER*)&startCount);
for (BigUnsigned p = 0; p < 10000; ++p)
{
if (LucasLehmerTest(p))
{
QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
uint64_t delta = lastCount - startCount;
textcolor(15);
cout << "p: " << p;
textcolor(2);
cout << "\ttime: " << 1000 * delta / frequency << " ms" << endl;
}
else
{
//cout << " ist keine Mersenne-Primzahl" << endl;
}
}
}
cpu frequency: 3127021 Hz
Lucas-Lehmer-Test (Mersenne-Zahlen)
p: 3 time: 0 ms
p: 5 time: 0 ms
p: 7 time: 0 ms
p: 13 time: 1 ms
p: 17 time: 1 ms
p: 19 time: 2 ms
p: 31 time: 2 ms
p: 61 time: 6 ms
p: 89 time: 14 ms
p: 107 time: 23 ms
p: 127 time: 39 ms
p: 521 time: 6140 ms
p: 607 time: 10858 ms
p: 1279 time: 191493 ms
p: 2203 time: 1617345 ms
p: 2281 time: 1855752 ms
p: 3217 time: 7159214 ms
p: 4253 time: 21421841 ms
p: 4423 time: 25019798 ms
Die Ergebnisse sind korrekt: http://www.mersenne.org/primes/
Der Lucas-Lehmer-Test ist zeitlich besser als der iterierende Miller-Rabin-Test mit vorgelagertem Divisions-Precheck, aber für sehr hohe Zahlenbereiche, wie bei Mersenne-Primzahlen aktuell notwendig, leider immer noch lahm.
Probleme: BigUnsigned, Multiplikation (==> Multiplikation mit schneller Fourier-Transformation, Schöhage und Strassen, 1971), noch single threaded
Eine deutlich spürbare Beschleunigung erzielt man durch einen Divisions-Precheck 3 ... 293 (dort liegt in etwa das Optimum) ab p = 1000:
#define NOMINMAX // due to conflict with max()
#include <windows.h> // warning: no C++ standard!
#include "MyBigUnsigned.h"
#include <iostream>
#include <cstdint>
#include <vector>
#include <cassert>
using namespace std;
const uint32_t DIVISOR_PRECHECK_MAX = 293;
void textcolor(unsigned short color = 15)
{
SetConsoleTextAttribute(::GetStdHandle(STD_OUTPUT_HANDLE), color);
}
uint64_t convertBigToU64(BigUnsigned num)
{
BigUnsigned ZweiHochVierUndSechzig = stringToBigUnsigned("18446744073709551616");
assert(num < ZweiHochVierUndSechzig);
uint64_t jHi = num.getBlock(1);
uint64_t jLo = num.getBlock(0);
return (jHi << 32 | jLo);
}
void convertU64ToBig(uint64_t num, BigUnsigned& b)
{
uint32_t numLo = (uint32_t)(num & 0xFFFFFFFF);
uint32_t numHi = (uint32_t)((num >> 32) & 0xFFFFFFFF);
b.setBlock(0, numLo);
b.setBlock(1, numHi);
}
BigUnsigned calculateMersenneNumber(BigUnsigned p)
{
BigUnsigned x = 2;
for (BigUnsigned i = 1; i < p; ++i)
{
x <<= 1;
}
return (x - 1);
}
bool LucasLehmerTest(BigUnsigned p)
{
BigUnsigned m = calculateMersenneNumber(p), s = 4;
for (BigUnsigned i = 2; i < p; ++i)
s = (s*s - 2) % m;
return (s == 0);
}
// division test with first prime numbers
bool divisionPrecheck(vector<bool>& primes, BigUnsigned p)
{
BigUnsigned m = calculateMersenneNumber(p);
for (BigUnsigned divisor = 3; divisor <= DIVISOR_PRECHECK_MAX; ++divisor)
{
if (primes[convertBigToU64(divisor)])
{
if (m % divisor == 0)
{
return false;
}
}
}
return true;
}
int main()
{
uint64_t startCount = 0, lastCount = 0, delta = 0, frequency;
QueryPerformanceFrequency((LARGE_INTEGER*)&frequency);
cout << "cpu frequency: " << frequency << " Hz" << endl << endl;
////////////////////////// Generate primes lookup table /////////////////////////////////////////////////////////////
uint64_t const primeLimit = DIVISOR_PRECHECK_MAX;
cout << "We generate vector<bool>(primes) up to: " << primeLimit << endl;
vector<bool> primes(primeLimit + 1); // calculated primes (low values)
vector<bool>& refPrimes = primes;
primes[0] = primes[1] = false;
for (uint64_t i = 2; i < primeLimit + 1; ++i) { primes[i] = true; }
// Erastothenes sieving loop
uint64_t iMax = (uint64_t)sqrt((double)primeLimit) + 1;
for (uint64_t i = 2; i < iMax; ++i)
{
if (primes[i])
{
for (uint64_t j = 2 * i; j < primeLimit + 1; j += i)
{
primes[j] = false;
}
}
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
BigUnsigned p;
cout << "Lucas-Lehmer-Test (Mersenne-Zahlen)" << endl; // https://de.wikipedia.org/wiki/Lucas-Lehmer-Test
QueryPerformanceCounter((LARGE_INTEGER*)&startCount);
for (BigUnsigned p = 0; p < 10000; ++p)
{
if (p > 1000)
{
if (divisionPrecheck(primes, p) == false)
{
continue;
}
}
if (LucasLehmerTest(p))
{
QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
uint64_t delta = lastCount - startCount;
textcolor(15);
cout << "p: " << p;
textcolor(2);
cout << "\ttime: " << 1000 * delta / frequency << " ms" << endl;
}
else
{
//cout << " ist keine Mersenne-Primzahl" << endl;
}
}
}
cpu frequency: 3127021 Hz
We generate vector<bool>(primes) up to: 293
Lucas-Lehmer-Test (Mersenne-Zahlen)
p: 3 time: 0 ms
p: 5 time: 0 ms
p: 7 time: 0 ms
p: 13 time: 1 ms
p: 17 time: 1 ms
p: 19 time: 1 ms
p: 31 time: 2 ms
p: 61 time: 6 ms
p: 89 time: 14 ms
p: 107 time: 23 ms
p: 127 time: 39 ms
p: 521 time: 5895 ms
p: 607 time: 10473 ms
p: 1279 time: 90122 ms
p: 2203 time: 326176 ms
p: 2281 time: 369233 ms
p: 3217 time: 1261672 ms
p: 4253 time: 3680716 ms
p: 4423 time: 4260730 ms
Die Zeiten in ms sind die gesamte vergangene Zeit ab Start des Tests. Der optimale max. Primzahl-Teiler für das jeweilige p müsste hier erst noch untersucht werden.
