Goldbach Vermutung

[Erhard Henkes]

Jetzt beschleunigen wir das aber noch an anderer Stelle, wie bereits von euch vorgeschlagen:

bool IsPrime(vector<bool>& primes, BigUnsigned number, BigUnsigned primeLimit)
{
	if (number <= primeLimit)
		return primes[number.toUnsignedLong()]; // lookup from bitset (in the RAM)

	if ((number % 3 == 0) || (number % 5 == 0) || (number % 7 == 0) || (number % 11 == 0) || (number % 13 == 0) ||
		(number % 17 == 0) || (number % 19 == 0) || (number % 23 == 0) || (number % 29 == 0) || (number % 31 == 0) ||
		(number % 37 == 0) || (number % 41 == 0) || (number % 43 == 0))
	{
		return false;
	}

	return IsPrimeMillerRabinOptimized(number);
}

Vergleich bei 10^24:

Ohne die Abfrage number % ... == 0 mit 3 ... 43:

`Now Goldbach first prime number pair will be detected

1000000000000000000000000 257 time: 1887 ms

1000000000000000000000002 349 time: 2449 ms

1000000000000000000000004 307 time: 2182 ms

1000000000000000000000006 263 time: 1908 ms

1000000000000000000000008 311 time: 2172 ms

1000000000000000000000010 3 time: 33 ms

1000000000000000000000012 5 time: 61 ms

1000000000000000000000014 7 time: 92 ms

1000000000000000000000016 487 time: 3157 ms

1000000000000000000000018 11 time: 122 ms

1000000000000000000000020 13 time: 152 ms`

Mit dieser Abfrage vor dem Milli-Miller-Rabin-Test:

`Now Goldbach first prime number pair will be detected

1000000000000000000000000 257 time: 520 ms

1000000000000000000000002 349 time: 931 ms

1000000000000000000000004 307 time: 561 ms

1000000000000000000000006 263 time: 478 ms

1000000000000000000000008 311 time: 798 ms

1000000000000000000000010 3 time: 33 ms

1000000000000000000000012 5 time: 64 ms

1000000000000000000000014 7 time: 65 ms

1000000000000000000000016 487 time: 636 ms

1000000000000000000000018 11 time: 35 ms

1000000000000000000000020 13 time: 64 ms

`

Eine spürbare Beschleunigung.

Nun fehlt noch die Lügnerliste als unordered_set<BigUnsigned>.

[volkard]

Wäre nett, wenn Du den Bezeichner IsPrime aufheben würdest für nicht-lügende Implementierungen.

[Bengo]

Probier mal aus, bei wie vielen Primzahlen, die du vorher prüfst das optimum ist. Hab das Gefühl dass 43 schon etwas zu viel ist.

Und noch mal ein kleiner neuer Input: AKS-Primzahltest, er ist immer richtig, aber etwas langsamer als MillerRabin.

[TGGC]

Ich wuerde sogar mal testen, einige der Divisionen vor dem Lookup zu machen. Memory Zugriff ist je nach System und Zustand der Caches extrem teuer. Ausserdem koennte man so die den vector<bool> um ein Vielfalches kuerzer machen, weil ich dort drin dann noch noch Werte speichern muss, die nicht durch die ersten n Primzahlen teilbar sind. Oder probier mal das Sieb des Eratosthenes fuer die ersten n Primzahlen vorzuberechnen, das koennte klein genug fuer den Cache sein.

[Erhard Henkes]

Danke für die Ideen und Hinweise. Wie ihr seht, befindet sich eine Primzahl bei der Überprüfung der Goldbach-Vermutung im Bereich < 1000. Dieser Zugriff ist bezüglich Geschwindigkeit kein bottleneck. Interessant wird es bei isPrime(Geradzahl minus erste Primzahl). Im obigen Fall benötigt man hierfür je zu prüfender komplementärer Primzahl knapp über 30 ms.

3,5,7, ...

1000000000000000000000010 3 t: 33 ms (test auf n-3)
1000000000000000000000012 5 t: 61 ms (test auf n-3, n-5)
1000000000000000000000014 7 t: 92 ms (test auf n-3, n-5, n-7)
1000000000000000000000018 11 t: 122 ms (test auf n-3, n-5, n-7, n-11)
1000000000000000000000020 13 t: 152 ms (test auf n-3, n-5, n-7, n-11, n-13)

mit Divisions-Precheck:
1000000000000000000000014 7 time: 65 ms

*Man erkennt sofort, dass hier einer von drei Miller-Rabin-Tests eingespart werden konnte.
333333333333333333333337 * 3 = 1000000000000000000000011
1000000000000000000000011 + 3 = 1000000000000000000000014
*

