Suche nach Mersenne-Primzahlen

[Erhard Henkes]

OK, dann nehmen wir das als bewiesen und checken nur prime p.

[SeppJ]

Erhard Henkes schrieb:

OK, dann nehmen wir das als bewiesen und checken nur prime p.

Wodurch dann automatisch die letzte Ziffer von 2^p-1 eine 1 oder 7 ist.

[Erhard Henkes]

Danke für die Kommentare und Hinweise.

Aktueller Code: http://pastebin.com/DjYAHdWX

Der 1-7-Test macht bisher nur Sinn bei p=2, also keine Beschleunigung. Den Divisions-Precheck habe ich auf 3 und 5 beschränkt. Für p werden nur Primes aus der Lookup-Tabelle verwendet. BigUnsigned könnte man noch gegen eine schnellere Unsigned-Klasse (z.B. gmp) tauschen.

cpu frequency: 3127021 Hz

We generate vector<bool>(primes) up to: 100000000
Lucas-Lehmer-Test (Mersenne-Zahlen)
 >>> 1-7-test makes sense for p = 2 <<<
p: 3    time: 5 ms
p: 5    time: 6 ms
p: 7    time: 6 ms
p: 13   time: 6 ms
p: 17   time: 7 ms
p: 19   time: 7 ms
p: 31   time: 8 ms
p: 61   time: 9 ms
p: 89   time: 12 ms
p: 107  time: 15 ms
p: 127  time: 18 ms
p: 521  time: 1766 ms
p: 607  time: 3116 ms
p: 1279 time: 46898 ms
p: 2203 time: 370715 ms
p: 2281 time: 429500 ms
work is done for 0 to 2399
p: 3217 time: 993202 ms
...

[SeppJ]

Eine Zahl, die auf 1 oder 7 endet, wird wohl kaum durch 5 teilbar sein. Und 2^p-1 kann auch nie durch 3 teilbar sein, denn wie im anderen Thread erklärt sind 2^n-1 für n gerade genau die Zahlen, die als Quersumme modulo 9 entweder 0, 3 oder 6 haben. Und alle ungeraden n haben 1, 7 oder 4, können also nicht durch 3 teilbar sein. Deine ganzen Vorschecks sind also allesamt sinnlos für p > 2.

[Erhard Henkes]

Ja, das macht keinen Sinn. Primzahl p als Input und Lucas-Lehmer-Test zum Sieben reichen.
Der Umweg über perfekte oder superperfekte Zahl dürfte durch die aufwändige Faktorisierung auch zu lahm sein. Eine Alternative bildet vlt. dieser "Elliptic-Test" http://www.math.rwth-aachen.de/~Gabriele.Nebe/Vorl/SEMCA/Mersenne.pdf

Das ist echt "rock bottom".

Momentan teste ich noch den Divisions-Precheck (7 ... 293) vor dem Lucas-Lehmer-Test. Zur Übersicht gebe ich den gefundenen Divisor für p aus:
http://pastebin.com/WFWeTVdY

Allerdings ergibt das nur wenige Funde. Dreht man DIVISOR_PRECHECK_MAX hoch, dann frisst das zu viel Zeit. Hier ein Beispiel: http://pastebin.com/WSRQhFGj
Könnte man vlt. in Abhängigkeit von p einstellen.

Aktueller Code: http://pastebin.com/cDWMMB7h

[Bengo]

Erhard Henkes schrieb:

Ja, wir bleiben nun hier.

Wenn p keine Primzahl ist, kann 2^p-1 keine Primzahl sein.

