Goldbach Vermutung

[Erhard Henkes]

Hab nur 2 Kerne

Dann stellst Du die numberOfThreads auf 2, 3 oder 4 (dein OS kümmert sich um das MT Scheduling) ein. Ich baue das noch um, sodass die richtige Threadanzahl verwendet wird und diese auch immer etwas zu schaffen haben.
Probiere momentan verschiedene Abläufe aus.

Wichtige Frage: Kann man die Gesamtauslastung (vom Prg aus) des gesamten Programms auf einen gewissen Wert (z.B. 50 oder 80%) einstellen?

[Erhard Henkes]

Hier ein erster Test für 10^240 (!):
`

Now Goldbach probable (no liar list available) first prime number pair will be detected:

...04 97

...16 109

...10 103

...08 101

...14 107

...20 113

...06 2357

...18 3527

...00 2351

...02 2039

...12 5279

total time: 3654 s`

Zeitbestimmend bei der Primzahl-Paarsuche ist die Zahl 10^240 + 12, da wir erst bei 5279 fündig werden. Hier könnten vlt. noch mehr Primzahlen (bisher: 3 ... 101) beim Pre-Divisions-Check helfen.

Man sieht, dass nun die Primzahlabstände (nach unten) von der untersuchten Geradzahl nun auch mehrfach im Bereich über 1000 liegen (bei 10^120 war es nur eine von elf Zahlen). Gibt es hier eine Regel für den Abstand zwischen zwei Primzahlen in Abhängigkeit vom untersuchten Zahlenbereich?

Dieser Zahlenbereich 10^240 stellt somit in etwa die Grenze dar für das aktuelle Programm. Ich gehe davon aus, dass Goldbach mit seiner Vermutung recht behalten wird. Das Programm ist auch ein hübscher CPU-Stresstest zum Testen neuer PCs. Einfach den Bereich X bis X + 100 (X: sehr große Zahl) berechnen lassen und Zahl der Threads entsprechend auf 51 einstellen.

[unskilled]

Wichtige Frage: Kann man die Gesamtauslastung (vom Prg aus) des gesamten Programms auf einen gewissen Wert (z.B. 50 oder 80%) einstellen?

macht man ja üblicherweiße über OS. alternativ kann man es ja mit nem sleep machen und zeit/tickcount nehmen

edit: so gehts ein wenig bequemer, wenn ichs richtig überflogen habe:
https://technet.microsoft.com/en-us/library/ff384148.aspx?f=255&MSPPError=-2147217396

[Erhard Henkes]

Ich habe doch noch mal 10^300 geprüft:

`..04 73

..10 79

..14 83

..02 71

..20 89

..06 353

..12 359

..00 347

..16 6427

..08 3767

..18 12149

total time: 14226 s

`

Da in dieser Region der Abstand zur nächsten Primzahl (nach unten) in einem Fall (10^300 + 18) auf fünfstellige Zahlen anwächst, braucht die Berechnung fast 4 Std. Das ruft nun wirklich neue Algos auf den Plan.

[sebi707]

Irgendwas stimmt mit deinem Multithreading nicht richtig. Ich hab die Diskussion hier nur so halb verfolgt, wollte aber mal das Programm ausprobieren. Neuere gcc Versionen (4.9) beschweren sich über diese Zeile:

t.emplace_back(PrimePairFound, primes, i+2*x, BigUnsigned((unsigned long)primeLimit), refSetLiar, retVal + x);

Mit einer älteren Version (4.4) compiliert es immerhin. Ausschnitt der Fehlermeldung:

main.cpp:336:125:   required from here
/opt/rh/devtoolset-3/root/usr/include/c++/4.9.2/functional:1665:61: error: no type named ‘type’ in ‘class std::result_of<void (*(std::vector<bool>, BigUnsigned, BigUnsigned, std::unordered_set<long unsigned int>, bool*))(std::vector<bool>&, BigUnsigned, BigUnsigned, std::unordered_set<long unsigned int>&, bool*)>’
       typedef typename result_of<_Callable(_Args...)>::type result_type;

Ich denke das Problem ist, dass deine PrimePairFound Funktion einige Parameter als Referenz nimmt. Wenn ich den Code so ändere compiliert es auch mit gcc 4.9:

t.emplace_back(PrimePairFound, std::ref(primes), i+2*x, BigUnsigned((unsigned long)primeLimit), std::ref(refSetLiar), retVal + x);

Aber selbst wenn es compiliert kriege ich einen Segfault in genau dieser Zeile.

[Erhard Henkes]

Im MS VS 2015 läuft es ohne Probleme. Daher ändere ich da gar nix.

[sebi707]