Die Milli-Rabin-Lügnerliste ist für den Bereich, den ich nun mittels BigInteger untersuche, keine Hilfe. Im oberen Bereich fehlen die Lügner-Werte und unten brauchen wir es nicht, weil wir (via Sieb des E.) zuverlässige Primzahlen im Speicher haben.

Ich werde eure Vorschläge bei der nächsten Version einfließen lassen. Problematisch war noch der Austausch zwischen uint64_t und BigInteger. Da fehlen Funktionen in dieser Klasse. Habe das versuchsweise ergänzt. das ist wichtig für die Lügnerliste, die bei Rückgabe von true bei Miller-Rabin in gewissen Zahlenbereichen noch in Aktion treten sollte.

uint64_t convertBigToU64 (BigUnsigned num)
{
    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);
}

int main()
{
	BigUnsigned startNumber = stringToBigUnsigned("1000000000000000000");
	BigUnsigned endNumber   = stringToBigUnsigned("1000000000000001000");

    uint64_t j = convertBigToU64(endNumber);
    cout << j << endl;

    BigUnsigned big;
    BigUnsigned& refBig = big;
    convertU64ToBig (j, refBig);
    cout << big << endl;

	wait();
}

Eine eigene auf die hier notwendigen Operationen abgespeckte MyBigInteger-Klasse, die intern mit Blöcken der Größe uint64_t (anstelle uint32_t) arbeitet, könnte hier beschleunigend wirken.

[Bengo]

Kein Plan, ob es schon gute Lügnerlisten gibt. (Wenn du suchst, besser nach Pseudoprimzahlen suchen ;))
Aber du kannst dir welche erstellen auch für den sehr hohen Bereich. Wenn du MillerRabin, mit 2,3,5,7 machst, sind die Lügner auch bis in den riesigen Bereich überschaubar. Die Liste kannst aber wahrscheinlich nur mit dem AKS-Test erstellen. Division ist wahrscheinlich zu langsam.

[Erhard Henkes]

Erster Test (Notebook):

`vor Lookup: nix

vor M.R.: 3 ... 43

1000000000000000000000010 3 t: 47 ms

1000000000000000000000012 5 t: 109 ms

1000000000000000000000014 7 t: 109 ms

1000000000000000000000016 487 t: 1014 ms

1000000000000000000000018 11 t: 47 ms

1000000000000000000000020 13 t: 109 ms

---

vor Lookup: nix

vor M.R.: 3 ... 17

1000000000000000000000010 3 t: 47 ms

1000000000000000000000012 5 t: 109 ms

1000000000000000000000014 7 t: 109 ms

1000000000000000000000016 487 t: 1389 ms <===

1000000000000000000000018 11 t: 47 ms

1000000000000000000000020 13 t: 109 ms

---

vor Lookup: 3 ... 17

vor M.R.: 3 ... 43

1000000000000000000000010 313 t: 998 ms

1000000000000000000000012 269 t: 624 ms

1000000000000000000000014 271 t: 1186 ms

1000000000000000000000016 487 t: 920 ms

1000000000000000000000018 401 t: 749 ms

1000000000000000000000020 277 t: 1919 ms

---------------------------------------------`

Die aktuelle Vorgehensweise erscheint bez. Geschwindigkeit nahe am Optimum.

Vielleicht sollte man eher noch mehr Primzahlen dazu nehmen?
Probieren wir 3 ... 101, dann steigt die Zeit bei keinem Divisorfund (z.B. hier bei n-3) an, überwiegend hilft es aber:

`

vor Lookup: nix

vor M.R.: 3 ... 101

1000000000000000000000010 3 t: 62 ms

1000000000000000000000012 5 t: 47 ms

1000000000000000000000014 7 t: 62 ms

1000000000000000000000016 487 t: 952 ms

1000000000000000000000018 11 t: 62 ms

1000000000000000000000020 13 t: 47 ms

`

Daher lasse ich es bei dieser Vorgehensweise:
Vor Lookup nichts, vor Miller-Rabin 3 ... 101.

[TGGC]

Wo ist der Sinn, wenn du nach dem Lookup nochmal das gleiche testest, wie davor schon? Und warum baust du es nicht direkt da ein, wo der % 2 Check ist?

[Erhard Henkes]

@TGGC: Ich werde das nochmal sauber messen. Das Optimum ist sicher auch nicht bei 101, sondern hängt von dem zu überprüfenden Zahlenbereich ab, aber der aktuelle Source ist IMHO recht nahe am Optimum.

Zunächst kommt noch die Miller-Rabin Liars List (2047 ... 18446744066047760377) mit 663.879.564 Bytes ins Spiel. Damit testen wir bei Rückgabe von true aus unserem Milli-Miller-Rabin-Test. Das werde ich beim Start oberhalb der letzten Zahl in der Datei (18446744066047760377) noch abschalten, da das Laden noch recht lange dauert.

