Goldbach Vermutung

volkard

Erhard Henkes schrieb:

@bengo, volkard: Ist jetzt der Milli-Miller-Rabin mit den Fehlertabellen (hässlich!) der richtige Weg ab 10^11 aufwärts?

Ja, bis 2^64 ist der am besten.

Wobei leichte Anpassungen meinem Stil entsprechen würden: Ich würde die 300M-Fehlertabelle noch zusammenschrumpfen, indem ich ihm einen Test eine weitere Basis gönnen würde. Aber das ist unwesentlich.

Erhard Henkes

Ich erwarte kein Gegenbeispiel zur Goldbachschen Vermutung, aber ohne Beweis kann man nicht zu 100% sicher sein.

200000000000 prime pairs: 281856501 calc. time: 99635 ms

Die Rechenzeiten für die Primzahl-Paare gehen jetzt linear mit der Eingangszahl hoch.

Daher macht die gesamte Primzahlenpaarzählung auf diese Weise keinen Sinn mehr, sondern es wird nach dem ersten Paar abgebrochen.

Ich denke momentan über eine Mischform nach: Soweit der implementierte RAM mitspielt, mit vector<bool> oder bitset (man kennt ja die genaue Größe) prüfen, darüber dann schneller Primzahltest, damit if (IsPrime(primes, a) && IsPrime(primes, i - a)) möglichst rasch durchlaufen wird. Dann Ausbau auf MT.

Erhard Henkes

Aktueller Stand:
Strategie und Zeiten finde ich akzeptabel in dem hohen Zahlenbereich

Fragen/Probleme:
Miller-Rabin klappt irgendwie nicht schnell bei hohen Zahlen,
habe 4 u. 5 Mrd. getestet, da geht es, also kein 2^32 problem
vlt. ist da noch ein fehler? (andere quellen checken)

Wollte modulo ersetzen, aber bei Zielplattform x64 nimmt er kein __asm.

// ConsoleApplication5.cpp : Definiert den Einstiegspunkt für die Konsolenanwendung.
// single-threaded prime number generator

#include "stdafx.h"
#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;
}

// modular exponentiation
uint64_t modulo(uint64_t base, uint64_t exponent, uint64_t mod)
{
	uint64_t x = 1;
	uint64_t y = base;
	while (exponent > 0)
	{
		if (exponent % 2 == 1)
		{
			x = (x * y) % mod;
		}			
		y = (y * y) % mod;
		exponent /= 2;
	}
	return x % mod;
}

// Miller-Rabin primality test, iteration signifies the accuracy
// http://www.sanfoundry.com/cpp-program-implement-miller-rabin-primality-test/
bool IsPrimeMillerRabinCalculation(uint64_t number)
{
	if (number == 2) return true;
	if (number % 2 == 0) return false;

	uint64_t s = number - 1;

	while (s % 2 == 0)
	{
		s /= 2;
	}

	uint64_t maxIteration = 1; // variable for accuracy
	for (uint64_t it=0; it<maxIteration; ++it)
	{
		uint64_t a = rand() % (number - 1) + 1; 
		uint64_t temp = s;
		uint64_t mod = modulo(a, temp, number);

		while ((temp != number - 1) && (mod != 1) && (mod != number - 1))
		{
			mod = mulmod(mod, mod, number);
			temp *= 2;
		}

		if (mod != number - 1 && temp % 2 == 0)
		{
			return false;
		}
	}
	return true;
}

bool IsPrimeByCalculation(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]; // from bitset in RAM
	else
		return IsPrimeByCalculation(number);
		//return IsPrimeMillerRabinCalculation(number);
}