Erhard Henkes schrieb:

Im MS VS 2015 läuft es ohne Probleme. Daher ändere ich da gar nix.

Nicht gerade die beste Einstellung. Ein Segfault kommt auch nicht einfach so. Irgendwas geht dort wohl ziemlich schief. Warum es bei VS2015 nicht passiert oder man nichts bemerkt weiß ich nicht. Vielleicht versuche ich nacher mal den Fehler einzugrenzen.

[sebi707]

Ok sorry. Das Problem saß mal wieder vorm Bildschirm. Seit wann muss man pthread dazu linken wenn man std::thread benutzt? Naja läuft jetzt jedenfalls wobei das mit dem std::ref trotzdem nötig war.

[Erhard Henkes]

OK, über den wrapper std::ref kann man reden. Das ist eher eine Verbesserung. Danke für den Hinweis.

t.emplace_back(PrimePairFound, std::ref(primes), i+2*x, BigUnsigned((unsigned long)primeLimit), std::ref(refSetLiar), retVal + x);

[Columbo]

sebi707 schrieb:

Ok sorry. Das Problem saß mal wieder vorm Bildschirm. Seit wann muss man pthread dazu linken wenn man std::thread benutzt?

Immer.

Naja läuft jetzt jedenfalls wobei das mit dem std::ref trotzdem nötig war.

Jo. Die result_of -Zeile hätte sonst nicht kompiliert, schließlich werden die Argumente mit declval simuliert, welches standardmäßig rvalues produziert.

[sebi707]

Arcoth schrieb:

sebi707 schrieb:

Ok sorry. Das Problem saß mal wieder vorm Bildschirm. Seit wann muss man pthread dazu linken wenn man std::thread benutzt?

Immer.

OK ich hab bisher unter Linux noch nichts mit std::thread gemacht. Entweder direkt POSIX Threads benutzt oder OpenMP. Unter Ubuntu kriegt man immerhin eine Exception die erklärt das Exceptions deaktiviert sind. Unter Redhat gibts nen Segfault 20 Funktionen tief in irgendwelchen std::thread Innereien.

[sebi707]

sebi707 schrieb:

Unter Ubuntu kriegt man immerhin eine Exception die erklärt das Exceptions deaktiviert sind.

Das Threads deaktiviert sind natürlich.

[Erhard Henkes]

Ich habe am Beispiel 1^120 einen erweiterten Divisions-Precheck (3 ... 199 anstelle 3 ... 101) ausprobiert:

// division test with first prime numbers
	if (
		/*(number % 2 == 0) || */
		(number % 3 == 0) || (number % 5 == 0) || (number % 7 == 0) || (number % 11 == 0) ||
		(number % 13 == 0) || (number % 17 == 0) || (number % 19 == 0) || (number % 23 == 0) ||
		(number % 29 == 0) || (number % 31 == 0) || (number % 37 == 0) || (number % 41 == 0) ||
		(number % 43 == 0) || (number % 47 == 0) || (number % 53 == 0) || (number % 59 == 0) ||
		(number % 61 == 0) || (number % 67 == 0) || (number % 71 == 0) || (number % 73 == 0) ||
		(number % 79 == 0) || (number % 83 == 0) || (number % 89 == 0) || (number % 97 == 0) ||
		(number % 101 == 0) || (number % 103 == 0) || (number % 107 == 0) || (number % 109 == 0) ||
		(number % 113 == 0) || (number % 127 == 0) || (number % 131 == 0) || (number % 137 == 0) ||
		(number % 139 == 0) || (number % 149 == 0) || (number % 151 == 0) || (number % 157 == 0) ||
		(number % 163 == 0) || (number % 167 == 0) || (number % 173 == 0) || (number % 179 == 0) ||
		(number % 181 == 0) || (number % 191 == 0) || (number % 193 == 0) || (number % 197 == 0) ||
		(number % 199 == 0) 
		)
	{
		return false;
	}

`Now Goldbach probable (no liar list available) first prime number pair will be detected:

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000018 191

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000008 181

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000006 179

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 173

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004 613

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000016 733

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000020 193

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010 619

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012 719

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002 709

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000014 1153

time: 184750 ms

total time: 184 s`

Eine deutliche Beschleunigung (184 s vs. 237 s mit 3 ... 101). Es verändert auch die Reihenfolge, in der die ersten Primzahlen als Teil einer Primzahlensumme auftauchen. Wo hier das Optimum in Abhängigkeit von der zu prüfenden Zahl liegt, ist eine interessante Frage.

Da wir in unser vector<bool> bereits die Primzahleninformationen im unteren Zahlenbereich gespeichert haben, können wir diese Information im Code nutzen:

// division test with first prime numbers
	const uint32_t DIVISOR_PRECHECK_MAX = 199;
	for (uint32_t divisor = 3; divisor <= DIVISOR_PRECHECK_MAX; ++divisor)
	{
		if (primes[divisor])
		{
			if (number % divisor == 0)
				return false;
		}
	}

Der Lookup-Vorgang kostet bei unserem Test im Bereich 1^120 in Summe ca. drei Sekunden (187 vs. 184 s). Man kann auf diese Weise aber leichter das Geschwindigkeits-Optimum (Welcher DIVISOR_PRECHECK_MAX in Abhängigkeit von startNumber?) finden.

Folgende kleine Serie habe ich getestet:
10^120:

DIVISOR PRECHECK von ... bis ...: Zeit zur Prüfung der 11 Ganzzahlen von 10^120 bis incl. 10^120 + 20 mit 11 parallelen Threads (langsamste Testschleife mit Milli-Rabin-isProbablyPrime ist zeitbestimmend):

`3 ... 101: 226 s

3 ... 199: 187 s

3 ... 997: 142 s

3 ... 9973: 98 s

3 ... 99991: 81 s

3 ... 999983: 79 s

3 ... 1999993: 84 s`

Fazit: Für Goldbach-Primzahlenpaar-Tests im Bereich 1^120 und höher kann ich 999983 als DIVISOR_PRECHECK_MAX vorteilhaft einsetzen.

[Erhard Henkes]

Ich sehe bezüglich weiterer Optimierungen folgende Felder:
- besseres BigInteger (vlt. doch gmp)
- voll ausgelastetes threading
- mulmod in asm <-- da hängt allerdings BigInteger dran
- die "loopende Tabelle im Cache" (habe die genaue Idee von TGGC noch nicht verstanden, Code wäre hilfreich)
- weitere mathematische Tricks

[sebi707]

Ich hab gerade mal was mit MPIR zusammengeklatsch. Das ist überhaupt kein Vergleich:

1...000    173
1...002    709
1...004    613
1...006    179
1...008    181
1...010    619
1...012    719
1...014   1153
1...016    733
1...018    191
1...020    193
time: 58 ms (Ja, wirklich ms)

Und das ist noch nichtmal Multithreaded... Hier der Code:

#include <gmp.h>
#include <gmpxx.h>
#include <iostream>
#include <iomanip>

using namespace std;

gmp_randstate_t rnd;

bool PrimePairFound(const mpz_class& i)
{
  bool retVal = false;

  mpz_class grenze = i / 2;
  mpz_class a = 0;
  for(a = 3; a <= grenze; a += 2)
  {
    mpz_class b = i - a;
    if(mpz_likely_prime_p(a.get_mpz_t(), rnd, 2) && mpz_likely_prime_p(b.get_mpz_t(), rnd, 2))
    {
      retVal = true;
      break;
    }
  }

  cout << i << "  " << setw(5) << a << " " << endl;
  return retVal;
}

int main()
{
  gmp_randinit_default(rnd);
  mpz_class startNumber("1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
  mpz_class endNumber("1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000020");

  clock_t start = clock();
  for(mpz_class i = startNumber; i <= endNumber; i += 2)
  {
    if(!PrimePairFound(i))
    {
      cout << "Counterevidence found for this number: " << i << endl;
    }
  }
  clock_t end = clock();
  cout << "time: " << 1000*(end - start) / CLOCKS_PER_SEC << " ms" << endl;
  cin.get();
}

Die Funktion für den Primzahltest hat zwar ein likely im Namen aber wenn ich es richtig verstanden habe ist es sehr sehr unwahrscheinlich eine Zahl als Primzahl zu erkennen wenn es eigentlich keine ist.

[Bengo]

Wikipedia sagt, der Test ist bei einer zufälligen Basis zu 25% nicht zuverlässig, skaliert aber sehr gut, nach 4 Tests sind nur noch 0,4%, nach 10 Tests, etwa 1ppm.

[sebi707]

Oder man kann auch die mpz_probab_prime_p Funktion nutzen, bei der man die Anzahl der Miller-Rabin-Test angeben kann. Wenn ich dort 100 Tests angebe braucht der Benchmark 112ms, dafür werden zusammengesetzte Zahlen mit einer Wahrscheinlichkeit von maximal 2^-200 ≈ 6,22*10^-61 als Primzahlen erkannt. Das ist schon ziemlich unwahrscheinlich...

[Erhard Henkes]

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.

Aktueller Sourcecode: http://codepad.org/puR6GPwG

[TGGC]

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.

[Finnegan]

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