Goldbach Vermutung

Erhard Henkes

Ich habe am Beispiel 1^120 einen erweiterten Divisions-Precheck (3 ... 199 anstelle 3 ... 101) ausprobiert:

// 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) || (number % 103 == 0) || (number % 107 == 0) || (number % 109 == 0) ||
		(number % 113 == 0) || (number % 127 == 0) || (number % 131 == 0) || (number % 137 == 0) ||
		(number % 139 == 0) || (number % 149 == 0) || (number % 151 == 0) || (number % 157 == 0) ||
		(number % 163 == 0) || (number % 167 == 0) || (number % 173 == 0) || (number % 179 == 0) ||
		(number % 181 == 0) || (number % 191 == 0) || (number % 193 == 0) || (number % 197 == 0) ||
		(number % 199 == 0) 
		)
	{
		return false;
	}

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

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000018 191

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000008 181

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000006 179

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 173

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004 613

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000016 733

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000020 193

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010 619

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012 719

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002 709

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000014 1153

time: 184750 ms

total time: 184 s`

Eine deutliche Beschleunigung (184 s vs. 237 s mit 3 ... 101). Es verändert auch die Reihenfolge, in der die ersten Primzahlen als Teil einer Primzahlensumme auftauchen. Wo hier das Optimum in Abhängigkeit von der zu prüfenden Zahl liegt, ist eine interessante Frage.

Da wir in unser vector<bool> bereits die Primzahleninformationen im unteren Zahlenbereich gespeichert haben, können wir diese Information im Code nutzen:

// division test with first prime numbers
	const uint32_t DIVISOR_PRECHECK_MAX = 199;
	for (uint32_t divisor = 3; divisor <= DIVISOR_PRECHECK_MAX; ++divisor)
	{
		if (primes[divisor])
		{
			if (number % divisor == 0)
				return false;
		}
	}

Der Lookup-Vorgang kostet bei unserem Test im Bereich 1^120 in Summe ca. drei Sekunden (187 vs. 184 s). Man kann auf diese Weise aber leichter das Geschwindigkeits-Optimum (Welcher DIVISOR_PRECHECK_MAX in Abhängigkeit von startNumber?) finden.

Folgende kleine Serie habe ich getestet:
10^120:

DIVISOR PRECHECK von ... bis ...: Zeit zur Prüfung der 11 Ganzzahlen von 10^120 bis incl. 10^120 + 20 mit 11 parallelen Threads (langsamste Testschleife mit Milli-Rabin-isProbablyPrime ist zeitbestimmend):

`3 ... 101: 226 s

3 ... 199: 187 s

3 ... 997: 142 s

3 ... 9973: 98 s

3 ... 99991: 81 s

3 ... 999983: 79 s