Ja, das kann man aussparen, wenn das bewiesen ist (habe es gelesen, man weiß aber oft nicht, was ist Beweis und was nur Vermutung. Umso besser, spart Zeit. Die p Primes kann man im vector<bool> primes alle vorher berechnen.

Der Satz ist ziemlich einfach zu beweisen. Ob die Umkehrung gilt, war bis 1536 unklar, und seit dem weiß man, dass sie nicht gilt.

Merke gerade, Umkehrung ist bei dem Satz nicht ganz eindeutig, ich meine:

Wenn p eine Primzahl ist, dann ist auch $2^p-1$ prim.

Wenn du auf Mersenne Primzahlen testen willst hilft das hier vielleicht weiter:

"Ist n eine ungerade Primzahl und q ein Primfaktor von Mn, so gilt q ≡ 1 (mod 2n) und q ≡ ±1 (mod 8). Beispiel: M11 = 2047 = 23·89 und 23 = 2·11+1; 89 = 4·2·11+1." Wikipedia Mersenne-Primzahl

[Erhard Henkes]

Mit meinem Programm schaffe ich es in wenigen Stunden gerade bis zu den Mersenne Primzahlen des Jahres 1963. Der Divisor-Precheck wird auch im Mersenne Forum als sinnvoll erachtet, habe aber bisher keine exakte Vorgehensweise gefunden, nur Forum-Gelaber.

@Bengo: Danke! Ich bin für jeden konstruktiven Hinweis dankbar bei diesem strohstrockenen mathematischen Thema.

Wenn p eine Primzahl ist, dann ist auch 2^p−1 prim.

Das gilt keinesfalls, sonst müsste man nicht den Lucas-Lehmer-Test http://page.math.tu-berlin.de/~kant/Mersenne/MersenneVortrag.pdf durchführen, der bei großen Zahlen die Schleife etwa p-mal durchläuft. Auf dem PC wird es so ab p=10000 ziemlich langatmig. Da fehlt ein genialer Precheck.

Als Alternative fand ich bisher nur den "Elliptic-Test" http://www.math.rwth-aachen.de/~Gabriele.Nebe/Vorl/SEMCA/Mersenne.pdf oder http://www.math.harvard.edu/~gross/preprints/Mersenne.pdf

[Erhard Henkes]

Durch das MT, das die Aufgabe für die 11 Threads in 11 Scheiben von je 2400 zu prüfende Mersenne-Zahlen aufteilt, kommt nun die Mersenneprimzahl für p=21701:

cpu frequency: 3127021 Hz
DIVISOR_PRECHECK_MAX: 293
number of threads: 11
Lucas Lehmer prime test for mersenne numbers
p: 3 time: 3 ms
p: 5 time: 4 ms
p: 7 time: 4 ms
p: 13 time: 5 ms
p: 17 time: 5 ms
p: 19 time: 6 ms

> divisionPrecheck finds 47 for p = 23 <<<
> divisionPrecheck finds 233 for p = 29 <<<
p: 31 time: 12 ms

> divisionPrecheck finds 223 for p = 37 <<<
p: 61 time: 27 ms

> divisionPrecheck finds 167 for p = 83 <<<
p: 89 time: 67 ms
p: 107 time: 113 ms
p: 127 time: 128 ms

> divisionPrecheck finds 263 for p = 131 <<<
p: 521 time: 2106 ms
p: 607 time: 3445 ms
p: 1279 time: 47946 ms
p: 2203 time: 369212 ms
p: 2281 time: 429714 ms
work is done for 0 to 2399
p: 3217 time: 997922 ms
p: 4253 time: 3862360 ms
p: 9689 time: 4496220 ms
p: 4423 time: 4509814 ms
work is done for 2400 to 4799
p: 9941 time: 15394338 ms
work is done for 4800 to 7199
p: 21701 time: 37340292 ms

Gefunden wurde diese Zahl 1978 auf einem "Supercomputer" CDC Cyber 174. Diese gigantische Zahl mit 6533 Dezimalstellen kann man heute - ohne große Mühe und besondere Kniffe - auf einem PC in knapp über 10 h finden.
http://www.ams.org/journals/mcom/1980-35-152/S0025-5718-1980-0583517-4/S0025-5718-1980-0583517-4.pdf (die damalige Rechenzeit betrug 7:40:20)
Interessant ist hier die Seite 2, auf der diese Tabelle von Wagstaff gezeigt wird, die für M(p) Faktoren der Formel 2*k*p + 1 aufzeigt. Die Tabelle zeigt p und k (jeweils bezogen auf den kleinsten Faktor).
Man könnte in diesem Sinne vlt. einen nutzbringenden Divisions-Precheck durchführen, in dem man k von 1 bis zu einem definierten k_max (z.B. 300) für Divisoren 2*k*p+1 durchlaufen lässt und auf Rest 0 testet.

Das wird einfach ausprobiert:

// division test with 2*k*p+1
BigUnsigned divisionPrecheck(BigUnsigned m, BigUnsigned p)
{
	for (BigUnsigned k = 1; k <= DIVISOR_PRECHECK_K_MAX; ++k)
	{		
		BigUnsigned Divisor = k * p;
		Divisor <<= 1; // 2 * k * p 
		Divisor += 1; // 2 * k * p +1 

		if (m % Divisor == 0)
		{
			return k;
		}		
	}

	return 0; // No divisor found
}

Da geht sofort die Post ab: http://pastebin.com/8PfgBpsP (DIVISOR_PRECHECK_K_MAX = 300)
Kleiner Ausschnitt:

>>> divisionPrecheck (2*k*p+1) finds k = 1 for p = 4391 <<<
>>> divisionPrecheck (2*k*p+1) finds k = 69 for p = 7411 <<<
>>> divisionPrecheck (2*k*p+1) finds k = 8 for p = 7417 <<<
p: 4423 time: 2874225 ms
>>> divisionPrecheck (2*k*p+1) finds k = 3 for p = 4441 <<<
>>> divisionPrecheck (2*k*p+1) finds k = 8 for p = 4447 <<<
>>> divisionPrecheck (2*k*p+1) finds k = 183 for p = 4457 <<<

Endlich ein brauchbarer Divisions-Precheck, der sinnlose Zeit im aufwändigen Lucas-Lehmer-Test spart (bei DIVISOR_PRECHECK_K_MAX = 300 ca. 25 - 35%).

Soll man das DIVISOR_PRECHECK_K_MAX konstant halten oder mit p hochziehen? Welcher Wert ist optimal? Zunächst ein direkter Vergleich mit dem doppelten Max-Wert für k:

Vergleich DIVISOR_PRECHECK_K_MAX: 300 vs. 600
http://pastebin.com/U4rNWGpm (DIVISOR_PRECHECK_K_MAX = 600)
p: 1279 time: 27160 ms vs. 25852 ms
p: 2203 time: 220527 ms vs. 213759 ms
p: 2281 time: 263063 ms vs. 258019 ms
p: 3217 time: 640619 ms vs. 600831 ms

Gewinner ist DIVISOR_PRECHECK_K_MAX: 600, was mich etwas wundert, denn es gibt viel mehr k unter 300 als darüber. Wir bleiben bei 600 und stoppen die Ausgaben der p u. k bei erfolgreichem Divisionscheck, damit das Programm auch schneller und unbeaufsichtigt laufen kann.

[Erhard Henkes]

Mittels vorgelagerter 2*k*p+1 Divisions- und anschließender Lucas Lehmer Prüfung lässt sich die Originalrechenzeit selbst mit der nicht optimalen Klasse BigUnsigned unterbieten:

cpu frequency: 3127021 Hz
DIVISOR_PRECHECK_K_MAX: 600
number of threads: 11
Lucas Lehmer prime test for mersenne numbers
p: 3    time: 4 ms
p: 5    time: 4 ms
p: 7    time: 5 ms
p: 13   time: 5 ms
p: 17   time: 6 ms
p: 19   time: 6 ms
p: 31   time: 7 ms
p: 61   time: 11 ms
p: 89   time: 16 ms
p: 107  time: 24 ms
p: 127  time: 32 ms
p: 521  time: 1032 ms
p: 607  time: 1648 ms
p: 1279 time: 23770 ms
p: 2203 time: 193230 ms
p: 2281 time: 233976 ms
work is done for 0 to 2399
p: 3217 time: 565830 ms
p: 4253 time: 2122783 ms
p: 4423 time: 2550675 ms
p: 9689 time: 3226074 ms
work is done for 2400 to 4799
p: 9941 time: 10907087 ms
work is done for 4800 to 7199
p: 21701        time: 20942929 ms
work is done for 7200 to 9599
p: 11213        time: 39180969 ms

Aktueller Code: http://pastebin.com/AizvmMFh

Der o.g. Divisions-Precheck reduziert die Rechenzeit für M(21701) von 37340 s auf 20943 s. Dies ist eine Ersparnis von 44%

-----------------------------------------------------

Nun steigen wir vorwitzig mit ein paar Veränderungen bei uint64_t/BigUnsigned, Bildschirmausgaben und im Threadbereich in die oberste Liga auf und testen damit die Laufzeit bei der höchsten z.Z. bekannten M(57885161):

#define NOMINMAX     // due to conflict with max()
#include <windows.h> // warning: no C++ standard!

#include "MyBigUnsigned.h" // http://henkessoft.de/Sonstiges/MyBigUnsigned.rar 
                           // based on https://mattmccutchen.net/bigint/
#include <iostream>
#include <cstdint>
#include <vector>
#include <cassert>
#include <thread>
#include <mutex>

using namespace std;

const uint32_t numberOfThreads = 11;
const uint32_t DIVISOR_PRECHECK_K_MAX = 600;
uint64_t const primeLimit = 200000000;
uint64_t frequency;

mutex mutex0;
mutex mutex1;

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(uint64_t p)
{
	BigUnsigned x = 2;

	for (uint64_t i = 1; i < p; ++i)
	{
		x <<= 1;
	}
	return (x - 1);
}

bool LucasLehmerTest(BigUnsigned m, uint64_t p, uint64_t startCount, uint64_t lastCount)
{
	BigUnsigned s = 4;

	for (uint64_t i = 2; i < p; ++i)
	{
		if(i%10 == 0)
		{ 
			QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
			uint64_t delta = lastCount - startCount;

			mutex1.lock();
			cout << "p = " << p << "  i = " << i;
			cout << "\ttime: " << 1000 * delta / frequency << " ms" << endl;
			mutex1.unlock();
		}	

		s = (s*s - 2) % m;
	}		

	return (s == 0);
}

// division test with 2*k*p+1
uint64_t divisionPrecheck(BigUnsigned m, uint64_t p)
{
	for (uint64_t k = 1; k <= DIVISOR_PRECHECK_K_MAX; ++k)
	{
		if (k % 100 == 0)
		{
			mutex1.lock();
			cout << "check: p = " << p << "  k = " << k << endl; 
			mutex1.unlock();
		}

		uint64_t Divisor = k * p;
		Divisor <<= 1; // 2 * k * p 
		Divisor += 1; // 2 * k * p +1 

		BigUnsigned BigDivisor(0);
		convertU64ToBig(Divisor, BigDivisor);

		if (m % BigDivisor == 0)
		{
			return k;
		}		
	}

	return 0; // No divisor found
}

void mersennePrimeCheck(vector<bool>& primes, uint64_t startCount, uint64_t lastCount, uint64_t startP, uint64_t limitP)
{
	for (uint64_t p = startP; p < limitP; ++p)
	{
		// primes lookup for p // In order for M(p) to be prime, p must itself be prime.
		if (!primes[p])
		{
			mutex0.lock();
			textcolor(3);
			cout << ">>> p = " << p << " no prime! <<<" << endl;
			mutex0.unlock();
			continue;
		}			

		// if p is prime, then we calculate the mersenne number
		BigUnsigned m = calculateMersenneNumber(p);

		// mersenne numbers with prime factors are no mersenne primes
		if (p > 19)
		{ 
			uint64_t divisor = divisionPrecheck(m,p);
			if (divisor !=0)
			{
				mutex0.lock();
				textcolor(3);
				cout << ">>> divisionPrecheck (2*k*p+1) finds k = " << divisor << " for p = " << p << " <<<" << endl;
				mutex0.unlock();
				continue;			
			}
		}		

		if (LucasLehmerTest(m, p, startCount, lastCount))
		{
			QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
			uint64_t delta = lastCount - startCount;

			mutex0.lock();
			textcolor(15);
			cout << "p: " << p; 
			textcolor(2);
			cout << "\ttime: " << 1000 * delta / frequency << " ms" << endl;
			mutex0.unlock();
		}
		else
		{
			cout << ">>> p = " << p << " LucasLehmerTest negative! <<<" << endl;
		}
	}

	mutex0.lock();
	textcolor(11);
	cout << "work is done for " << startP << " to " << limitP-1 << endl;
	mutex0.unlock();
}

int main()
{
	uint64_t startCount = 0, lastCount = 0, delta = 0;
	QueryPerformanceFrequency((LARGE_INTEGER*)&frequency);
	cout << "cpu frequency: " << frequency << " Hz" << endl;
	cout << "DIVISOR_PRECHECK_K_MAX: " << DIVISOR_PRECHECK_K_MAX << endl;
	cout << "number of threads: " << numberOfThreads << endl;

	////////////////////////// Generate primes lookup table /////////////////////////////////////////////////////////////

	vector<bool> primes(primeLimit + 1); // calculated primes (low values)

	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;
			}
		}
	}

	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	cout << "Lucas Lehmer prime test for mersenne numbers" << endl; // https://de.wikipedia.org/wiki/Lucas-Lehmer-Test

	QueryPerformanceCounter((LARGE_INTEGER*)&startCount);

	vector<thread> t;
	uint64_t deltaP = 2;
	uint64_t startP = 110502;  //57885160;
	uint64_t limitP = startP + deltaP;

	for (int x = 0; x < numberOfThreads; ++x)
	{
		t.emplace_back(mersennePrimeCheck, primes, startCount, lastCount, startP + x * deltaP, limitP + x * deltaP);
	}

	for (int x = 0; x < numberOfThreads; ++x)
	{
		t[x].join(); // main has to wait for the results of the threads. 
	}
}