@Volkard: Gibt es Liars oberhalb 2^64 zum Download?

Hier zunächst mein aktueller Sourcecode, damit wir harmonisiert bleiben:

#include <iostream>
#include <iomanip>
#include <cstdint>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <vector>
#include <unordered_set>
#include <fstream>
#include "BigIntegerLibrary.hh" // https://mattmccutchen.net/bigint/

using namespace std;

void wait()
{
	cout << "Press any key to continue." << endl;
	cin.clear();
	cin.ignore(numeric_limits<streamsize>::max(), '\n');
	cin.get();
}

unordered_set<uint64_t> setLiar;
uint64_t lastLiar = 0;

uint64_t convertBigToU64(BigUnsigned num)
{
	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);
}

//calculates (a * b) % c
BigUnsigned mulmod(BigUnsigned a, BigUnsigned b, BigUnsigned mod)
{
	BigUnsigned x = 0;
	BigUnsigned y = a % mod;

	while (b > 0)
	{
		if (b % 2 == 1)
		{
			x = (x + y) % mod;
		}
		y = (y * 2) % mod;
		b /= 2;
	}
	return x % mod;
}

BigUnsigned powmod_Volkard(BigUnsigned base, BigUnsigned exp, BigUnsigned modul)
{
	//rechnet 'base hoch exp mod modul'
	BigUnsigned a1 = base, z1 = exp, x = 1, m = modul;
	while (z1 != 0)
	{
		while ((z1 % 2) == 0)
		{
			z1 /= 2;
			a1 = mulmod(a1, a1, m);
		}
		--z1;
		x = mulmod(x, a1, m);
	}
	return x;
}

bool IsProbablyPrime_MillerRabinOptimized(BigUnsigned n)
{
	/*
	example:
	step: "calculate s and d"
	n = 1001
	n-1 = 1000
	s = 3
	d = 125
	1000 / 2^3 = 125

	step: "test with base a"
	ad = a^d mod n  (cf. function powmod)
	...
	*/

	if (n == 2)
		return true;
    // "s" is the logarithm to base base 2 of the biggest number 2^s dividing n-1
    // d = (n-1)/(2^s)
    // calculate d and s
	BigUnsigned d = n - 1;
	BigUnsigned s = 0;
	while ((d & 1) == 0)
    {
        ++s; d >>= 1;
    }

    // We do not use random base a = {2 ... n-1}, but fix it to a = 2
    BigUnsigned a = 2; // base ("witness")
	BigUnsigned ad = (powmod_Volkard(a%n, d, n)); // (a^d) mod n

    if (ad == 1)
        return true; // a^d mod n == 1
	if ((s > 0) && (ad == (n - 1)))
		return true; // a^d mod n == -1 (tests (a^d)^(2^0))

    // test for r = {0 ... s-1}
	for (BigUnsigned r = 1; r < s; ++r)
	{
		// ad is currently ad^(2^(r-1)), thus we square it to get ad^(2^r))
		ad = mulmod(ad, ad, n);
		if (ad == (n - 1))
            return true; // tests (a^d)^(2^(r+1)))
	}
	return 0; // false, composite
}

bool IsPrimeDivisionTest(BigUnsigned number)
{
	if (number<2)    return false;
	if (number == 2)   return true;
	if (number % 2 == 0) return false;
	for (BigUnsigned i = 3; i*i <= number; i += 2)
	{
		if (number%i == 0) return false;
	}
	return true;
}

bool IsPrime(vector<bool>& primes, BigUnsigned number, BigUnsigned primeLimit, unordered_set<uint64_t>& setLiar)
{
	// lookup from bitset (in the RAM)
	if (number <= primeLimit)
		return primes[number.toUnsignedLong()]; 

	// division test with first prime numbers
	if ( 
		 /*(number % 2 == 0) || */ 
		 (number % 3 == 0) || (number % 5 == 0) || (number % 7 == 0) || (number % 11 == 0) || 
		 (number % 13 == 0) || (number % 17 == 0) || (number % 19 == 0) || (number % 23 == 0) || 
		 (number % 29 == 0) || (number % 31 == 0) || (number % 37 == 0) || (number % 41 == 0) || 
		 (number % 43 == 0) || (number % 47 == 0) || (number % 53 == 0) || (number % 59 == 0) ||
		 (number % 61 == 0) || (number % 67 == 0) || (number % 71 == 0) || (number % 73 == 0) || 
		 (number % 79 == 0) || (number % 83 == 0) || (number % 89 == 0) || (number % 97 == 0) || 
		 (number % 101 == 0)
	   )
	{
		return false;
	}

	// Miller-Rabin, not perfect, but fast
	if (!IsProbablyPrime_MillerRabinOptimized(number))
	{
		return false; // 
	}		
	else
	{		
		BigUnsigned bigLastLiar;
		BigUnsigned& refBigLastLiar = bigLastLiar;
		convertU64ToBig(lastLiar, refBigLastLiar);

		if (number > bigLastLiar)
		{
			return true;
		}			
		else
		{
			if (setLiar.find(convertBigToU64(number)) != setLiar.cend())
			{
				// cout << "ALARM: Miller-Rabin-Liar found: " << number << "  ";
				// wait();
				return false; // function tells the truth by failure list (0 ... 2^64), above lastLiar there will be "Prime Hell" again.
			}								
			return true;
		}			
	}
}

