Goldbach Vermutung

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.

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

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.

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

Erhard Henkes

Das hier läuft schon deutlich fixer:

BigUnsigned mulmod(BigUnsigned a, BigUnsigned b, BigUnsigned mod)
{
	/////////////
	QueryPerformanceCounter((LARGE_INTEGER*)&startCount);
	/////////////
	BigUnsigned x = 0;
	BigUnsigned y = a % mod;

	while (b > 0)
	{
		if (((uint32_t)b.getBlock(0)) & 1)
		{
			x += y;
			x %= mod;			
		}

		y <<= 1;
		y %= mod;

		b >>= 1;		
	}
	/////////////
	QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
	sumTime_Function += (lastCount - startCount);
	/////////////
	return x;
}

Finnegan

Erhard Henkes schrieb:

Das hier läuft schon deutlich fixer:

Cool-O-Matik!

Gibt wahrscheinlich noch ne Reihe andere Stellen im Code, die davon profitieren würden. Ist nur ein bissl blöd dass das der Code-Lesbarkeit schadet. Das ist für mich eher so C-Hacker-Stil von jemandem der grad von Assembler kommt

Erhard Henkes

Deine Variante ist 8% schneller, daher baue ich sie gerne ein. Insgesamt ein großer Schritt.

Also neuer Sourcecode: http://codepad.org/ReYuzGjk (nur mulmod)

Erhard Henkes

... auch an anderen Stellen verändert: http://codepad.org/f5zbHebm

@Finnegan: Danke!

mulmod vorher: 70% jetzt: 45 % (bei 10^54)

Es lohnt sich also, weiter an dem Interface BigInteger und mulmod zu tunen.

Finnegan

Erhard Henkes schrieb:

mulmod vorher: 70% jetzt: 45 % (bei 10^54)

Gern. Schön wenn man so ne Funktion hat, die im Profiler leuchtet wie ein Weihnachtsbaum. Da ist das Optimieren leicht.
Leider hat man auch oft mit Programmen zu tun wo alle Funktionen irgendwie gleich "langsam" sind, das ist immer ein Krampf.
...aber bei so ner schönen fetten spanischen Galleone schraub ich mir glatt mein Holzbein an und leg' ne Überstunde ein

Erhard Henkes

Vielleicht könnte man auch eine eigene MyBigInteger-Klasse auf Basis uint64_t anstelle uint32_t basteln, nur mit den Funktionen, die man hier braucht. Dann könnte man das auf Geschwindigkeit genau dafür trimmen.

Jodocus

Was ist das Ziel des ganzen hier? Willst du so weit kommen, wie du nur kannst? Dann würde ich jetzt eher aufhören, zu mikrooptimieren, sondern mir die BigInt-Implementierung ansehen. BigInts gibt's viele, aber richtig gute nicht. Wenn du's selber machen willst, dann würde ich sie auf doppelter Wortbreite basieren lassen und ebenfalls im Zweierkomplement darstellen.

Die Addition muss mit ziemlicher Sicherheit in Assembler gemacht werden, da es dort Mnemonics wie ADC gibt, die einem das explizite Carryen sparen (sofern der Compiler es nicht schon so generiert).

Danach Multiplikation: da kannst du mit Karazuba oder gleich mit Toom-Cook anfangen. Oder du machst vorher etwas Mathe, lädst dir die Fast Fourier Transform runter und implementierst Schönhage-Strassen. Die Algorithmen lassen sich auch allesamt parallelisieren + cache-optimieren.

DJohn

Hallo

Erhard Henkes schrieb:

Ich sehe bezüglich weiterer Optimierungen folgende Felder:
...

Wenn ich den Code richtig verstehe, beginnt Ihr bei jeder Zahl wieder vorn mit den Test potentieller Primzahlpärchen. Bei den großen Zahlen wie 10^20 oder sogar 10^240 sind die gefundenen großen Primzahlen doch viel zu wertvoll, um sie gleich wieder wegzuwerfen.