Durch unser Prime Lookup für p sieht man rasch diese Ausgabe:

cpu frequency: 3127021 Hz
DIVISOR_PRECHECK_K_MAX: 600
number of threads: 11
Lucas Lehmer prime test for mersenne numbers
>>> p = 57885160 no prime! <<<
>>> p = 57885162 no prime! <<<
>>> p = 57885163 no prime! <<<
work is done for 57885162 to 57885163
>>> p = 57885164 no prime! <<<
>>> p = 57885165 no prime! <<<
work is done for 57885164 to 57885165
>>> p = 57885166 no prime! <<<
>>> p = 57885168 no prime! <<<
>>> p = 57885169 no prime! <<<
work is done for 57885168 to 57885169
>>> p = 57885170 no prime! <<<
>>> p = 57885171 no prime! <<<
work is done for 57885170 to 57885171
>>> p = 57885172 no prime! <<<
>>> p = 57885174 no prime! <<<
>>> p = 57885175 no prime! <<<
work is done for 57885174 to 57885175
>>> p = 57885176 no prime! <<<
>>> p = 57885177 no prime! <<<
work is done for 57885176 to 57885177
>>> p = 57885178 no prime! <<<
>>> p = 57885179 no prime! <<<
work is done for 57885178 to 57885179
>>> p = 57885180 no prime! <<<
>>> p = 57885181 no prime! <<<
work is done for 57885180 to 57885181

