Goldbach Vermutung

volkard

Erhard Henkes schrieb:

Sorry, ihr habt es geschafft. Bin draußen. Es sei denn, jemand schreibt diese beiden Funktionen in C++ vollständig hin.

Sorry. Das hatten wir Dir aus Versehen verheimlicht.
Alle Konsorten des MilliMillers leben von der unglaublich schnellen mod_pow von Bengo.
2^7647768764054048747840%67283 brauchen da nicht 7647768764054048747840 teure Divisionen(oder Mod-Rechnungen), sondern nur log₂(7647768764054048747840)=73 davon. Deswegen geht da echt die Post ab, wenns Richtung 64-Bit geht, er wird ums Verrecken nicht sprürbar langsamer, wärend alle anderen Ecken des Untersuchungsprogramms sich quadratische Laufzeit gönnen und man an jenen Ecken nicht sieht, wie man gen 64Bit ernsthaft kommen kann.

Erhard Henkes

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.

volkard

Erhard Henkes schrieb:

Interessant ist, wenn man bei den von Miller-Rabin gefunden Primes noch einen double check macht mit dem Divisions-Test:

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

Genau.
Wobei ich es eher sorum mag:

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

//sobald ein Ergebnis feststeht, sofort(!) mit dem neuen Wissen rausreturnen.

Und jetzt erschlägt Dich einfach, daß es so unglaublucg viele Primzahlen doch noch da oben gibt. Da schneckt der IsPrimeDivisionTest total.

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

	return IsPrimeDivisionTest(number); // double check!
}

Und dann noch

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

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.

Daten bis 100 Mio. reichen im bitset.
Das powmod von Arcoth führt zu besseren Zahlen
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.