3 ... 1999993: 84 s`

Fazit: Für Goldbach-Primzahlenpaar-Tests im Bereich 1^120 und höher kann ich 999983 als DIVISOR_PRECHECK_MAX vorteilhaft einsetzen.

Erhard Henkes

Ich sehe bezüglich weiterer Optimierungen folgende Felder:
- besseres BigInteger (vlt. doch gmp)
- voll ausgelastetes threading
- mulmod in asm <-- da hängt allerdings BigInteger dran
- die "loopende Tabelle im Cache" (habe die genaue Idee von TGGC noch nicht verstanden, Code wäre hilfreich)
- weitere mathematische Tricks

Ich hab gerade mal was mit MPIR zusammengeklatsch. Das ist überhaupt kein Vergleich:

1...000    173
1...002    709
1...004    613
1...006    179
1...008    181
1...010    619
1...012    719
1...014   1153
1...016    733
1...018    191
1...020    193
time: 58 ms (Ja, wirklich ms)

Und das ist noch nichtmal Multithreaded... Hier der Code:

#include <gmp.h>
#include <gmpxx.h>
#include <iostream>
#include <iomanip>

using namespace std;

gmp_randstate_t rnd;

bool PrimePairFound(const mpz_class& i)
{
  bool retVal = false;

  mpz_class grenze = i / 2;
  mpz_class a = 0;
  for(a = 3; a <= grenze; a += 2)
  {
    mpz_class b = i - a;
    if(mpz_likely_prime_p(a.get_mpz_t(), rnd, 2) && mpz_likely_prime_p(b.get_mpz_t(), rnd, 2))
    {
      retVal = true;
      break;
    }
  }

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

int main()
{
  gmp_randinit_default(rnd);
  mpz_class startNumber("1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
  mpz_class endNumber("1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000020");

  clock_t start = clock();
  for(mpz_class i = startNumber; i <= endNumber; i += 2)
  {
    if(!PrimePairFound(i))
    {
      cout << "Counterevidence found for this number: " << i << endl;
    }
  }
  clock_t end = clock();
  cout << "time: " << 1000*(end - start) / CLOCKS_PER_SEC << " ms" << endl;
  cin.get();
}

Die Funktion für den Primzahltest hat zwar ein likely im Namen aber wenn ich es richtig verstanden habe ist es sehr sehr unwahrscheinlich eine Zahl als Primzahl zu erkennen wenn es eigentlich keine ist.

Bengo

Wikipedia sagt, der Test ist bei einer zufälligen Basis zu 25% nicht zuverlässig, skaliert aber sehr gut, nach 4 Tests sind nur noch 0,4%, nach 10 Tests, etwa 1ppm.

Oder man kann auch die mpz_probab_prime_p Funktion nutzen, bei der man die Anzahl der Miller-Rabin-Test angeben kann. Wenn ich dort 100 Tests angebe braucht der Benchmark 112ms, dafür werden zusammengesetzte Zahlen mit einer Wahrscheinlichkeit von maximal 2^-200 ≈ 6,22*10^-61 als Primzahlen erkannt. Das ist schon ziemlich unwahrscheinlich...

Erhard Henkes

Schwächen bez. Geschwindigkeit sind die auf uint32_t Blöcken aufbauende Klasse BigInteger und die C-Funktion mulmod. Der Code bewegt sich nach meinen Messungen mit QueryPerformanceCounter zu ca. 70% in mulmod, das mit dieser Klasse hantiert.

Aktueller Sourcecode: http://codepad.org/puR6GPwG

TGGC

Erhard Henkes schrieb:

die "loopende Tabelle im Cache" (habe die genaue Idee von TGGC noch nicht verstanden, Code wäre hilfreich)

Ist eigentlich ganz simpel. Angenommen du willst Teilbarkeit durch 2 testen mit einer Tabelle:

bool table2[2] = {true, false};

isDivBy2(uint i)
{
    return table2[i%2];
}

Sollte klar sein, wie sich das fortsetzt fuer andere Zahlen

isDivBy2or3or5(uint i)
{
    return table2[i%2] || table2[i%3] || table2[i%5];
}

Jetzt kommt der Witz, das auch die Ausgabe dieser Funktion sich nach 2 * 3 * 5 Schritten wiederholt, man kann das also wieder in eine neue Tabelle fuellen, danach hat man:

isDivBy2or3or5Combined(uint i)
{
    return table2or3or5[i%(2*3*5)];
}

Das kann man beliebig erweitern und so auf Kosten von Speicher Divisionen sparen. Im Grunde Sieb des Erastothenes, in der aber nur die ersten n Primzahlen eingetragen sind. Bei vielen Zahlen steht mit dem Test ziemlich schnell fest, das ein Teiler in den ersten n Primzahen vorhanden ist, d.h. ein Eintrag in deinem riesigen bit vector ist dann gar nicht mehr noetig. table2or3or5 haette ja eine Laenge von 30. Wenn dort z.b. nur 10 mal true drin steht, so brauch der vector auch nur fuer diese 10 Werte ein Bit enthalten, die anderen 20 Bit brauchst du gar nicht erst speichern.

Finnegan

Erhard Henkes schrieb:

Schwächen bez. Geschwindigkeit sind die auf uint32_t Blöcken aufbauende Klasse BigInteger und die C-Funktion mulmod. Der Code bewegt sich nach meinen Messungen mit QueryPerformanceCounter zu ca. 70% in mulmod, das mit dieser Klasse hantiert.

Da fällt mir grad auf: Hab mir mulmod nicht soo genau angesehen, aber so wie ich das sehe wird dort lediglich in Zeile 8 das x inklusive Modulo neu gesetzt und das mod bleibt unverändert.
Braucht man da das letzte Modulo beim return überhaupt noch? Denke das ist als BigInt noch ein bisschen teurer als ein einzelner 32bit int , oder?

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;

Ebenso: Ich weiss nicht wie gut die BigInt-Lib optimiert ist, ich weiss dass man sowas bei C++-Primitiven nicht machen braucht, aber das hier könnte was anderes sein:

1. Test auf "b ist ungerade" in Zeile 6 kann man eventuell vereinfachen zu "Ist niedrigstwertiges Bit 1?". Wenn die Lib nach Adam Riese rechnet dann ist das ne volle Division durch 2 die er mit dem Muldulo sonst macht.
2. Division durch 2 kann man als Bit-Shift um 1 nach rechts darstellen (eventuell flotter, bei nem C++-Compiler und Primitiven würde ich erwarten dass er das optimiert, aber bei den BigInts bin ich mir nicht so sicher).
3. Multiplikation mit 2 ebenfalls: Bit-Shift links.
4. Weiss nicht ob die temporären BigInts die da erzeugt werden nicht teurer sind, als wenn man die Operationen in-place macht, also z.B. statt x = (x + y) % mod;

x += y;
  x %= mod;

Wie gesagt, bei Compiler-Primitiven wäre das eher sinnlos bis kontraproduktiv, aber bei Bigints könnte es was bringen... muss mir mal anschauen wie die Lib implementiert ist.

Gruss,
Finnegan

P.S.: bigint-Lib kurz überflogen:
1. Keine flache Datenstruktur. Bigints werden auf dem Heap erzeugt -> lokale Bigints nicht im "heißen" Cache (Stack).
2. Keine Move-Konstruktoren/Assignment -> Temporäre Bigints erfordern je ein zusätzliches new + memcpy + delete[] -> in-place-Operationen, also via +=, %=, et al. müssten eigentlich flotter sein.
3. Keine Modulo/Divisions/Multiplikations-Optimierung für "Power-Of-Two" -> Wenn geht Shift-Operatoren verwenden, besonders bei div/mul mit 2. Und Test auf ungerade sollte tatsächlich via getBit()
schneller sein. Ich schau mal ob ich kurz einen Alternativvorschlag zusammenschustere

P.P.S.: Könntest es mal damit probieren (hoffe ich habe nix verbockt!). Habe versucht so weit wie möglich die relativ teuren Kopien zu vermeiden (pre-C++11-oldschool):

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

    while (b > 0)
    {
        if (b.getBit(0))
        {
            x += y;
            x %= mod;
        }
        y <<= 1;
        y %= mod;
        b >>= 1;
    }
    return x;
}

Wäre zwar vielleicht besser direkt ie Bigint-Lib anzupassen und ihr das "Moven" beizubringen und ein paar andere Sachen, aber das würde wohl zu weit ausholen

Erhard Henkes

Das hier läuft schon deutlich fixer:

BigUnsigned mulmod(BigUnsigned a, BigUnsigned b, BigUnsigned mod)
{
	/////////////
	QueryPerformanceCounter((LARGE_INTEGER*)&startCount);
	/////////////
	BigUnsigned x = 0;
	BigUnsigned y = a % mod;

	while (b > 0)
	{
		if (((uint32_t)b.getBlock(0)) & 1)
		{
			x += y;
			x %= mod;			
		}

		y <<= 1;
		y %= mod;

		b >>= 1;		
	}
	/////////////
	QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
	sumTime_Function += (lastCount - startCount);
	/////////////
	return x;
}

Finnegan

Erhard Henkes schrieb:

Das hier läuft schon deutlich fixer:

Cool-O-Matik!

Gibt wahrscheinlich noch ne Reihe andere Stellen im Code, die davon profitieren würden. Ist nur ein bissl blöd dass das der Code-Lesbarkeit schadet. Das ist für mich eher so C-Hacker-Stil von jemandem der grad von Assembler kommt

Erhard Henkes

Deine Variante ist 8% schneller, daher baue ich sie gerne ein. Insgesamt ein großer Schritt.

Also neuer Sourcecode: http://codepad.org/ReYuzGjk (nur mulmod)

Erhard Henkes

... auch an anderen Stellen verändert: http://codepad.org/f5zbHebm

@Finnegan: Danke!

mulmod vorher: 70% jetzt: 45 % (bei 10^54)

Es lohnt sich also, weiter an dem Interface BigInteger und mulmod zu tunen.

Finnegan

Erhard Henkes schrieb:

mulmod vorher: 70% jetzt: 45 % (bei 10^54)

Gern. Schön wenn man so ne Funktion hat, die im Profiler leuchtet wie ein Weihnachtsbaum. Da ist das Optimieren leicht.
Leider hat man auch oft mit Programmen zu tun wo alle Funktionen irgendwie gleich "langsam" sind, das ist immer ein Krampf.
...aber bei so ner schönen fetten spanischen Galleone schraub ich mir glatt mein Holzbein an und leg' ne Überstunde ein

Erhard Henkes

Vielleicht könnte man auch eine eigene MyBigInteger-Klasse auf Basis uint64_t anstelle uint32_t basteln, nur mit den Funktionen, die man hier braucht. Dann könnte man das auf Geschwindigkeit genau dafür trimmen.

Jodocus

Was ist das Ziel des ganzen hier? Willst du so weit kommen, wie du nur kannst? Dann würde ich jetzt eher aufhören, zu mikrooptimieren, sondern mir die BigInt-Implementierung ansehen. BigInts gibt's viele, aber richtig gute nicht. Wenn du's selber machen willst, dann würde ich sie auf doppelter Wortbreite basieren lassen und ebenfalls im Zweierkomplement darstellen.

Die Addition muss mit ziemlicher Sicherheit in Assembler gemacht werden, da es dort Mnemonics wie ADC gibt, die einem das explizite Carryen sparen (sofern der Compiler es nicht schon so generiert).

Danach Multiplikation: da kannst du mit Karazuba oder gleich mit Toom-Cook anfangen. Oder du machst vorher etwas Mathe, lädst dir die Fast Fourier Transform runter und implementierst Schönhage-Strassen. Die Algorithmen lassen sich auch allesamt parallelisieren + cache-optimieren.

DJohn

Hallo

Erhard Henkes schrieb:

Ich sehe bezüglich weiterer Optimierungen folgende Felder:
...

Wenn ich den Code richtig verstehe, beginnt Ihr bei jeder Zahl wieder vorn mit den Test potentieller Primzahlpärchen. Bei den großen Zahlen wie 10^20 oder sogar 10^240 sind die gefundenen großen Primzahlen doch viel zu wertvoll, um sie gleich wieder wegzuwerfen.

1.) Wenn ich z.B. für x1 = 13 + (x1-13) gefunden habe, weil x1-13 eine Primzahl ist, dann kann ich für alle Primzahlen größer 13 aus meiner Primzahltabelle weiter PrimePairs bilden und die aus der Liste der zu prüfenden Zahlen streichen.
Also:
x3 = (x1+4) = 17 + (X1-13)
x4 = (x1+6) = 19 + (X1-13)
u.s.w.
2.) für x2 = (x1+2) passt das noch nicht, aber auch hier gibt es eine Optimierung:
beim Prüfen von x1 habe ich für x1-3, x1-5, x1-7, x1-11 schon festgestellt, dass das keine Primzahlen sind. Für x2 verschiebt sich das um 2, somit kann ich mir ein paar Primzahltests gleich sparen: von x2-5, x2-7, x2-13 weiß ich schon, dass das keine Primzahlen sind. Hierfür kann man sich intern eine Tabelle halten für alle ungeraden Zahlen von endNumber-3 (Index 0) bis startNumber-n mit jeweils eine Flag für "Prime", "not Prime", oder "Unbestimmt" (evtl. noch ein "wird gerade berechnet" fürs parallelisieren). Natürlich muss die Tabelle irgendwie begrenzt werden, aber ich schätze mal, dass man die meisten PrimPairs findet, bevor die Tabelle unhandlich groß wird.

Finnegan

Erhard Henkes schrieb:

Vielleicht könnte man auch eine eigene MyBigInteger-Klasse auf Basis uint64_t anstelle uint32_t basteln, nur mit den Funktionen, die man hier braucht. Dann könnte man das auf Geschwindigkeit genau dafür trimmen.

Ich würde so frech sein und behaupten, dass der Compiler das schon macht soweit er kann. Zumindest bei dem was ist schon in von Visual C++ generiertem Maschinencode gesehen habe, würde ich sogar davon ausgehen, dass er, wenn z.B. die entsprechende Architektur aktiviert ist (bspw. -arch:avx oder -arch:avx2 - SSE2 ist per default an), er sogar schon entsprechende Vektor-Instruktionen erzeugt... aber probieren geht über studieren. Hab zwar jetzt keine Zeit mehr dafür, aber ein kurzes Mini-Testprogramm und dann im Debug-Modus mit Code-Optimierungen und mit Breakpoint auf die Operation auf "Show Disassembly" gehen, sollte Aufschluss geben (xmm/ymm-Register und vpadd*-Instruktionen wären z.B. bei Addition ein Hinweis... wichtig bei so einem Test: Variablen mit denen man sowas testet sollten nicht durch irgendwelche Konstanten/Literale gefüllt werden, sonst rechnet der das Ergebnis einfach beim Kompilieren aus und baut es direkt in die EXE ein [Stichw. Konstantenfaltung] ... ich hab mir da immer mit dem volatile -Keyword ausgeholfen wenn ich wirkliche Instruktionen sehen wollte).

Falls man es wirklich selbst implementieren will, könnte es allerdings mehr bringen, dem CPU-Cache besser entgegenzukommen. Da muss man allerdings einiges an Hirnschmalz reinstecken, auch wenn die Grundregel simpel ist: Möglichst linear auf den Speicher zugreifen und keine wilden Pointer-Sprünge quer durch den RAM. Daher habe ich auch vorhin bemerkt dass die BigInts alle auf dem Heap liegen (theoretisch kreuz- und quer verteilt). Je nach deren grösse macht es Sinn die in einem Array auf dem Stack zu haben, dann sind z.B. die lokalen BigInts auf denen man gerade arbeitet im (Stack-)RAM "nah" beieinander und man kann besser von den CPU Cache-Leveln profitieren. Die std::strings haben wegen sowas z.B. die SSO (Small String Optimization). Kurze Strings liegen auf dem Stack, erst wenn sie länger werden wird der Heap genommen. Eine richtig gute BigInt-Library würde das wahrscheinlich ähnlich machen.

Finnegan

P.S.: Noch eine einfachere Möglichkeit, mit der man noch ein bisschen rausquetschen könnte: Whole Program Optimization und Link Time Code Generation. Wenn die cpp-Dateien der BigInt-Lib Teil des Projekts sind kann er dann auch aus diesen Funktionen inlinen oder die Nutzung der CPU-Register für das ganze Programm optimieren. Es gibt Fälle wo sowas durchaus nochmal um die 10% bringen kann, wenn mn Glück hat

Erhard Henkes

Bei den großen Zahlen wie 10^20 oder sogar 10^240 sind die gefundenen großen Primzahlen doch viel zu wertvoll, um sie gleich wieder wegzuwerfen.

Ja, das ist bei großen Zahlen völlig richtig.

Was ist das Ziel des ganzen hier? Willst du so weit kommen, wie du nur kannst? Dann würde ich jetzt eher aufhören, zu mikrooptimieren, sondern mir die BigInt-Implementierung ansehen.

Verbesserte Abläufe, mulmod und BigUnsigned sind die Brennpunkte. So etwas in der Art habe ich bisher noch nicht unternommen. Ich finde es interessant zu sehen, was man heute auf dem eigenen PC untersuchen kann und welche Optimierungsstrategien dabei wichtig sind.

Eine Widerlegung durch Finden einer Anti-Goldbach-Zahl erwarte ich nicht mehr, aber genau das versucht das Programm.

Erhard Henkes

Ich habe den Vorschlag von DJohn, hohe Primzahlen zu speichern und vor dem Miller-Rabin per Lookup zu checken, mal auf die Schnelle eingefügt:
Den Erfolg kann man durch Aktivierung von // cout << " LH "; in isPrime(...) überprüfen.

#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 <vector>
#include <unordered_set>
#include <fstream>
#include <thread>
#include <mutex>
#include "BigIntegerLibrary.hh" // https://mattmccutchen.net/bigint/

using namespace std;

///#define MILLER_RABIN_LIAR_CHECK
#define FUNCTION_TIME_ANALYSIS

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

BigUnsigned hiPrimes[100000];   // calculated primes (high values)
uint64_t countHiNumbersToBeTested = 0;
uint64_t countSuccessfulDivisionPrecheck = 0;
uint64_t countSavedMR = 0;
uint64_t countMR = 0;

uint64_t frequency;

#ifdef FUNCTION_TIME_ANALYSIS
	uint64_t startCount=0, lastCount=0;
	uint64_t sumTime_Function = 0;
#endif

void setConsole()
{
	_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)
{
	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 mutex0;
mutex mutex1;
mutex mutex2;
mutex mutex3;
mutex mutex4;

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

BigUnsigned mulmod(BigUnsigned a, BigUnsigned b, const BigUnsigned& mod)
{
#ifdef FUNCTION_TIME_ANALYSIS
	QueryPerformanceCounter((LARGE_INTEGER*)&startCount);
#endif

	BigUnsigned x = 0;
	BigUnsigned& y = a;

	while (b > 0)
	{
		if (b.getBit(0))
		{
			x += y;
			x %= mod;
		}
		y <<= 1;
		y %= mod;
		b >>= 1;
	}

#ifdef FUNCTION_TIME_ANALYSIS
	QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
	sumTime_Function += (lastCount - startCount);
#endif

	return x;
}

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.getBit(0))
		{
			z1 >>= 1;
			a1 = mulmod(a1, a1, m);
		}
		--z1;
		x = mulmod(x, a1, m);
	}
	return x;
}

bool IsProbablyPrime_MillerRabinOptimized(BigUnsigned n)
{
    mutex4.lock();
    countMR++;
    mutex4.unlock();

	/*
	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) - low primes
	if (number <= primeLimit)
		return primes[number.toUnsignedLong()];

	// lookup from bitset (in the RAM) - high primes
	if ((endNumber - number) <= 100000)
	{
		uint32_t element = (endNumber - number).toUnsignedLong();
		if (hiPrimes[element] == 1)
		{
			// cout << " LH Prime: endN - " << element << " ";
			mutex3.lock();
			countSavedMR++;
			mutex3.unlock();
			return true;
		}
		if (hiPrimes[element] == 0)
		{
			// cout << " LH No-Prime: endN - " << element << " ";
			mutex3.lock();
			countSavedMR++;
			mutex3.unlock();
			return false;
		}
	}

	// division test with first prime numbers
	const uint32_t DIVISOR_PRECHECK_MAX = 1999993;
	for (uint32_t divisor = 3; divisor <= DIVISOR_PRECHECK_MAX; ++divisor)
	{
		if (primes[divisor])
		{
			if (number % divisor == 0)
			{
				if ((endNumber - number) <= 100000)
				{
					uint32_t element = (endNumber - number).toUnsignedLong();
					hiPrimes[element] = 0;
				}
				mutex2.lock();
				countSuccessfulDivisionPrecheck++;
				mutex2.unlock();
				return false;
			}
		}
	}

	// Miller-Rabin, not perfect, but fast
	if (!IsProbablyPrime_MillerRabinOptimized(number))
	{
		if ((endNumber - number) <= 100000)
		{
			uint32_t element = (endNumber - number).toUnsignedLong();
			hiPrimes[element] = 0;
		}
		return false;
	}
	else
	{
		if (lastLiar == 0)
		{
			return true;
		}

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

		if (number > bigLastLiar)
		{
			if ((endNumber - number) <= 100000)
			{
				uint32_t element = (endNumber - number).toUnsignedLong();
				hiPrimes[element] = 1;
			}
			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 >> 1);
	BigUnsigned a = 0;
	for (a = 3; a <= grenze; a += 2)
	{
	    bool retValLo = IsPrime(primes, a, primeLimit, setLiar);
	    if (retValLo == false)
	    {
	        continue;
	    }

        mutex1.lock();
	    countHiNumbersToBeTested++;
	    mutex1.unlock();

	    bool retValHi = IsPrime(primes, i - a, primeLimit, setLiar);

		if (retValLo && retValHi)
		{
			*retVal = true;
			break;
		}
	}

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

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

	QueryPerformanceFrequency((LARGE_INTEGER*)&frequency);
	cout << "cpu frequency: " << frequency << " Hz" << endl << endl;

	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); // calculated primes (low values)

    // calculated primes (high values)
	for (int element = 0; element < 100000; ++element)
	{
		if ((element & 1))
			hiPrimes[element] = 2; // not determined
		else
			hiPrimes[element] = 0; // no prime
	}
	// endNumber (index 0) to endnumber-100000 (index 100000)
	// 0: no prime
	// 1: (probably) prime
	// 2: not determined
	// 3: just being calculated

	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

	uint64_t zeit1, zeit2, zeit2_old, delta;
	QueryPerformanceCounter((LARGE_INTEGER*)&zeit1);
	zeit2 = zeit1;
	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, std::ref(primes), i+2*x, BigUnsigned((unsigned long)primeLimit), std::ref(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;
		QueryPerformanceCounter((LARGE_INTEGER*)&zeit2);
		delta = (zeit2 - zeit2_old);
		sumTime += delta;
		cout << "time: " << 1000 * delta / frequency << " ms" << endl << endl;
	}

	uint64_t totalTimeMS = 1000 * sumTime / frequency;
	cout << "total time: " << totalTimeMS << " ms" << endl << endl;

#ifdef FUNCTION_TIME_ANALYSIS
	uint64_t totalTimeAtFunctionMS = 1000 * sumTime_Function / frequency;
	cout << "total time at mulmod: " << totalTimeAtFunctionMS << " ms (" << 100 * totalTimeAtFunctionMS / totalTimeMS << " %)" << endl << endl;
#endif

    cout << "High Numbers to be tested: " << countHiNumbersToBeTested << endl;
    cout << setprecision(3);
    cout << "Number of successful Division Prechecks: " << countSuccessfulDivisionPrecheck << " (" << 100 * (double)countSuccessfulDivisionPrecheck/(double)countHiNumbersToBeTested << " %)" << endl;
    cout << "Number of MR tests: " << countMR << " (" << 100 * (double)countMR/(double)countHiNumbersToBeTested << " %)" << endl;
    cout << "Number of saved MR tests by high prime lookup: " << countSavedMR << " (" << 100 * (double)countSavedMR/(double)countHiNumbersToBeTested << " %)" << endl;

	wait();
}

`from 10^60 to 10^60 + 100:

High Numbers to be tested: 1954

Number of successful Division Prechecks: 547 (28

Number of MR tests: 104 (5.32

Number of saved MR tests by high prime lookup: 1303 (66.7

`

`from 2 * 10^60 to 2 * 10^60 + 1000:

High Numbers to be tested: 46574

Number of successful Division Prechecks: 1951 (4.19

Number of MR tests: 752 (1.61

Number of saved MR tests by high prime lookup: 43871 (94.2

`

Die Idee des "high prime lookup" spart gerade beim Testen vieler Zahlen in einem hohen Zahlenbereich eine Vielzahl an Miller-Rabin-Tests.

Thanks to DJohn.

Erhard Henkes

Für weitere Optimierungen nehme ich den Goldbach-Test (Geradzahl, einfaches Primzahlenpaar) bei 2*10^120 bis 2*10^120 + 1000 als Basis.
`

total time: 426900 ms

total time at mulmod: 248637 ms (58

High Numbers to be tested: 98055

Number of successful Division Prechecks: 2867 (2.92

Number of MR tests: 788 (0.804

Number of saved MR tests by high prime lookup: 94400 (96.3 %)`
komplett: http://pastebin.com/ziNH4sna

Zum eventuellen Herumschrauben an der freien BigInteger Library habe ich diese auf eine einfache h/cpp-Kombination konzentriert und auf das Wesentliche reduziert: http://henkessoft.de/Sonstiges/MyBigUnsigned.rar
Ich bin allerdings nicht sicher, ob das Sinn macht.