P = 57885161, 57885167 und 57885173 sind prim und daher noch im Rennen im Precheck und ggf. im Lucas-Lehmer-Test, der nun in der inneren p-Schleife fast 58 Millionen (!) mal durchlaufen wird. Bin gespannt, ob es damit überhaupt noch geht in akzeptabler Zeit. Ansonsten: Andere BigUnsigned-Lib und beschleunigte Multiplikation. Wir sind am Punkt.

Kann man vergessen! daher zurück zu "kleineren" Zahlen, nämlich M(110503), damit wir etwas Zuschauer-Feeling erleben am Bildschirm und ein Gespür für den zeitlichen Ablauf erhalten:

cpu frequency: 3127021 Hz
DIVISOR_PRECHECK_K_MAX: 600
number of threads: 11
Lucas Lehmer prime test for mersenne numbers
>>> p = 110502 no prime! <<<
>>> p = 110504 no prime! <<<
>>> p = 110505 no prime! <<<
work is done for 110504 to 110505
>>> p = 110506 no prime! <<<
>>> p = 110507 no prime! <<<
work is done for 110506 to 110507
>>> p = 110508 no prime! <<<
>>> p = 110509 no prime! <<<
work is done for 110508 to 110509
>>> p = 110510 no prime! <<<
>>> p = 110511 no prime! <<<
work is done for 110510 to 110511
>>> p = 110512 no prime! <<<
>>> p = 110513 no prime! <<<
work is done for 110512 to 110513
>>> p = 110514 no prime! <<<
>>> p = 110515 no prime! <<<
work is done for 110514 to 110515
>>> p = 110516 no prime! <<<
>>> p = 110517 no prime! <<<
work is done for 110516 to 110517
>>> p = 110518 no prime! <<<
>>> p = 110519 no prime! <<<
work is done for 110518 to 110519
>>> p = 110520 no prime! <<<
>>> p = 110521 no prime! <<<
work is done for 110520 to 110521
>>> p = 110522 no prime! <<<
>>> p = 110523 no prime! <<<
work is done for 110522 to 110523
check: p = 110503  k = 100
...
check: p = 110503  k = 600
p = 110503  i = 10      time: 50826 ms
p = 110503  i = 20      time: 57404 ms
p = 110503  i = 30      time: 86831 ms
p = 110503  i = 40      time: 116222 ms
p = 110503  i = 50      time: 145534 ms
...
p = 110503  i = 140     time: 407994 ms
p = 110503  i = 150     time: 437177 ms
...
p = 110503  i = 230     time: 671057 ms
p = 110503  i = 240     time: 700347 ms
...
p = 110503  i = 620     time: 1806152 ms
p = 110503  i = 630     time: 1834953 ms

