Goldbach Vermutung

Bengo

Denke dem MillerRabin Test sollte man sowieso vorher einen Divisionstest mit den Zahlen 2,3,5,7,11 voranstellen, das geht schnell und kann schon einen Großteil der zusammengesetzten Zahlen vorab erkennen. Aber an einer BigInterger Lib wirst du nicht vorbeikommen, wenn du die Grenze 2^32 bzw. 2^64 überschreiten willst. Hab selbst noch keine benutzt, kann dir daher da auch keine gute Empfehlung machen.

Bengo

Hab mal die mzp_class der gmp genutzt:

#include <iostream>
#include <cmath>
#include <gmpxx.h>

using namespace std;

mpz_class mulmod(mpz_class a, mpz_class b, mpz_class mod)
{
    mpz_class x = 0;
    mpz_class y = a % mod;

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

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

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

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

    return res;
}

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

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

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

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

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

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

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

int main(int argc, char** argv) {

    for (int i = 2; i < 10000; ++i) {
        cout << i<<":"<< (int)(IsPrimeMillerRabinOptimized(i, 2)==IsPrimeDivisionTest(i))<< endl;
    }

}

Bei der Zahl 100000000000000000039, gibt der Miller zur Basis 2, aus, dass es eine Primzahl sein könnte. Die einfache Division kommt auf meinem Computer nicht mehr zum Ende. Und die zahl ist tatsächlich prim.

Auch die zahl 1000000000000000000000000000, macht keine Probleme.

Die Zahl 1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000019, ist prim zur Basis 2. Wird auch von http://primzahlen.zeta24.com/de/online_primzahltest.php bestätigt. Der Test dauert etwa eine halbe Sekunde.

Erhard Henkes

Bengo: Ja, Du hast Recht. Wir müssen eine Stufe höher klettern. Probiere bitte https://mattmccutchen.net/bigint/ C++ Big Integer Library als Vergleich dazu aus, welches schneller ist (GMP soll schneller sein. Die Frage ist um wieviel). Schaffe es momentan nicht gmpxx/gmp unter Windows zum Laufen zu bekommen.

Recht einfach für Dich:

//#include <gmpxx.h>
#include "BigIntegerLibrary.hh"

#define mpz_class BigUnsigned

Auf diese Weise findet man alle Lügner (falsche Primes nach Miller-Rabin) bis 10^6:

int main(int argc, char** argv)
{

     for (int i = 3; i < 1000000; ++i)
     {
         if (((int)IsPrimeMillerRabinOptimized(i, 2) == 1 && (int)IsPrimeDivisionTest(i) == 0))
             cout << i << " Luegner" << endl; // Rabin-Miller Lügner
         if (((int)IsPrimeMillerRabinOptimized(i, 2) == 0 && (int)IsPrimeDivisionTest(i) == 1))
             cout << i << " Ueberseher" << endl; // Rabin-Miller Primzahlen-Überseher
     }     

     wait();
 }

`15841 Luegner

29341 Luegner

42799 Luegner

49141 Luegner

52633 Luegner

65281 Luegner

74665 Luegner

80581 Luegner

85489 Luegner

88357 Luegner

90751 Luegner

104653 Luegner

130561 Luegner

196093 Luegner

220729 Luegner

233017 Luegner

252601 Luegner

253241 Luegner

256999 Luegner

271951 Luegner

280601 Luegner

314821 Luegner

357761 Luegner

390937 Luegner

458989 Luegner

476971 Luegner

486737 Luegner

489997 Luegner

514447 Luegner

580337 Luegner

635401 Luegner

647089 Luegner

741751 Luegner

800605 Luegner

818201 Luegner

838861 Luegner

873181 Luegner

877099 Luegner

916327 Luegner

976873 Luegner

983401 Luegner`

Bengo

Erhard Henkes schrieb:

Schaffe es momentan nicht gmpxx/gmp unter Windows zum Laufen zu bekommen.

Da gibt es eine ganz einfache Lösung: Linux benutzen! Ich brauchte exakt nur sudo apt-get install libgmp3-dev eingeben (Ubuntu), Enter drücken und im Code includieren, und die gmp Library in der IDE hinzufügen.

