Goldbach Vermutung

[volkard]

Hab da noch was in der Krimskramskiste gefunden.
Vielleicht besser, wenn es Bengo vorher mal testet, es ist aus einem Experimetierverzeichnis und Bengo scheint erdbebenfester zu sein. Wenns klappt, haste den Millimillerrabin instant.

UInt64 mulModUnguarded(UInt64 a,UInt64 b,UInt64 n){
	assert(a<n);
	assert(b<n);
//Beware: if (a*b)/n > 2^32, there will be an FP exception
//cout<<a<<' '<<b<<' '<<n<<'\n';
    UInt64 d, dummy; /* d will get a*b mod n */
    asm ("mulq %3\n\t" /* mul a*b -> rdx:rax */
         "divq %4\n\t" /* (a*b)/n -> quot in rax remainder in rdx */
         :"=a"(dummy), "=&d"(d) /* output */
         :"a"(a), "rm"(b), "rm"(n) /* input */
         :"cc" /* mulq and divq can set conditions */
        );
    return d;
}

[Erhard Henkes]

Momentan führt das aktuelle mod_pow zu Problemen beim Miller-Rabin. Von daher ist Volkards Vorschlag interessant. Dort spielt die Musik.

1. Daten bis 100 Mio. reichen im bitset.
2. Das powmod von Arcoth führt zu besseren Zahlen
3. Ich poste meinen gesamten Code, der zur Zeit gut läuft:

#include <iostream>
#include <iostream>
#include <iomanip>
#include <cstdint>
#include <cstdlib>
#include <algorithm>
#include <cmath> 
#include <ctime>
#include <vector>

using namespace std;

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