Man sieht, dass die Zeit innerhalb des Lucas-Lehmer-Tests zunehmend ist, sich dann bei ca. 29 s pro 10 Schritte einpegelt. Bei einer Zahl wie p = 110503 (Walter Colquitt & Luke Welsh, NEC SX-2, Jan 1988) beissen wir uns bereits zeitlich die Zähne aus. Das wird ca. 4 Tage dauern, bis man i auf p hochgezählt hat. Das ist also aktuell der Grenzbereich für die Bestimmung, wenn wir wissen(!), wo die Zahl sich befindet, sie also schon kennen.

Noch ein ernüchternder zeitlicher Check (2*k*p+1 - Divisions-Precheck ist aus Zeitersparnis ausgeschaltet) für die Bestimmung von M(1257787), Slowinski & Gage, 1996:

p = 1257787 i = 2 s(len) = 1    time: 76741 ms
p = 1257787 i = 3 s(len) = 1    time: 76743 ms
p = 1257787 i = 4 s(len) = 1    time: 76744 ms
p = 1257787 i = 5 s(len) = 1    time: 76745 ms
p = 1257787 i = 6 s(len) = 2    time: 76746 ms
p = 1257787 i = 7 s(len) = 4    time: 76746 ms
p = 1257787 i = 8 s(len) = 8    time: 76747 ms
p = 1257787 i = 9 s(len) = 16   time: 76748 ms
p = 1257787 i = 10 s(len) = 31  time: 76749 ms
p = 1257787 i = 11 s(len) = 61  time: 76750 ms
p = 1257787 i = 12 s(len) = 122 time: 76751 ms
p = 1257787 i = 13 s(len) = 244 time: 76754 ms
p = 1257787 i = 14 s(len) = 487 time: 76762 ms
p = 1257787 i = 15 s(len) = 973 time: 76789 ms
p = 1257787 i = 16 s(len) = 1946        time: 76896 ms
p = 1257787 i = 17 s(len) = 3892        time: 77318 ms
p = 1257787 i = 18 s(len) = 7783        time: 78990 ms
p = 1257787 i = 19 s(len) = 15565       time: 85604 ms
p = 1257787 i = 20 s(len) = 31130       time: 112214 ms
p = 1257787 i = 21 s(len) = 39306       time: 345827 ms
p = 1257787 i = 22 s(len) = 39306       time: 737912 ms
p = 1257787 i = 23 s(len) = 39306       time: 1131007 ms
p = 1257787 i = 24 s(len) = 39306       time: 1525815 ms
p = 1257787 i = 25 s(len) = 39306       time: 1921041 ms