Eigentlich sollte es jetzt keine Fälle mehr geben, in denen zu unrecht false zurückgegeben wird. Also ich hab keine mehr gefunden.

Erhard Henkes schrieb:

Schaffe es momentan nicht gmpxx/gmp unter Windows zum Laufen zu bekommen.

GMP ist mal wieder so ein Projekt was die Existenz von Windows ignoriert. Unter Windows fährt man mit MPIR (http://mpir.org/index.html#about) besser. Das Interface ist sogar kompatibel zu GMP.

volkard

Erhard Henkes schrieb:

Ich weiß nicht genau, wo es anfängt, aber es ist unterhalb 2^64 / 2.

Ja, ab dort irgendwo hab ich lauter Überseher. Und Überseher darf es gar nicht geben.
Habs so gefixt.

#include <iostream>
#include <cassert>
#include <cstdint>
using namespace std;

typedef unsigned int UInt32;
static_assert(sizeof(UInt32)==4,"wrong size");

typedef unsigned long long UInt64;
static_assert(sizeof(UInt64)==8,"wrong size");

typedef UInt64 u64;

UInt64 mulmod(UInt64 a, UInt64 b, UInt64 m)
{
	UInt64 r;
	asm
	("mulq %2;"
	 "divq %3;"
	 : "=&d"(r), "+%a"(a)
	 : "rm"(b), "rm"(m)
	 : "cc"
	);
	return r;
}
/*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;
}*/
UInt64 powmod(UInt64 base,UInt64 exp,UInt64 modul){
    //rechnet 'base hoch exp mod modul'
    UInt64 a1=base,z1=exp,x=1,m=modul;
    while(z1!=0){
        while((z1%2)==0){
            z1/=2;
            a1=mulmod(a1,a1,m);
        }
        --z1;
        x=mulmod(x,a1,m);
    }
    return x;
}
bool isSPRP(UInt64 n,UInt64 a)
{
	if(a%n==0) return true;
  UInt64 d=n-1;
  UInt64 ad;
  UInt64 s=0;

  // break down n-1 into d*(2^s). Linux ffs() can be better
  while (0==(d&1)) {
    ++s; d>>=1;
  }

  ad=(powmod(a%n, d, n)); // (a^d) mod n

  if (ad==1) return 1; // 1 == a^d mod n
  if (s>0 && ad==(n-1)) return 1; // -1 == a^d mod n (tests (a^d)^(2^0))
  for (UInt64 r=1; r<s; ++r) {
    // ad is currently ad^(2^(r-1)) so we square it to get ad^(2^r));
    ad=mulmod(ad,ad,n);
    if (ad==(n-1)) return 1; //tests (a^d)^(2^(r+1)))
  }
  return 0; // false, composite
}
/*bool IsPrimeMillerRabinOptimized(uint64_t number, uint64_t base){
	return isSPRP(number,base);
}*/

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;
	for(uint64_t i=3; i*i<=number; i+=2)
	{
		if(number%i == 0) return false;
	}
	return true;
}

int main()
{
	for(uint64_t i = uint64_t(1)<<32; ; ++i)
	{
		if(IsPrimeMillerRabinOptimized(i, 2) && !IsPrimeDivisionTest(i))
			cout << i << " Luegner" << endl; // Rabin-Miller Lügner
		if(!IsPrimeMillerRabinOptimized(i, 2) && IsPrimeDivisionTest(i))
			cout << i << " Ueberseher" << endl; // Rabin-Miller Primzahlen-Überseher
	}
	wait();
}

Erhard Henkes

Also der reihe nach: Die Tests bis 2^32 mit uint64_t ergeben allesamt keine falschen Rückgaben von false. Oberhalb von 2^32 tauchen falsche false auf. Daher benötigt man schnelle BigInteger. gmp soll wohl das schnelleste sein. Ich versuche es mal einzubauen.

Ah, neuer Post von Volkard.

Bengo

Erhard Henkes schrieb:

Also der reihe nach: Die Tests bis 2^32 mit uint64_t ergeben allesamt keine falschen Rückgaben von false. Oberhalb von 2^32 tauchen falsche false auf. Daher benötigt man schnelle BigInteger. gmp soll wohl das schnelleste sein. Ich versuche es mal einzubauen.

genau, da 2^32 im Quadrat 2^64 ergibt. Wenn es noch größer wird, dann gibt es Fehler durch einen overflow. Mit einer BigInteger Library kann ein overflow nicht passieren und auch die Werte oberhalb von 2^32 sind korrekt.

OT aber an denjenigen, der diesen Code ursprünglich geschrieben hat:

#include <cstdint>
// [...] 
typedef unsigned int UInt32;
static_assert(sizeof(UInt32)==4,"wrong size");

typedef unsigned long long UInt64;
static_assert(sizeof(UInt64)==8,"wrong size");

Warum läd man den cstdint Header, hat einen Compiler der offensichtlich C++11 unterstützt, und nutzt dennoch nicht die uint32_t und uint64_t Typen?

volkard

sebi707 schrieb:

OT aber an denjenigen, der diesen Code ursprünglich geschrieben hat:
#include <cstdint>
// [...] 
typedef unsigned int UInt32;
static_assert(sizeof(UInt32)==4,"wrong size");
 
typedef unsigned long long UInt64;
static_assert(sizeof(UInt64)==8,"wrong size");
Warum läd man den cstdint Header, hat einen Compiler der offensichtlich C++11 unterstützt, und nutzt dennoch nicht die uint32_t und uint64_t Typen?

*meld*
Ursprünglich war das Teil einer lib, die nur ganz wenige Standard-Headers oder gar keine benutzt.
Jemand hier benutzte u64, also schnell "typedef UInt64 u64;" nachgezogen.
Die meisten benutzen uint64_t, also das include nachgezogen.

Jetzt kann ich nach Belieben von beiden anderen CopyPasten und es geht.

Nach dem include hätte ich meine typedef umbauen sollen, stimmt, so isses gar häßlich.

Erhard Henkes

Ich brauche Hilfe, da MS VS den asm code nicht annimmt bei x64.

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

Könnt ihr das in C übersetzen? Ich schaffe das leider nicht.

Mit

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

Gibt es massenweise Überseher ab 2^32. Möchte aber nicht auf x86 zurück, weil x64 so schön schnell ist.

Columbo

Ich schreibe gerade u.a. mulmod für u128 (dazu kommt wahrscheinlich bald ein Thread). Dasselbe geht mit u64 . Folgendes habe ich u.a. als non-assembler Variante stehen:

u64 mulmod_peasant(u64 a, u64 b, u64 m) {
	u64 res = 0;
	for (b %= m; a != 0; a >>= 1) {
		if (a & 1) {
			if (res >= m - b)
				res -= m;
			res += b;
		}
		b += b - (b >= m / 2) * m;
	}
	return res;
}

Das ist jetzt aber auch guarded.

Bengo

Es kommt so oder so ab 2^32 zu Fehlern, da bei MillerRabin Zahlen entstehen, die bis zum Quadrat der Testzahl gehen, und da es ohne BigInteger eben nur bis 2^64 geht, wirst du immer auf dieses Problem stoßen.
Mit Assemblercode kann man vielleicht bei der modpow Funktion und dem temp = temp *temp % number, ohne dieses Quadrat auskommen. Bringt dir aber nicht so viel, weil du sowieso früher oder später BigInteger nutzen wirst

Erhard Henkes

Ja, mit mulmod_peasant kommt es auch sofort zu Übersehern. Also ran an BigInt.

volkard

Ich finde, IsPrime muss erstmal beschleunigt werden.

int main() {
	ifstream luegner("2-SPRP-2-to-64.txt");//https://miller-rabin.appspot.com/
	UInt64 n;
	size_t count=0;
	while(luegner>>n) {
		++count;
		cout<<count<<' '<<n<<endl;
	}
}

Ausgabe letzte Zeile:
31894014 18446744066047760377
Leider mag man so eine große Lügner-Tabelle nicht im RAM halten.

while(luegner>>n) {
		if(isSPRP(n,3)) continue;
		++count;
		cout<<count<<' '<<n<<endl;
	}