bool PrimePairFound(vector<bool>& primes, BigUnsigned i, const BigUnsigned primeLimit, unordered_set<uint64_t>& setLiar)
{
	bool retVal = false;

	BigUnsigned grenze = i / 2;
	BigUnsigned a = 0;
	for (a = 3; a <= grenze; a += 2)
	{
		if (IsPrime(primes, a, primeLimit, setLiar) && IsPrime(primes, i - a, primeLimit, setLiar))
		{
			retVal = true;
			break;
		}
	}

	cout << i << "  " << setw(5) << a << " ";
	return retVal;
}

int main()
{
	BigUnsigned startNumber = stringToBigUnsigned("1000000000000000000000000");
	BigUnsigned endNumber   = stringToBigUnsigned("1000000000000000000000050");

	uint64_t const primeLimit = 2000000; //200000000000;
	cout << "We generate vector<bool>(primes) up to: " << primeLimit << endl; // (this amount depends on the capability of your RAM)

	vector<bool> primes(primeLimit); // processed numbers;
	time_t startTime, endTime;		 // measuring execution time
	uint64_t global_count = 0;		 // total number of primes
	time(&startTime);

	// initialize the number array
	primes[0] = false;
	primes[1] = false;
	for (uint64_t i = 2; i<primeLimit + 1; ++i)
	{
		primes[i] = true;
	}

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

	time(&endTime);

	// count the number of primes
	for (uint64_t i = 0; i<primeLimit + 1; ++i)
	{
		if (primes[i])
			global_count++;
	}

	cout << "In " << difftime(endTime, startTime) << " s we found " << global_count
		<< " prime numbers (2 to " << primeLimit << " )." << endl << endl;

	// Prepare Miller-Rabin "Liar-Liste" and detect Goldbach prime pairs with Miller-Rabin test
	ifstream fileLiar("F:\\Daten\\Miller_Rabin_Liar\\2-SPRP-2-to-64.txt", fstream::in); // https://miller-rabin.appspot.com/
	cout << "Miller-Rabin Liar file starts to be read from file." << endl;
	uint64_t count=0;
	uint64_t u64LiarNumber=0;
	while (!fileLiar.eof())
	{
		fileLiar >> u64LiarNumber;
		// cout << u64LiarNumber << "\t";
		setLiar.insert(u64LiarNumber);
		count++;
	}
	fileLiar.close();
	lastLiar = u64LiarNumber;
	cout << "Miller-Rabin Liar file with " << count << " liars is read." << endl << "Last liar: " << u64LiarNumber << endl << endl;

	unordered_set<uint64_t>& refSetLiar = setLiar;

	cout << "Now Goldbach first prime number pair will be detected:" << endl;

	clock_t zeit1, zeit2, zeit2_old;
	zeit2 = zeit1 = clock();

	for (BigUnsigned i = startNumber; i <= endNumber; i += 2)
	{
		if (!PrimePairFound(primes, i, BigUnsigned((unsigned long)primeLimit), refSetLiar))
		{
			cout << "Counterevidence found for this number: " << i << endl;
			wait();
		}

		zeit2_old = zeit2;
		zeit2 = clock();

		cout << "time: " << 1000 * (zeit2 - zeit2_old) / CLOCKS_PER_SEC << " ms" << endl;
	}
	wait();
}

IMHO wäre nun MT oder verstärkte Nutzung des RAM dran.