Das Problem besteht darin, dass s*s eine gewaltige Größe annimmt. An einem gewissen i wird diese Größe - bestimmt durch BigUnsigned::getLength() - konstant. Die Zeitdifferenz für einen Schritt im Lucas-Lehmer-Test beträgt dann 393 s. Bei p-2 Schritten wären dies 1257787 * 393 s = ca. 15,7 Jahre!

Klares KO in dieser Liga für meinen PC mit meiner geschätzten BigUnsigned-Klasse und klassischer Multiplikation.

https://primes.utm.edu/notes/1257787.html (September 1996)

The proof of this 378,632 digit number's primality (using the traditional Lucas-Lehmer test) took about 6 hours on one CPU of a CRAY T94 super computer.

Beeindruckend.

Mal sehen, ob da nicht noch etwas geht.

[Erhard Henkes]

Momentan bringt das mich noch nicht weiter, aber Bengo hat in einer Literaturstelle noch folgenden Precheck auf der Ebene p gefunden:

bool precheck_mod4_is_3_and_2p_plus_1_isPrime (vector<bool>& primes, uint64_t p)
{
    return ( (p % 4 == 3) && (primes[2 * p + 1]));
}

//...
        if (p > 3)
        {
            if (precheck_mod4_is_3_and_2p_plus_1_isPrime (primes, p))
            {              
                cout << ">>> p = " << p << " p_mod4_is_3_and_2p_plus_1_isPrime! <<<" << endl;              
                continue;
            }
        }

Im unteren p-Bereich p<10000 sieht man keinen großen Vorteil (so wie beim Divisor-Check), aber bei großen p sollte sich alles positiv auswirken, was M(p) nicht anfassen muss und trotzdem aussortieren kann. Da gibt es bisher zwei Tests auf Ebene p: 1) isPrime(p) 2) precheck_mod4_is_3_and_2p_plus_1_isPrime

Aktueller Code (incl. dieses Prechecks): http://pastebin.com/5BjQgiUN

Vielen Dank an Bengo für die praktische Unterstützung!