1501720 18446732893888604471
Jup, dahin geht der Weg.

while(luegner>>n) {
		if(n%3==0) continue;
		if(!isSPRP(n,3)) continue;
		++count;
		cout<<count<<' '<<n<<endl;
	}

1501720 18446732893888604471
Bringt nix.

while(luegner>>n) {
		if(n%3==0) continue;
		if(n%5==0) continue;
		if(n%7==0) continue;
		if(n%11==0) continue;
		if(n%13==0) continue;
		if(n%17==0) continue;
		if(n%19==0) continue;
		if(n%23==0) continue;
		if(n%29==0) continue;
		if(n%31==0) continue;
		if(n%37==0) continue;
		if(n%41==0) continue;
		if(n%43==0) continue;
		if(!isSPRP(n,3)) continue;
		++count;
		cout<<count<<' '<<n<<endl;
	}

1501438 18446732893888604471
Bringt nix.

131157 18446602774641402961
Ok, wenig genug, um sie im RAM zu halten.

Und dann im Array der fantastisch guten Lügner mit binary search schauen.

Oder die harte Tour:

while(luegner>>n) {
		if(!isSPRP(n,3)) continue;
		if(!isSPRP(n,5)) continue;
		if(!isSPRP(n,7)) continue;
		if(!isSPRP(n,11)) continue;
		if(!isSPRP(n,13)) continue;
		if(!isSPRP(n,17)) continue;
		if(!isSPRP(n,19)) continue;
		++count;
		cout<<count<<' '<<n<<endl;
	}

1 341550071728321
2 84983557412237221
3 230245660726188031
4 1134931906634489281
5 1144336081150073701
6 1167748053436849501
7 1646697619851137101
8 3825123056546413051
9 4265186605968234451
10 5474093792130026911
11 7033671664103127781
12 7361235187296010651
13 8276442534101054431
14 14688059738864848381
15 16043083915816662841

So in dieser Richtung sollte man gehen und einen deterministischen Primzahlentester bauen, der bis 2^64 sicher nicht lügt.

volkard

volkard schrieb:

Erhard Henkes schrieb:

Ich weiß nicht genau, wo es anfängt, aber es ist unterhalb 2^64 / 2.

Ja, ab dort irgendwo hab ich lauter Überseher. Und Überseher darf es gar nicht geben.
Habs so gefixt.

/*bool IsPrimeMillerRabinOptimized(uint64_t number, uint64_t base){
	return isSPRP(number,base);
}*/

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

Sorry.
Die obere Version, die auskommentierte war ok. Die Untere hat Überseher. Liegt gar nicht an mulmod.

Columbo

PS: ICC würde ich definitiv mal probieren. Ist echt der Hammer was der noch manchmal rausholt.

Erhard Henkes

Danke! Wichtiger Hinweis.

Bei mir sieht es mit der BigInteger Library (nicht so schnell wie gmp, das schaffe ich aber auch noch) momentan so aus:

#include <iostream>
#include "BigIntegerLibrary.hh"

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 
BigUnsigned mulmod(BigUnsigned a, BigUnsigned b, BigUnsigned mod)
{
	BigUnsigned x = 0;
	BigUnsigned y = a % mod;

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

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

bool isSPRP(BigUnsigned n, BigUnsigned a)
{
	if (a%n == 0) 
		return true;
	BigUnsigned d = n - 1;
	BigUnsigned ad;
	BigUnsigned s = 0;

	// break down n-1 into d*(2^s). Linux ffs() can be better
	while ((d & 1) == 0 )
	{
		++s; d >>= 1;
	}

	ad = (powmod_Volkard(a%n, d, n)); // (a^d) mod n

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

bool IsPrimeMillerRabinOptimized(BigUnsigned number, BigUnsigned base)
{
	return isSPRP(number,base);
}

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

int main()
{	

	for (BigUnsigned i(4294967295); ; ++i)
	{
		if (i % 1000 == 0)
			cout << "i = " << i << endl;

		if (IsPrimeMillerRabinOptimized(i, 2) && !IsPrimeDivisionTest(i))
			cout << i << " Luegner" << endl; // Miller-Rabin Lügner
		if (!IsPrimeMillerRabinOptimized(i, 2) && IsPrimeDivisionTest(i))
			cout << i << " Ueberseher" << endl; // Miller-Rabin Primzahlen-Überseher
	}
	wait();
}

Im Bereich hinter 2^32 keine Überseher. Der Ansatz mit BigInteger ist folglich hilfreich.

Nun noch mit:

for (BigUnsigned i= stringToBigUnsigned("9223372036854776000"); ; ++i)

Auch keine keine Überseher in diesem Bereich, in dem es vor Kurzem richtig schlimm wurde. Das ist fein. Dann können wir diesen Divisionstest wieder abschalten und die Goldbach-Vermutung weiter strapazieren.

Erhard Henkes

#include <iostream>
#include <iomanip>
#include <cstdint>
#include <cstdlib>
#include <algorithm>
#include <cmath> 
#include <ctime>
#include <vector>
#include "BigIntegerLibrary.hh"

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 
BigUnsigned mulmod(BigUnsigned a, BigUnsigned b, BigUnsigned mod)
{
	BigUnsigned x = 0;
	BigUnsigned y = a % mod;

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

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

bool isSPRP(BigUnsigned n, BigUnsigned a)
{
	if (a%n == 0) 
		return true;
	BigUnsigned d = n - 1;
	BigUnsigned ad;
	BigUnsigned s = 0;

	// break down n-1 into d*(2^s). Linux ffs() can be better
	while ((d & 1) == 0 )
	{
		++s; d >>= 1;
	}

	ad = (powmod_Volkard(a%n, d, n)); // (a^d) mod n

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

bool IsPrimeMillerRabinOptimized(BigUnsigned number, BigUnsigned base)
{
	return isSPRP(number,base);
}

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)
{
	if (number <= primeLimit)
		return primes[number.toUnsignedLong()]; // lookup from bitset (in the RAM)
	return IsPrimeMillerRabinOptimized(number, 2);							   
}

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

int main()
{
	BigUnsigned startNumber = stringToBigUnsigned("1000000000000000000000");
	BigUnsigned endNumber   = stringToBigUnsigned("1000000000000000010000");

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

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

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

	// sieving loop
	uint64_t iMax = (uint64_t)sqrt((double)primeLimit) + 1;
	for (uint64_t i = 2; i<iMax; ++i)
	{
		if (primes[i])
		{
			for (uint64_t j = 2 * i; j<primeLimit + 1; j += i)
			{
				primes[j] = false;
			}
		}
	}

	time(&endTime);

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

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

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

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

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

		zeit2_old = zeit2;
		zeit2 = clock();

		if (i % 10 == 0)
			cout << "time: " << 1000 * (zeit2 - zeit2_old) / CLOCKS_PER_SEC << " ms" << endl;
	}
	wait();	
}

`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

1000000000000000000000 101 time: 641 ms

1000000000000000000010 181 time: 1035 ms

1000000000000000000020 191 time: 1053 ms

1000000000000000000030 131 time: 752 ms

1000000000000000000040 211 time: 1107 ms

1000000000000000000050 151 time: 854 ms

1000000000000000000060 173 time: 942 ms

1000000000000000000070 241 time: 1279 ms

1000000000000000000080 181 time: 983 ms

1000000000000000000090 191 time: 1005 ms

1000000000000000000100 271 time: 1392 ms

1000000000000000000110 211 time: 1096 ms`

Nicht schlecht für den Bereich einer Trilliarde. Ich gebe nur jede fünfte berechnete Geradzahl aus. Man sieht, dass bei 10^21 die Primzahlendichte noch sehr hoch ist. Der Abstand ist immer noch kleiner 1000.

Erhard Henkes

Selbst bei einer Vigintillion (10^120) macht der Milli-Miller-Rabin noch Freude, und immer noch diese hohe Primzahlendichte, also die erste Zahl des Paares unter 10^3. Beeindruckend.

`Now Goldbach first prime number pair will be detected

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 173 time: 133982 ms

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002 709 time: 433121 ms

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004 613 time: 378876 ms`