`We generate vector<bool>(primes) up to: 2000000

In 0 s we found 148933 prime numbers (2 to 2000000 ).

Miller-Rabin Liar file starts to be read from file.

Miller-Rabin Liar file with 31894015 liars is read.

Last liar: 18446744066047760377

Now Goldbach first prime number pair will be detected:

1000000000000000000000000 257 time: 448 ms

1000000000000000000000002 349 time: 758 ms

1000000000000000000000004 307 time: 409 ms

1000000000000000000000006 263 time: 374 ms

1000000000000000000000008 311 time: 643 ms

1000000000000000000000010 3 time: 32 ms

1000000000000000000000012 5 time: 31 ms

1000000000000000000000014 7 time: 32 ms

1000000000000000000000016 487 time: 480 ms

1000000000000000000000018 11 time: 32 ms

1000000000000000000000020 13 time: 31 ms

1000000000000000000000022 613 time: 643 ms

1000000000000000000000024 17 time: 61 ms

1000000000000000000000026 19 time: 62 ms

1000000000000000000000028 331 time: 343 ms

1000000000000000000000030 23 time: 61 ms

1000000000000000000000032 379 time: 852 ms

1000000000000000000000034 337 time: 477 ms

1000000000000000000000036 29 time: 92 ms

1000000000000000000000038 31 time: 91 ms

1000000000000000000000040 631 time: 809 ms

1000000000000000000000042 389 time: 505 ms

1000000000000000000000044 37 time: 122 ms

1000000000000000000000046 349 time: 403 ms

1000000000000000000000048 41 time: 93 ms

1000000000000000000000050 43 time: 153 ms`

Im Bereich 10^24 geht es also recht flott ab. Hier mein Benchmark bei 10^120:

`

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 173 time: 30615 ms

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002 709 time: 126508 ms

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004 613 time: 71922 ms

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000006 179 time: 34387 ms

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000008 181 time: 37669 ms

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010 619 time: 92918 ms`

[Erhard Henkes]

@Volkard, Bengo, Arcoth: Habe ich alle wichtigen Ideen zur Beschleunigung inkorporiert oder habe ich noch etwas Wichtiges übersehen? Was könnte man sinnvoll parallel laufen lassen? Zehn Zahlen gleichzeitig testen?

Die Klasse BigInteger bietet auch einen Algo für modexp:

BigUnsigned modexp(const BigInteger &base, const BigUnsigned &exponent,	const BigUnsigned &modulus) 
{
	BigUnsigned ans = 1, base2 = (base % modulus).getMagnitude();
	BigUnsigned::Index i = exponent.bitLength();
	// For each bit of the exponent, most to least significant...
	while (i > 0) 
	{
		i--;
		// Square.
		ans *= ans;
		ans %= modulus;
		// And multiply if the bit is a 1.
		if (exponent.getBit(i)) 
		{
			ans *= base2;
			ans %= modulus;
		}
	}
	return ans;
}

[volkard]

Erhard Henkes schrieb:

@Volkard: Gibt es Liars oberhalb 2^64 zum Download?

Ich kenne keine.

[TGGC]

Wie gesagt, wuerde ich es nochmal mit einer loopenden Tabelle pobieren, die in den Cache passt um die Teilbarkeit durch die ersten Primzahlen zu pruefen. Das wuerde dann mehrere Divisionen durch eine Division und ein Lookup ersetzen. Gleichzeitig kann man das dann eben auch nutzen um den primes vector viel kleiner zu machen.

[Erhard Henkes]

Ab dem Bereich 10^240 sind wir im Bereich ca. 1h/Geradzahl angelangt:

10^240, erste Primzahl 2351, t: 2215 s
10^241, erste Primzahl 2069, t: 3929 s

Hier findet man eine erste zum Paar passende Primzahl erst über 1000.

[Bengo]

Das sind etwa 800bit, scheint mir so, als ob man das noch sehr optimieren kann, da mein Computer auch 4096 bit Primzahlen generieren kann.(In weniger als einer Stunde) bei SSL

[volkard]

Bengo schrieb:

Das sind etwa 800bit, scheint mir so, als ob man das noch sehr optimieren kann, da mein Computer auch 4096 bit Primzahlen generieren kann.(In weniger als einer Stunde) bei SSL

Welchen Zweck hätte eine 4096-Bit-Primzahl bei SSL? Sicher, daß Du da 4096-Bit-Primzahlen berechnen läßt?

Naja, SEHR optimieren nenne ich das nicht. Klar handoptimiertes asm an manche Stellen, MMX/SSE, Multitreading, damit kann man noch was holen. Mit der riesigen Gefahr, dass man sich dabei verzettelt und einen läppischen Faktor 10 kriegt und dafür Einbußen an der Einfachheit und Wartbarkeit hat und deswegen am Algo spart und billiardenfach lahmer ist.

[Erhard Henkes]

Bitte nicht übersehen, dass wir bei 10^241 mit der ersten Primzahl 2069 in den 3929 Sekunden bei allen Primzahlen von 3 bis 2069 alle entsprechenden Partner "Geradzahl minus erste Primzahl" auf ihre Primzahlfähigkeit prüfen. Im Lookup findet da schon lange nichts mehr, man ist somit auf den kleinen Divisionstest (z.Z. 3 ... 101) angewiesen und auf den Milli-Miller-Rabin. Den Lügnertest kann man natürlich gar nicht erst laden und abschalten über 2^64.