Jetzt fehlt noch die geniale BigUnsigned-Lib und die blitzschnelle Multiplikation (in C++ code).

[Erhard Henkes]

Die bisher genutzte Klasse BigUnsigned ist zwar sehr schön einsehbar und nachvollziehbar, aber sie ist einfach zu langsam für Mersennezahlen.

Daher verwende ich ab sofort das auf GMP basierende MPIR 2.7.0 (mpir.dll wurde von mir mit MS VS 2015 für x64 geschwindigkeits-optimiert erstellt):

Aktueller Code:

#define NOMINMAX     // due to conflict with max()
#include <windows.h> // warning: no C++ standard!

#include <mpirxx.h>
#include <iostream>
#include <cstdint>
#include <vector>
#include <cassert>
#include <thread>
#include <mutex>

///#define DETAILS
///#define OUTPUT_OF_PRECHECKS
///#define NEGATIVE_LL

using namespace std;

const uint32_t numberOfThreads = 11;
const uint32_t DIVISOR_PRECHECK_K_MAX = 600;
uint64_t const primeLimit = 200000000;
uint64_t frequency;

mutex mutex0;
mutex mutex1;
mutex mutex2;
mutex mutex3;

void textcolor(unsigned short color = 15)
{
	SetConsoleTextAttribute(::GetStdHandle(STD_OUTPUT_HANDLE), color);
}

mpz_class calculateMersenneNumber(uint64_t p)
{
	mpz_class x = 2;

	for (uint64_t i = 1; i < p; ++i)
	{
		x <<= 1;
	}
	return (x - 1);
}

bool LucasLehmerTest(mpz_class m, uint64_t p, uint64_t startCount, uint64_t lastCount)
{
	mpz_class s = 4;

	for (uint64_t i = 2; i < p; ++i)
	{
		s = (s*s - 2) % m;

#ifdef DETAILS 
		QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
		uint64_t delta = lastCount - startCount;

		mutex1.lock();
		cout << "p = " << p << " i = " << i << " s(len) = " << s.getLength();
		cout << "\ttime: " << 1000 * delta / frequency << " ms" << endl;
		mutex1.unlock();
#endif
	}

	return (s == 0);
}

// division test with 2*k*p+1
uint64_t divisionPrecheck(mpz_class m, uint64_t p, uint64_t startCount, uint64_t lastCount)
{
	for (uint64_t k = 1; k <= DIVISOR_PRECHECK_K_MAX; ++k)
	{
		mpz_class Divisor = k * p;
		Divisor <<= 1; // 2 * k * p 
		Divisor += 1; // 2 * k * p +1 

#ifdef DETAILS 
		QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
		uint64_t delta = lastCount - startCount;

		mutex2.lock();
		cout << "check: p = " << p << "  k = " << k;
		cout << "\ttime: " << 1000 * delta / frequency << " ms" << endl;
		mutex2.unlock();
#endif

		if (m % Divisor == 0)
		{
			return k;
		}
	}

	return 0; // No divisor found
}

bool precheck_mod4_is_3_and_2p_plus_1_isPrime (vector<bool>& primes, uint64_t p)
{
	return ((p % 4 == 3) && (primes[2 * p + 1]));
}

void mersennePrimeCheck(vector<bool>& primes, uint64_t startCount, uint64_t lastCount, uint64_t startP, uint64_t limitP)
{
	for (uint64_t p = startP; p < limitP; ++p)
	{
		// primes lookup for p // In order for M(p) to be prime, p must itself be prime.
		if (!primes[p])
		{
#ifdef DETAILS 
			mutex0.lock();
			textcolor(3);
			cout << ">>> p = " << p << " no prime! <<<" << endl;
			mutex0.unlock();
#endif
			continue;
		}

		// precheck: P = 3 (mod 4) && isPrime (2*p+1) ?
		if (precheck_mod4_is_3_and_2p_plus_1_isPrime (primes, p))
		{
#ifdef OUTPUT_OF_PRECHECKS
			mutex3.lock();
			textcolor(5);
			cout << ">>> p = " << p << " p_mod4_is_3_and_2p_plus_1_isPrime! <<<" << endl;
			mutex3.unlock();
#endif
			continue;
		}

		// if p is prime, then we calculate the mersenne number
		mpz_class m = calculateMersenneNumber(p);

		// mersenne numbers with (2*k*p+1) factors are no mersenne primes
		if (p > 19)
		{
			uint64_t divisor = divisionPrecheck(m, p, startCount, lastCount);
			if (divisor != 0)
			{
#ifdef OUTPUT_OF_PRECHECKS
				mutex0.lock();
				textcolor(3);
				cout << ">>> divisionPrecheck (2*k*p+1) finds k = " << divisor << " for p = " << p << " <<<" << endl;
				mutex0.unlock();
#endif
				continue;
			}
		}

		if (LucasLehmerTest(m, p, startCount, lastCount))
		{
			QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
			uint64_t delta = lastCount - startCount;

			mutex0.lock();
			textcolor(15);
			cout << "p: " << p;
			textcolor(2);
			cout << "\ttime: " << 1000 * delta / frequency << " ms" << endl;
			mutex0.unlock();
		}
#ifdef NEGATIVE_LL
		else
		{
			cout << ">>> p = " << p << " LucasLehmerTest negative! <<<" << endl;
		}
#endif
	}

	mutex0.lock();
	textcolor(11);
	cout << "work is done for " << startP << " to " << limitP - 1 << endl;
	mutex0.unlock();
}