bool ZweiPrimzahlenGefundenInPrimeTable(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 (!ZweiPrimzahlenGefundenInPrimeTable(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();
}

`

Start number to be executed: 18000000000000000000

Last number to be executed: 18000000000000000010

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

18000000000000000000 79 17999999999999999921 time: 20872 ms

18000000000000000002 211 17999999999999999791 time: 21012 ms

18000000000000000004 83 17999999999999999921 time: 16571 ms

18000000000000000006 173 17999999999999999833 time: 27670 ms

18000000000000000008 1747 17999999999999998261 time: 17035 ms

18000000000000000010 89 17999999999999999921 time: 16691 ms`

Bengo

Denke deine modulo Funktion ist zu ineffizient.

Du kannst deinen exponenten in so viele Teile aufteilen, dass ein Teil noch so groß ist, wie du es noch berechnen kannst. Dann ein Teil modulo rechnen. Das Ergebnis noch hoch der Anzahl der Teile und dann noch mal modulo. Das sollte wesentlich schneller gehen.

Bei mulmod ist es ähnlich

Und ich bin mir hier zwar nicht vollständig sicher, aber glaube du nutzt den Trick mit den zweierpotenzen falsch.
Die Idee ist, man nimmt a^d und rechnet mod n;
Dann rechnen man das ergebnis mal 2 und rechnen wieder mod n;
und so weiter, bis 1 bzw -1 rausgekommen ist.

Ok eigentlich nicht mal 2, sondern hoch 2. Wenn man es abstract genug sieht, ist das sowieso alles das gleiche

Bengo

Ich würde das so implementieren:

miller(int zahl, int basis) {
  int d = zahl -1;
  while (d % 2 == 0) {
    d = d/2;

  }
  int temp = basis^d;
  if (temp % zahl == 1) {
    return ist prim;
  }

  while(temp <= basis^zahl) {
     temp = temp ^2;
     if (temp % zahl == zahl -1) {
        return ist prim;
     }
  }
  return ist nicht prim;

}

Gut das ^ gibt es jetzt nicht, das müsste man eventuell noch mal implementieren, insbeondere, dass es keinen überlauf gibt.
Und die Bedinung in der zweiten while, müsste man auch noch mal optimieren, ist aber ganz einfach zu machen.

Erhard Henkes

@Bengo: Ich habe das auf die Schnelle ausprobiert, geht bei kleinen Zahlen gut, aber bei großen Zahlen gibt es durch die Potenzen vermutlich Probleme mit den Grenzen des Zahlenraums 0 ... 2^64-1. Welche pow(...) sollte man da einsetzen? double pow(double, uint..._t) anstelle uint...? Echtes Potenzproblem.

uint64_t pow_n(uint64_t x, uint64_t n)
{
	uint64_t y=1;
	for (uint64_t i=0; i<n; ++i)
		y *= x;
	return y;
}

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

	uint64_t temp = pow_n(base, d);

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

	while (temp <= pow_n(base, number)) 
	{
		temp *= temp;

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

Bengo

Du hast ja im Prinzip nur eine Potenzrechung die kritisch ist, nämlich das basis^d; Das müsstest du im Notfall eben iterativ machen, du kannst zwischendurch immer wieder modulo zahl rechnen.
Kannst auch https://en.wikipedia.org/wiki/Exponentiation_by_squaring nehmen um es zu beschleunigen.

Da fällt mir ein, dass du statt temp *= temp, temp = (temp * temp) % zahl; machen must. und dann if(temp == zahl - 1). Das geht auch wieder ein Stück schneller.

Statt der zweiten while Schleife, macht man dann doch eher eine for schleife die so oft durchgeht, wie von der ersten while Schleife durch 2 geteilt wurde.

unskilled

Bengo schrieb:

temp = (temp * temp) % zahl;

temp2 = temp%zahl;
temp2 *= temp2;
temp = temp2%zahl;

ist ein div mehr (der imho nicht wirklich notwendig ist) aber dafür kein überlauf.

Erhard Henkes

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

Bengo

unskilled schrieb:

Bengo schrieb:

temp = (temp * temp) % zahl;

temp2 = temp%zahl;
temp2 *= temp2;
temp = temp2%zahl;

ist ein div mehr (der imho nicht wirklich notwendig ist) aber dafür kein überlauf.

Der Überlauf kann in dem Fall nicht passieren, da im vorherigen Schleifendurchlauf schon mod zahl gerechnet wurde.

Also die pow Funktion kann man so implementieren:

uint64_t mod_pow(uint64_t base, uint64_t expo, uint64_t mod) {
    uint64_t x = 1;
    while (expo > 1) {
        if (expo % 2 == 1) {
            x *= base;
            expo--;
        }
        x *= x;
        expo /= 2;
        x = x % mod;
    }
    return (x*base) % mod;
}

und die andere Funktion:

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 = mod_pow(base, d, number);

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

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

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

Erhard Henkes

Bengo: Danke! Echt klasse. Läuft bis 1 Mrd. gut
bei größer 2^32 nicht mehr

uint64_t temp = mod_pow(base, d, number); <== kein mod als parameter?

`Start number to be executed: 4000000000

Last number to be executed: 4000000010

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

In 9 s we found 50847534 prime numbers (2 to 1000000000 ).

Now Goldbach first prime number pair will be detected

4000000000 241903391 3758096609 time: 4851 ms

4000000002 233510801 3766489201 time: 4692 ms

4000000004 95035861 3904964143 time: 2012 ms

4000000006 40126469 3959873537 time: 896 ms

4000000008 241903399 3758096609 time: 4856 ms

4000000010 107522011 3892477999 time: 2254 ms`

aber:
`Start number to be executed: 5000000000

Last number to be executed: 5000000010

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

In 10 s we found 50847534 prime numbers (2 to 1000000000 ).

Now Goldbach first prime number pair will be detected`

nix passiert mehr ...

Columbo

So geht das.

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

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

		x = sqrmod(x, mod);
	}

	return res;
}

Erhard Henkes

sehr gut, jetzt gehen auch 5 Mrd.

`Start number to be executed: 5000000000

Last number to be executed: 5000000010

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

In 9 s we found 50847534 prime numbers (2 to 1000000000 ).

Now Goldbach first prime number pair will be detected

5000000000 97 4999999903 time: 1 ms

5000000002 251 4999999751 time: 2 ms

5000000004 67 4999999937 time: 1 ms

5000000006 103 4999999903 time: 2 ms

5000000008 71 4999999937 time: 1 ms

5000000010 73 4999999937 time: 1 ms

`

auch 10^18 noch beeindruckend schnell:

`

Start number to be executed: 1000000000000000000

Last number to be executed: 1000000000000000010

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

In 9 s we found 50847534 prime numbers (2 to 1000000000 ).

Now Goldbach first prime number pair will be detected

1000000000000000000 137 999999999999999863 time: 4 ms

1000000000000000002 139 999999999999999863 time: 5 ms

1000000000000000004 37 999999999999999967 time: 3 ms

1000000000000000006 503 999999999999999503 time: 9 ms

1000000000000000008 41 999999999999999967 time: 4 ms

1000000000000000010 43 999999999999999967 time: 3 ms`

Bei 10^19 ist der Ofen allerdings aus. Da kommt nix mehr.
Wohl zu nah an 18.446.744.073.709.551.616

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

`Now Goldbach first prime number pair will be detected

1000000000000000000 137 999999999999999863 time: 4230 ms

1000000000000000002 139 999999999999999863 time: 4253 ms

1000000000000000004 37 999999999999999967 time: 4194 ms

1000000000000000006 503 999999999999999503 time: 4228 ms

1000000000000000008 41 999999999999999967 time: 4249 ms

1000000000000000010 43 999999999999999967 time: 4229 ms`

Die Resultate waren also alle OK, kostet gleich 4 sec pro Ganzzahl.

Bengo

Hab ein paar blöde Fehler gemacht. Also Arcoth hat völlig recht, meine implementierung dieser Funktion ist völliger schwachsinn und funktionert nur manchmal richtig.

Bei der Miller Funktion hab ich auch was vergessen, muss so heißen:

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 = mod_pow(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;
}

überarbeitet: (funktionert jetzt auch, habs nachgeprüft, auch die miller funktion)

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 = x * aktpot;
        }
        aktpot = (aktpot * aktpot) % mod;
        expo = expo / 2;
    }
   return x % mod;
}

Bengo

Erhard Henkes schrieb:

Bei 10^19 ist der Ofen allerdings aus. Da kommt nix mehr.
Wohl zu nah an 18.446.744.073.709.551.616

mmh schwach abgeschätzt würde ich sagen, number^2 muss noch kleiner als 18.446.744.073.709.551.616 sein, sonst können falsche Ergenisse kommen.

Denke jetzt muss eine big integer lib ran

Erhard Henkes

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

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