Ich werde im nächsten Schritt das C++ 11 std::thread einbauen und verschiedene parallele Strategien testen.

Beispiele:
- mehrere Geradzahlen parallel auf Goldbach-Vermutung prüfen
- mehrere erste Primzahlen parallel auf ihren möglichen Primzahl-Partner prüfen (hier findet man aber nicht mehr die erste kleine Primzahl, sondern irgendein Paar)

[Bengo]

volkard schrieb:

Bengo schrieb:

Das sind etwa 800bit, scheint mir so, als ob man das noch sehr optimieren kann, da mein Computer auch 4096 bit Primzahlen generieren kann.(In weniger als einer Stunde) bei SSL

Welchen Zweck hätte eine 4096-Bit-Primzahl bei SSL? Sicher, daß Du da 4096-Bit-Primzahlen berechnen läßt?

Naja, SEHR optimieren nenne ich das nicht. Klar handoptimiertes asm an manche Stellen, MMX/SSE, Multitreading, damit kann man noch was holen. Mit der riesigen Gefahr, dass man sich dabei verzettelt und einen läppischen Faktor 10 kriegt und dafür Einbußen an der Einfachheit und Wartbarkeit hat und deswegen am Algo spart und billiardenfach lahmer ist.

Stimmt hast recht, das RSA Modul hat 4096bit, die Primzahlen dann aber immernoch mehr als 2048bit.

[Erhard Henkes]

Die Thread-Unterstützung habe ich mal eingebaut (noch nicht ideal) und im Bereich 10^120 bis 10^120 + 20 getestet (11 Geradzahlen, 11 Threads). Das läuft gut, man er hält die Ergebnisse in der gleichen Zeit wie vorher die langsamste Zahl in diesem Bereich. 11 Threads sind bei mir mit 12 virtuellen Kernen optimal. Hier das Programm, falls es jemand ausprobieren/optimieren will:

#include "stdafx.h"  // MS VC++  

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

#include <iostream>
#include <iomanip>
#include <cstdint>
#include <cstdlib>
#include <cassert>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <vector>
#include <unordered_set>
#include <fstream>
#include <thread>
#include <mutex>
#include "BigIntegerLibrary.hh" // https://mattmccutchen.net/bigint/

//#define MILLER_RABIN_LIAR_CHECK

using namespace std;

void setConsole() // windows.h
{
	_CONSOLE_SCREEN_BUFFER_INFO info;
	HANDLE handle = GetStdHandle(STD_OUTPUT_HANDLE);
	GetConsoleScreenBufferInfo(handle, &info);
	COORD c;
	c.X = std::max<SHORT>(info.dwSize.X, 150);
	c.Y = std::max<SHORT>(info.dwSize.Y, 1000);
	SetConsoleScreenBufferSize(handle, c);

	SMALL_RECT RECT;
	RECT.Left = 0;
	RECT.Top = 0;
	RECT.Bottom = std::max(info.srWindow.Bottom - info.srWindow.Top, 40 - 1);
	RECT.Right = std::max(c.X - 1, 100 - 1);
	SetConsoleWindowInfo(handle, TRUE, &RECT);
}

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

void wait()
{
	cout << "Press any key to continue." << endl;
	cin.clear();
	cin.ignore(numeric_limits<streamsize>::max(), '\n');
	cin.get();
}

mutex mutex_; // for cout 

unordered_set<uint64_t> setLiar;
uint64_t lastLiar = 0;

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

//calculates (a * b) % c
BigUnsigned mulmod(BigUnsigned a, BigUnsigned b, BigUnsigned mod)
{
	BigUnsigned x = 0;
	BigUnsigned y = a % mod;

	while (b > 0)
	{
		if (b % 2 == 1)
		{
			x = (x + y) % mod;
		}
		y = (y * 2) % mod;
		b /= 2;
	}
	return x % mod;
}

BigUnsigned powmod_Volkard(BigUnsigned base, BigUnsigned exp, BigUnsigned modul)
{
	//rechnet 'base hoch exp mod modul'
	BigUnsigned a1 = base, z1 = exp, x = 1, m = modul;
	while (z1 != 0)
	{
		while ((z1 % 2) == 0)
		{
			z1 /= 2;
			a1 = mulmod(a1, a1, m);
		}
		--z1;
		x = mulmod(x, a1, m);
	}
	return x;
}