int main()
{
	uint64_t startCount = 0, lastCount = 0, delta = 0;
	QueryPerformanceFrequency((LARGE_INTEGER*)&frequency);
	cout << "cpu frequency: " << frequency << " Hz" << endl;
	cout << "DIVISOR_PRECHECK_K_MAX: " << DIVISOR_PRECHECK_K_MAX << endl;
	cout << "number of threads: " << numberOfThreads << endl;

	////////////////////////// Generate primes lookup table /////////////////////////////////////////////////////////////

	vector<bool> primes(primeLimit + 1); // calculated primes (low values)

	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;
			}
		}
	}

	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	cout << "Lucas Lehmer prime test for mersenne numbers" << endl; // https://de.wikipedia.org/wiki/Lucas-Lehmer-Test

	QueryPerformanceCounter((LARGE_INTEGER*)&startCount);

	vector<thread> t;
	uint64_t deltaP = 2400;
	uint64_t startP = 0;  // 57885160; 
	uint64_t limitP = startP + deltaP;

	for (int x = 0; x < numberOfThreads; ++x)
	{
		t.emplace_back(mersennePrimeCheck, primes, startCount, lastCount, startP + x * deltaP, limitP + x * deltaP);
	}

	for (int x = 0; x < numberOfThreads; ++x)
	{
		t[x].join(); // main has to wait for the results of the threads. 
	}
}

Das Ergebnis ist gegenüber der nun umgehend in den Ruhestand verabschiedeten Klasse BigUnsigned mehr als überzeugend:

cpu frequency: 3127021 Hz
DIVISOR_PRECHECK_K_MAX: 600
number of threads: 11
Lucas Lehmer prime test for mersenne numbers
p: 5    time: 7 ms
p: 7    time: 8 ms
p: 13   time: 8 ms
p: 17   time: 8 ms
p: 19   time: 9 ms
p: 31   time: 9 ms
p: 61   time: 10 ms
p: 89   time: 11 ms
p: 107  time: 12 ms
p: 127  time: 13 ms
p: 521  time: 29 ms
p: 607  time: 35 ms
p: 1279 time: 203 ms
p: 2203 time: 1395 ms
p: 2281 time: 1687 ms
p: 3217 time: 3600 ms
p: 4253 time: 12141 ms
p: 4423 time: 14403 ms
p: 9689 time: 15894 ms
p: 9941 time: 53003 ms
p: 21701        time: 79777 ms
p: 11213        time: 185587 ms
p: 19937        time: 368181 ms
p: 23209        time: 1153250 ms

Man kann MPIR für geschwindigkeitskritische Aufgaben wirklich empfehlen. BigUnsigned macht nur Sinn, wenn man selbst in der Klasse herum werkeln will. das kann man bei GMP vergessen.

Es fehlt noch die rasante verbesserte FT-Multiplikation (für LL Test), oder ist diese in GMP schon eingebaut? Könnte ich mir bei diesen rasanten Zeiten vorstellen. EDIT: Ist bereits eingebaut bei GMP.

Nun können wir uns mit dem PC wieder an die wirklich großen Zahlen heran wagen.

Für den Bereich ab p = 100000 ist das Programm ideal. Im Bereich p = 750000 wird es allerdings ebenfalls zäh (ca. 21 h pro Zahl pro Thread mit dem Lucas-Lehmer-Test). Es wird Zeit, dass der Lucas-Lehmer-Test durch etwas genialeres ersetzt wird, aber GIMPS verwendet ihn ebenfalls.