//calculates (a * b) % c taking into account that a * b might overflow
uint64_t mulmod(uint64_t a, uint64_t b, uint64_t mod)
{
	uint64_t x = 0;
	uint64_t y = a % mod;

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

uint64_t powmod(uint64_t x, uint64_t exp, uint64_t mod)
{
	uint64_t res = 1;

	for (; exp; exp >>= 1)
	{
		if (exp & 1)
			res = mulmod(res, x, mod);

		x = mulmod(x, x, mod);
	}

	return res;
}

// lame duck!
/* 
uint64_t mod_pow(uint64_t base, uint64_t expo, uint64_t mod) 
{
	uint64_t x = 1;
	uint64_t aktpot = base;
	while (expo > 0) 
	{
		if (expo % 2 == 1) 
		{
			x *= aktpot;
		}
		aktpot = (aktpot * aktpot) % mod;
		expo /= 2;
	}
	return x % mod;
}
*/

bool IsPrimeMillerRabinOptimized(uint64_t number, uint64_t base)
{
	uint64_t d = number - 1;
	int counter = 0;
	while (d % 2 == 0)
	{
		d /= 2;
		++counter;
	}

	uint64_t temp = powmod(base, d, number);

	if (temp == 1 || temp == number - 1) //Hier
	{
		return true;
	}

	for (int i = 0; i < counter; ++i)
	{
		temp = (temp * temp) % number;

		if (temp == number - 1)
		{
			return true;
		}
	}
	return false;
}

bool IsPrimeDivisionTest(uint64_t number)
{
	if (number<2)    return false;
	if (number == 2)   return true;
	if (number % 2 == 0) return false;

	uint64_t maxI = uint64_t((sqrt((double)number)+1));

	for (uint64_t i=3; i<=maxI; i+=2)
	{
		if (number%i == 0) return false;
	}
	return true;
}

bool IsPrime(vector<bool>& primes, uint64_t number, uint64_t primeLimit)
{
	if (number <= primeLimit)
		return primes[number]; // lookup from bitset (in the RAM)
	if (!IsPrimeMillerRabinOptimized(number, 2))
		return false;
	return IsPrimeDivisionTest(number);
}

bool PrimePairFound(vector<bool>& primes, uint64_t i, const uint64_t primeLimit)
{
	bool retVal = false;
	uint64_t grenze = i/2;
	uint64_t a = 0;
	for (a=3; a<=grenze; a+=2)
	{
		if (IsPrime(primes, a, primeLimit) && IsPrime(primes, i-a, primeLimit))
		{
			retVal = true;
			break;
		}
	}
	cout << i << "  " << a << " " << i-a << " ";
	return retVal;
}

int main()
{
	uint64_t startNumber, endNumber;
	cout << "Start number to be executed: ";
	cin >> startNumber;
	cout << "Last number to be executed:  ";
	cin >> endNumber;

	uint64_t const primeLimitByRAM = 100000000; //200000000000;
	const uint64_t primeLimit = min(primeLimitByRAM, endNumber);

	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;

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

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

	for (uint64_t i = startNumber; i <= endNumber; i += 2) 
	{
		if (!PrimePairFound(primes, i, primeLimit))
		{
			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();
}
}

Damit schnackelt es wieder. Doppelcheck mit Divisionstest kann man natürlich weglassen, wenn man Goldbach widerlegen will, was man wohl vergessen kann.

Hier mal mit 2 Trillionen ohne double check mit Divisionstest:
`Start number to be executed: 2000000000000000000

Last number to be executed: 2000000000000000010

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

In 0 s we found 5761455 prime numbers (2 to 100000000 ).

Now Goldbach first prime number pair will be detected

2000000000000000000 601 1999999999999999399 time: 10 ms

2000000000000000002 131 1999999999999999871 time: 4 ms

2000000000000000004 593 1999999999999999411 time: 9 ms

2000000000000000006 397 1999999999999999609 time: 8 ms

2000000000000000008 137 1999999999999999871 time: 4 ms

2000000000000000010 139 1999999999999999871 time: 4 ms`

Arcoths powmod ist klasse. bengo hatte zwischendrin aber auch einen guten Code gepostet, dann aber wieder kaputt gemacht. Deshalb verwende ich Arcoths code.

[Columbo]

@volkard: mulmod sieht so aus:

u64 mulmod(u64 a, u64 b, u64 m)
{
	u64 r;
	asm
	( "mulq %2;"
	  "divq %3;"
	: "=&d" (r), "+%a" (a)
	: "rm"  (b), "rm"  (m)
	: "cc"
	);
	return r;
}

…und funktioniert wunderbar bei allen Werten.

[Erhard Henkes]

@Arcoth: nein, dein Code ist klasse. Liefert ideal den Partner zu niedrigen Primzahlen. Dadurch geht das rasch. Da ist echt viel Musik zum Optimieren drinnen.

[volkard]

Arcoth schrieb:

@volkard: mulmod sieht so aus:

u64 mulmod(u64 a, u64 b, u64 m)
{
	u64 r;
	asm
	( "mulq %2;"
	  "divq %3;"
	: "=&d" (r), "+%a" (a)
	: "rm"  (b), "rm"  (m)
	: "cc"
	);
	return r;
}

…und funktioniert wunderbar bei allen Werten.

thx.

[volkard]

Arcoth schrieb:

…und funktioniert wunderbar bei allen Werten.

Nee, dafür hab ich doch den Kommentar reingeschrieben.

int main()
{
	UInt64 a=1;
	UInt64 b=1;
	for(;;){
		a=a*a+1;
		b=b*b+1;
		b=b*b+1;
		cout << a<<' '<<b<<' '<<10 << endl;
		cout << mulmod(a,b,10) << endl;
    }
    return 0;
}

2 5 10
0
5 677 10
5
26 210066388901 10
6
677 6521269100592521125 10

Process returned -1 (0xFFFFFFFF)   execution time : 0.003 s
Press ENTER to continue.

Klappt nur, wenn man vorher mathematische Überlegungen macht, warum's klappt. Deswegen bei mir der lange nameUnguarded. Den dürfen schlaue Funktionen wie powmod benutzen, weil sie garantieren, dass alles ok ist.

Um genau zu sein isses dann auchnur die powmodUnguarded, die wiederum aus dem Godbach-Such-Programm weiß, daß sie nicht mit ungültigen Parametern aufgerufen wird.

[Erhard Henkes]

Bei mir spielt der compiler nicht mit bei x64 (nur bei x32). Intrinsics kann ich noch nicht. Außerdem kommt nun die 2^64 Barriere.

[Bengo]

Erhard Henkes schrieb:

Bengo: bitte prüfe nochmal den Miller-Code. Es geht gar nix mehr. Offenbar hast Du irgendwas gekillt, habe leider den funktionierenden Code überschrieben.

hab noch mal geguckt, der Code funktionert bei mir einwandfrei. Hab ihn auch noch mal aus dem Forum kopiert. Die mod_pow würde ich aber von Arcoth nehmen, die mulmod Funktion kannte ich noch nicht, macht die Sache aber glaube etwas schneller. (Kannst ja mal vergleichen welche schneller ist)

[Columbo]

@volkard: Das ist merkwürdig, denn ich habe dieses mulmod so schon seit Ewigkeiten verwendet, und es funktionierte Prima. Aber tatsächlich!

[Erhard Henkes]

Bengo: zur Kärung, weil hier momentan einiges parallel läuft, ich verwende momentan folgendes für Ziel x64, was bestens läuft:

//calculates (a * b) % c taking into account that a * b might overflow
uint64_t mulmod(uint64_t a, uint64_t b, uint64_t mod)
{
	uint64_t x = 0;
	uint64_t y = a % mod;

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

uint64_t powmod(uint64_t x, uint64_t exp, uint64_t mod)
{
	uint64_t res = 1;

	for (; exp; exp >>= 1)
	{
		if (exp & 1)
			res = mulmod(res, x, mod);

		x = mulmod(x, x, mod);
	}

	return res;
}

Das klappt hervorragend mit Bengos Miller-Rabin.

Volkards asm mulmod bekomme ich mit MS VS unter x64 nicht zum Laufen, nur mit x86, und da ist der asm code langsamer als mit der C-Variante mit x64.

Zum Vergleich (s.o.):
`Now Goldbach first prime number pair will be detected

2000000000000000000 601 1999999999999999399 time: 41 ms

2000000000000000002 131 1999999999999999871 time: 13 ms

2000000000000000004 593 1999999999999999411 time: 41 ms

2000000000000000006 397 1999999999999999609 time: 30 ms

2000000000000000008 137 1999999999999999871 time: 13 ms

2000000000000000010 139 1999999999999999871 time: 14 ms`

[Bengo]

Erhard Henkes schrieb:

Der Weg bisher ist nicht schlecht, habe einiges gelernt. Nur diesen Miller-Rabin-Algo habe ich noch nicht verstanden, sonst würde ich den Fehler, der jetzt hinein geraten ist, selbst finden. So brauche ich wieder Bengo.

Also ich stütze mich dabei auf den Wikipedia artikel:

https://de.wikipedia.org/wiki/Miller-Rabin-Test

Es geht darum zu zeigen, dass a^{n-1} \equiv 1 \mod n gilt.
Um die Berechnungen etwas abzukürzen, teilt man aus n-1 erst die größte Zweierpotenz raus, die zahl die wir erhalten nennen wir d. Wenn a^d \equiv 1 \mod n, dann ist unsere Bedinung erfüllt. 1^irgentwas, ist schießlich immer wieder 1.

Nun quadriert man das ganze so lange bis man -1 erreicht. Weil dann wird es am Ende auch wieder eine 1. (eine 1 kann hier nicht vorkommen, da 1 nur aus -1 erreicht werden kann, und das erste Glied war ja keine 1)

Ich habe gerade die Befürchtung, dass dieser ganze Milleralgo bei der mod_pow Funktion gar keinen richtigen Vorteil mehr bringt.

[volkard]

Erhard Henkes schrieb:

Bei mir spielt der compiler nicht mit bei x64 (nur bei x32).

Welcher BS, welcher Compiler, welcher Prozessor?

Erhard Henkes schrieb:

Außerdem kommt nun die 2^64 Barriere.

Jo, die kommt gnadenlos. Selbst wenn jeder Primzahlentest nur einen Takt brauchen würde. Willst ja RIESIGE Zahlenbereiche untersuchen.

Das macht aber auch voll Spaß. Nachdem die Primzahlen sauflott sind, diesbezüglich würde ich am Ball bleiben.

Danach machste wieder pro Woche so viel Beschleunigung durch verbessertes Programm, dass immer es nichtig war, auch nur anzufangen, Ergebnisse zu speichern. Fast endloser Spaß!

[volkard]

Arcoth schrieb:

@volkard: Das ist merkwürdig, denn ich habe dieses mulmod so schon seit Ewigkeiten verwendet, und es funktionierte Prima.

Hab ich auch ewig verwendet. Bis ich mal angefangen hatte, Unit-Tests zu bauen.
Es ist ja auch perfekt hierfür.

[Erhard Henkes]

Danke für die Erklärungen. Echt tolle Unterstützung hier. Dickes Lob an euch!

Nochmal zur praktischen Klarstellung:

x86 mit dem C-code für mulmod:
`

Now Goldbach first prime number pair will be detected

2000000000000000000 601 1999999999999999399 time: 10 ms

2000000000000000002 131 1999999999999999871 time: 4 ms

2000000000000000004 593 1999999999999999411 time: 9 ms

2000000000000000006 397 1999999999999999609 time: 8 ms

2000000000000000008 137 1999999999999999871 time: 4 ms

2000000000000000010 139 1999999999999999871 time: 4 ms`

x32 mit dem asm code für mulmod:
`

Now Goldbach first prime number pair will be detected

2000000000000000000 601 1999999999999999399 time: 41 ms

2000000000000000002 131 1999999999999999871 time: 13 ms

2000000000000000004 593 1999999999999999411 time: 41 ms

2000000000000000006 397 1999999999999999609 time: 30 ms

2000000000000000008 137 1999999999999999871 time: 13 ms

2000000000000000010 139 1999999999999999871 time: 14 ms`

Die Zielplattform x86 bringt hier also wirklich Vorteile.

Wenn man das in x64 asm umwandeln könnte, wäre das sicher interessant.

Mich treibt inzwischen aber eher die Frage um, welche Lib man nimmt, um in den Bereich 10^19 usw. vorzudringen. Da wird asm code für u64 sicher stören. C ist da wohl portabler.

[volkard]

Bengo schrieb:

Ich habe gerade die Befürchtung, dass dieser ganze Milleralgo bei der mod_pow Funktion gar keinen richtigen Vorteil mehr bringt.

Doch.
Millimiller (also strake Pseudoprimzahlen) ist schneller und viel stärker als nur Pseudoprimzahlen wie in https://www.c-plusplus.net/forum/333757-19

Gegen 2^64 kackt TrialDivision natürlich total ab. Aber auch Eratostenes kackt wonderbar ab. Bei Miller muss man immer schauen, was gemeint ist. Deswegen hab ich Millimillerrabin als Wort erfunden, damit klar ist, daß man hier nur einen SPRP-Test macht, statt bis unten hin zu millern und dennoch nicht ganz sicher ist. Nur der eine SPRP-Test macht nicht sicher, dazu doch die Lügner-Liste.

[Erhard Henkes]

Welches BS, welcher Compiler, welcher Prozessor?

Windows 7

MS VS 2015 Community

Intel Core i7 3930K (Sandy Bridge-E) mit 6 echten und 12 virtuellen Kernen (cpu noch bei 8%)

3,2 GHz (könnte man sicher tunen)
32 GB RAM (bitset bis 200 Mrd. passt hier rein, Berechnung 52 min., Laden fehlt noch, z.Z. nicht notwendig)

MMX, SSE, SSE2, SSE3, SSSE3, SSE 4.1 u. 4.2, EM64T, VT-x, AES, AVX (hier habe ich AVX2 eingestellt bei Optimierung)

[volkard]

Erhard Henkes schrieb:

Nochmal zur praktischen Klarstellung:
x86 mit dem C-code für mulmod:
`

Now Goldbach first prime number pair will be detected

2000000000000000000 601 1999999999999999399 time: 10 ms

2000000000000000002 131 1999999999999999871 time: 4 ms

2000000000000000004 593 1999999999999999411 time: 9 ms

2000000000000000006 397 1999999999999999609 time: 8 ms

2000000000000000008 137 1999999999999999871 time: 4 ms

2000000000000000010 139 1999999999999999871 time: 4 ms`

x32 mit dem asm code für mulmod:
`

Now Goldbach first prime number pair will be detected

2000000000000000000 601 1999999999999999399 time: 41 ms

2000000000000000002 131 1999999999999999871 time: 13 ms

2000000000000000004 593 1999999999999999411 time: 41 ms

2000000000000000006 397 1999999999999999609 time: 30 ms

2000000000000000008 137 1999999999999999871 time: 13 ms

2000000000000000010 139 1999999999999999871 time: 14 ms`

Die Zielplattform x86 bringt hier also wirklich Vorteile.

Junge, junge, Du hast gerde für 2000000000000000000 augenblickliche Ergebnisse. Gestern hatten wir noch auf die Primzahlen bis 10000000000 einen kaffetrinkenwarten müssen bevor die Untersuchung anfing.

[Bengo]

Hab mal die jetige MillierRabin (mit basis 2) funktion mit der divisionprime Funktion gleichgesetzt. Es gibt bei mir aber ein paar mehr unstimmigkeiten, als es eigentlich geben sollte:
2027, 2039, 2047 (ok das ist normal), und dann gibt es noch ein paar mehr

[Bengo]

War ein Fehler in meiner implementierung der modpow Funktion.

Mit der C-Variante unter 64bit funktionert es so wie es soll.

übrigens kann man for(int i = 0;...) durch for(int i = 1; ...) ersetzen.
Eigentlich sogar muss.

[Erhard Henkes]

Junge, junge, Du hast gerde für 2000000000000000000 augenblickliche Ergebnisse. Gestern hatten wir noch auf die Primzahlen bis 10000000000 einen kaffetrinkenwarten müssen bevor die Untersuchung anfing.

Ja, begeistert mich auch sehr. Hätte ich nicht erwartet, dass man in dem Bereich noch so gut arbeiten kann, und wir haben noch Multithreading in der Hinterhand und x64 Assembler-Routinen (das sollte man aber vermeiden, solange man die Vorgehensweise und die Algos variiert/optimiert).

for(int i = 1; ...)

Habe ich abgeändert. Vergleich:

`Start number to be executed: 2000000000000000000

Last number to be executed: 2000000000000000010

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

In 1 s we found 5761455 prime numbers (2 to 100000000 ).

Now Goldbach first prime number pair will be detected

2000000000000000000 601 1999999999999999399 time: 9 ms

2000000000000000002 131 1999999999999999871 time: 4 ms

2000000000000000004 593 1999999999999999411 time: 10 ms

2000000000000000006 397 1999999999999999609 time: 7 ms

2000000000000000008 137 1999999999999999871 time: 4 ms

2000000000000000010 139 1999999999999999871 time: 4 ms`

Ich habe getestet, wie weit wir uns mit dem code an 2^64 heran arbeiten können. 10^19 klappt bisher nicht, aber 8 Trillionen (9 geht auch blitzschnell) geht locker:

`Start number to be executed: 8000000000000000000

Last number to be executed: 8000000000000000010

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

In 0 s we found 5761455 prime numbers (2 to 100000000 ).

Now Goldbach first prime number pair will be detected

8000000000000000000 97 7999999999999999903 time: 3 ms

8000000000000000002 41 7999999999999999961 time: 3 ms

8000000000000000004 43 7999999999999999961 time: 3 ms

8000000000000000006 103 7999999999999999903 time: 3 ms

8000000000000000008 47 7999999999999999961 time: 3 ms

8000000000000000010 107 7999999999999999903 time: 4 ms`

Ich lasse das Programm nun von 9 nach 10 Trillionen "durchsausen", um festzustellen, wo es aushakt.