bool IsProbablyPrime_MillerRabinOptimized(BigUnsigned n)
{
	/*
	example:
	step: "calculate s and d"
	n = 1001
	n-1 = 1000
	s = 3
	d = 125
	1000 / 2^3 = 125

	step: "test with base a"
	ad = a^d mod n  (cf. function powmod)
	...
	*/

	if (n == 2)
		return true;
	// "s" is the logarithm to base base 2 of the biggest number 2^s dividing n-1
	// d = (n-1)/(2^s)
	// calculate d and s
	BigUnsigned d = n - 1;
	BigUnsigned s = 0;
	while ((d & 1) == 0)
	{
		++s; d >>= 1;
	}

	// We do not use random base a = {2 ... n-1}, but fix it to a = 2
	BigUnsigned a = 2; // base ("witness")
	BigUnsigned ad = (powmod_Volkard(a%n, d, n)); // (a^d) mod n

	if (ad == 1)
		return true; // a^d mod n == 1
	if ((s > 0) && (ad == (n - 1)))
		return true; // a^d mod n == -1 (tests (a^d)^(2^0))

					 // test for r = {0 ... s-1}
	for (BigUnsigned r = 1; r < s; ++r)
	{
		// ad is currently ad^(2^(r-1)), thus we square it to get ad^(2^r))
		ad = mulmod(ad, ad, n);
		if (ad == (n - 1))
			return true; // tests (a^d)^(2^(r+1)))
	}
	return 0; // false, composite
}

bool IsPrimeDivisionTest(BigUnsigned number)
{
	if (number<2)    return false;
	if (number == 2)   return true;
	if (number % 2 == 0) return false;
	for (BigUnsigned i = 3; i*i <= number; i += 2)
	{
		if (number%i == 0) return false;
	}
	return true;
}

bool IsPrime(vector<bool>& primes, BigUnsigned number, BigUnsigned primeLimit, unordered_set<uint64_t>& setLiar)
{
	// lookup from bitset (in the RAM)
	if (number <= primeLimit)
		return primes[number.toUnsignedLong()];

	// division test with first prime numbers
	if (
		/*(number % 2 == 0) || */
		(number % 3 == 0) || (number % 5 == 0) || (number % 7 == 0) || (number % 11 == 0) ||
		(number % 13 == 0) || (number % 17 == 0) || (number % 19 == 0) || (number % 23 == 0) ||
		(number % 29 == 0) || (number % 31 == 0) || (number % 37 == 0) || (number % 41 == 0) ||
		(number % 43 == 0) || (number % 47 == 0) || (number % 53 == 0) || (number % 59 == 0) ||
		(number % 61 == 0) || (number % 67 == 0) || (number % 71 == 0) || (number % 73 == 0) ||
		(number % 79 == 0) || (number % 83 == 0) || (number % 89 == 0) || (number % 97 == 0) ||
		(number % 101 == 0)
		)
	{
		return false;
	}

	// Miller-Rabin, not perfect, but fast
	if (!IsProbablyPrime_MillerRabinOptimized(number))
	{
		return false; //
	}
	else
	{
		if (lastLiar == 0)
		{
			return true;
		}

		BigUnsigned bigLastLiar;
		BigUnsigned& refBigLastLiar = bigLastLiar;
		convertU64ToBig(lastLiar, refBigLastLiar);

		if (number > bigLastLiar)
		{
			return true;
		}
		else
		{
			if (setLiar.find(convertBigToU64(number)) != setLiar.cend())
			{
				// cout << "ALARM: Miller-Rabin-Liar found: " << number << "  ";
				// wait();
				return false; // function tells the truth by failure list (0 ... 2^64), above lastLiar there will be "Prime Hell" again.
			}
			return true;
		}
	}
}

void PrimePairFound(vector<bool>& primes, BigUnsigned i, const BigUnsigned primeLimit, unordered_set<uint64_t>& setLiar, bool* retVal)
{
	*retVal = false;

	BigUnsigned grenze = i / 2;
	BigUnsigned a = 0;
	for (a = 3; a <= grenze; a += 2)
	{
		if (IsPrime(primes, a, primeLimit, setLiar) && IsPrime(primes, i - a, primeLimit, setLiar))
		{
			*retVal = true;
			break;
		}
	}

	mutex_.lock();
	cout << i << "  " << setw(5) << a << endl;
	mutex_.unlock();
}