1.) Wenn ich z.B. für x1 = 13 + (x1-13) gefunden habe, weil x1-13 eine Primzahl ist, dann kann ich für alle Primzahlen größer 13 aus meiner Primzahltabelle weiter PrimePairs bilden und die aus der Liste der zu prüfenden Zahlen streichen.
Also:
x3 = (x1+4) = 17 + (X1-13)
x4 = (x1+6) = 19 + (X1-13)
u.s.w.
2.) für x2 = (x1+2) passt das noch nicht, aber auch hier gibt es eine Optimierung:
beim Prüfen von x1 habe ich für x1-3, x1-5, x1-7, x1-11 schon festgestellt, dass das keine Primzahlen sind. Für x2 verschiebt sich das um 2, somit kann ich mir ein paar Primzahltests gleich sparen: von x2-5, x2-7, x2-13 weiß ich schon, dass das keine Primzahlen sind. Hierfür kann man sich intern eine Tabelle halten für alle ungeraden Zahlen von endNumber-3 (Index 0) bis startNumber-n mit jeweils eine Flag für "Prime", "not Prime", oder "Unbestimmt" (evtl. noch ein "wird gerade berechnet" fürs parallelisieren). Natürlich muss die Tabelle irgendwie begrenzt werden, aber ich schätze mal, dass man die meisten PrimPairs findet, bevor die Tabelle unhandlich groß wird.

Finnegan

Erhard Henkes schrieb:

Vielleicht könnte man auch eine eigene MyBigInteger-Klasse auf Basis uint64_t anstelle uint32_t basteln, nur mit den Funktionen, die man hier braucht. Dann könnte man das auf Geschwindigkeit genau dafür trimmen.

Ich würde so frech sein und behaupten, dass der Compiler das schon macht soweit er kann. Zumindest bei dem was ist schon in von Visual C++ generiertem Maschinencode gesehen habe, würde ich sogar davon ausgehen, dass er, wenn z.B. die entsprechende Architektur aktiviert ist (bspw. -arch:avx oder -arch:avx2 - SSE2 ist per default an), er sogar schon entsprechende Vektor-Instruktionen erzeugt... aber probieren geht über studieren. Hab zwar jetzt keine Zeit mehr dafür, aber ein kurzes Mini-Testprogramm und dann im Debug-Modus mit Code-Optimierungen und mit Breakpoint auf die Operation auf "Show Disassembly" gehen, sollte Aufschluss geben (xmm/ymm-Register und vpadd*-Instruktionen wären z.B. bei Addition ein Hinweis... wichtig bei so einem Test: Variablen mit denen man sowas testet sollten nicht durch irgendwelche Konstanten/Literale gefüllt werden, sonst rechnet der das Ergebnis einfach beim Kompilieren aus und baut es direkt in die EXE ein [Stichw. Konstantenfaltung] ... ich hab mir da immer mit dem volatile -Keyword ausgeholfen wenn ich wirkliche Instruktionen sehen wollte).

Falls man es wirklich selbst implementieren will, könnte es allerdings mehr bringen, dem CPU-Cache besser entgegenzukommen. Da muss man allerdings einiges an Hirnschmalz reinstecken, auch wenn die Grundregel simpel ist: Möglichst linear auf den Speicher zugreifen und keine wilden Pointer-Sprünge quer durch den RAM. Daher habe ich auch vorhin bemerkt dass die BigInts alle auf dem Heap liegen (theoretisch kreuz- und quer verteilt). Je nach deren grösse macht es Sinn die in einem Array auf dem Stack zu haben, dann sind z.B. die lokalen BigInts auf denen man gerade arbeitet im (Stack-)RAM "nah" beieinander und man kann besser von den CPU Cache-Leveln profitieren. Die std::strings haben wegen sowas z.B. die SSO (Small String Optimization). Kurze Strings liegen auf dem Stack, erst wenn sie länger werden wird der Heap genommen. Eine richtig gute BigInt-Library würde das wahrscheinlich ähnlich machen.

Finnegan

P.S.: Noch eine einfachere Möglichkeit, mit der man noch ein bisschen rausquetschen könnte: Whole Program Optimization und Link Time Code Generation. Wenn die cpp-Dateien der BigInt-Lib Teil des Projekts sind kann er dann auch aus diesen Funktionen inlinen oder die Nutzung der CPU-Register für das ganze Programm optimieren. Es gibt Fälle wo sowas durchaus nochmal um die 10% bringen kann, wenn mn Glück hat