3. Goldbach Vermutung

Goldbach Vermutung


  • 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.

    0
  • F

    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 😃

    0
  • E

    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;
}
    0
  • F

    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 🙄

    0
  • E

    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)

    0
  • E

    ... 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.

    0
  • F

    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 😃

    0
  • E

    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.

    0
  • J

    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.

    0
  • D

    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.

    0
  • F

    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

    0
  • E

    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. 😃 😉

    0
  • E

    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. 👍

    0
  • E

    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.

    0
  • V

    Erhard Henkes schrieb:

    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. 👍

    Auch wenn ich mich wieder unbeliebt mache, schau mal nach "fensterweise" auf der ersten Seite dieses Threads. Ich fürchte, dahin geht es zunächst. Wir werden festestellen, daß bis zu grotesken Größen die beiden äußeren Fester reichen; evenuelle Zweifler kann man dann ja noch extra untersuchen.

    Ich denke gerade darüber nach, statt der äußeren Fenster nur ein einziges inneres Fenster um sqrt(n) zu führen. Dann könnte man bis 2^128 untersuchen obwohl man nur isPrime(UInt64) hat.

    An einem schnellen isPrime(UInt64) arbeite ich seit einer Woche, alse genau die Gegenrichtung von Dir.

    0
  • E

    Mich überrascht bei diesem Thema (noch immer unbewiesene Vermutung des Herren Goldbach) die Vielfalt der Möglichkeiten bei der Herangehensweise und dass man mit einem leistungsstarken PC in wirklich beeindruckende Zahlenbereiche vordringen kann. Ich bin gespannt auf Volkards Neuerungen.

    0
  • B

    Der Satz ist zwar unbewiesen, aber man kann sehr einfach beweisen, dass es verdammt unwahrscheinlich wäre, wenn er nicht wahr wäre 😃

    0
  • ?

    Erhard Henkes schrieb:

    [...] und dass man mit einem leistungsstarken PC in wirklich beeindruckende Zahlenbereiche vordringen kann. Ich bin gespannt auf Volkards Neuerungen.

    Wobei man bei der ganzen Sache nicht vergessen sollte, das man zwar im Bereich 10^120 einzelne Zahlen in ein paar Millisekunden testen kann aber nicht alle Zahlen bis dahin. Selbst wenn du die gesamte Rechenkapazität der Erde nutzen würdest, kriegst du bis zum Weltuntergang nicht alle Zahlen getestet. Nur so als Vergleich: Das Universum ist ca. 4*10^17 Sekunden alt.

    0
  • E

    sebi707: Ja, das ist richtig. Egal, wie toll ein Algorithmus oder Rechner ist, einfach einige hundert Zehnerpotenzen weiter, und schon ist Schluss. Alle Zahlen bis 2^64 zu checken erfordert umfangreiches verteiltes Rechnen. Bei Collatz habe ich da mal vor ein paar Jahren Teil genommen. Da gabs 20000er Pakete. http://www.ericr.nl/wondrous/

    Hier ist ein interessanter Artikel zum Berechnen der Goldbach Vermutung (schon 15 Jahre alt): http://www.ams.org/journals/mcom/2001-70-236/S0025-5718-00-01290-4/S0025-5718-00-01290-4.pdf

    0
  • E

    Ich bin mit dem aktuellen Programm wieder etwas weiter ins "Zahlen-Universum" vorgestoßen und habe den Bereich von 10^360 bis 10^360 + 100 durchgerechnet:

    `total time: 9302703 ms

    total time at mulmod: 8612644 ms (92 😵

    High Numbers to be tested: 28593

    Number of successful Division Prechecks: 8801 (30.8 😵

    Number of MR tests: 758 (2.65 😵

    Number of saved MR tests by high prime lookup: 19034 (66.6 %)`

    Alle Ergebnisse: http://pastebin.com/Ftgq6eMt

    Die höchste Spanne nach unten zur nächsten Primzahl ist dort 34259.

    0
Anmelden zum Antworten