int main()
{
	setConsole(); // windows.h // big console for large numbers! ^^
	textcolor(14); // windows.h // 10 = light green, 14 = light yellow

	BigUnsigned startNumber = stringToBigUnsigned("1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
	BigUnsigned endNumber = stringToBigUnsigned("1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000020");
	const int numberOfThreads = 11; // count of even numbers from start to end = (end - start)/2 + 1

	uint64_t const primeLimit = 2000000; //200000000000;
	cout << "We generate vector<bool>(primes) up to: " << primeLimit << endl; // (this amount depends on the capability of your RAM)

	vector<bool> primes(primeLimit); // processed numbers;
	time_t startTime, endTime;       // measuring execution time
	uint64_t global_count = 0;       // total number of primes
	time(&startTime);

	// initialize the number array
	primes[0] = false;
	primes[1] = false;
	for (uint64_t i = 2; i<primeLimit + 1; ++i)
	{
		primes[i] = true;
	}

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

	time(&endTime);

	// count the number of primes
	for (uint64_t i = 0; i<primeLimit + 1; ++i)
	{
		if (primes[i])
			global_count++;
	}

	cout << "In " << difftime(endTime, startTime) << " s we found " << global_count
		<< " prime numbers (2 to " << primeLimit << " )." << endl << endl;

#ifdef MILLER_RABIN_LIAR_CHECK
	// Prepare Miller-Rabin "Liar-Liste" and detect Goldbach prime pairs with Miller-Rabin test
	ifstream fileLiar("F:\\Daten\\Miller_Rabin_Liar\\2-SPRP-2-to-64.txt", fstream::in); // https://miller-rabin.appspot.com/
	cout << "Miller-Rabin Liar file starts to be read from file." << endl;
	uint64_t count = 0;
	uint64_t u64LiarNumber = 0;
	while (!fileLiar.eof())
	{
		fileLiar >> u64LiarNumber;
		// cout << u64LiarNumber << "\t";
		setLiar.insert(u64LiarNumber);
		count++;
	}
	fileLiar.close();
	lastLiar = u64LiarNumber;
	cout << "Miller-Rabin Liar file with " << count << " liars is read." << endl << "Last liar: " << u64LiarNumber << endl << endl;
#endif // MILLER_RABIN_LIAR_CHECK

	unordered_set<uint64_t>& refSetLiar = setLiar;

#ifdef MILLER_RABIN_LIAR_CHECK
	cout << "Now Goldbach first prime number pair will be detected:" << endl;
#else
	cout << "Now Goldbach probable (no liar list available) first prime number pair will be detected:" << endl;
#endif

	clock_t zeit1, zeit2, zeit2_old, delta;
	zeit2 = zeit1 = clock();
	uint64_t sumTime = 0;

	for (BigUnsigned i = startNumber; i <= endNumber; i += (2 * numberOfThreads))
	{		
		bool retVal[numberOfThreads];		
		vector<thread> t;

		for (int x = 0; x<numberOfThreads; ++x)
		{
			retVal[x] = false;		
			if ((i+2*x) <= endNumber)
				t.emplace_back(PrimePairFound, primes, i+2*x, BigUnsigned((unsigned long)primeLimit), refSetLiar, retVal + x);
		}

		for (int x = 0; x < numberOfThreads; ++x)
		{
			if ((i + 2 * x) <= endNumber)
			{
				t[x].join(); // we have to wait for the result of the thread!
				if (retVal[x] == false)
					cout << "Counterevidence found for this number: " << i + 2 * x << endl;							
			}				
		}

		zeit2_old = zeit2;
		zeit2 = clock();
		delta = zeit2 - zeit2_old;
		sumTime += delta;
		cout << "time: " << 1000 * delta / CLOCKS_PER_SEC << " ms" << endl << endl;
	}
	cout << "total time: " << sumTime/1000 << " s" << endl << endl;
	wait();
}

`We generate vector<bool>(primes) up to: 2000000

In 0 s we found 148933 prime numbers (2 to 2000000 ).

Now Goldbach probable (no liar list available) first prime number pair will be detected:

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000018 191

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 173

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000006 179

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000008 181

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000016 733

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004 613

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000020 193

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010 619

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012 719

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002 709

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000014 1153

time: 237940 ms

total time: 237 s`

[Bengo]

Erhard Henkes schrieb:

Die Thread-Unterstützung habe ich mal eingebaut (noch nicht ideal) und im Bereich 10^120 bis 10^120 + 20 getestet (11 Geradzahlen, 11 Threads). Das läuft gut, man er hält die Ergebnisse in der gleichen Zeit wie vorher die langsamste Zahl in diesem Bereich. 11 Threads sind bei mir mit 12 virtuellen Kernen optimal. Hier das Programm, falls es jemand ausprobieren/optimieren will:

Hab nur 2 Kerne

[Erhard Henkes]

Hab nur 2 Kerne

Dann stellst Du die numberOfThreads auf 2, 3 oder 4 (dein OS kümmert sich um das MT Scheduling) ein. Ich baue das noch um, sodass die richtige Threadanzahl verwendet wird und diese auch immer etwas zu schaffen haben.
Probiere momentan verschiedene Abläufe aus.

Wichtige Frage: Kann man die Gesamtauslastung (vom Prg aus) des gesamten Programms auf einen gewissen Wert (z.B. 50 oder 80%) einstellen?