Goldbach Vermutung
-
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.
-
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 dieuint32_t
unduint64_t
Typen?
-
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 dieuint32_t
unduint64_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.
-
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.
-
Ich schreibe gerade u.a.
mulmod
füru128
(dazu kommt wahrscheinlich bald ein Thread). Dasselbe geht mitu64
. 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.
-
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
-
Ja, mit mulmod_peasant kommt es auch sofort zu Übersehern. Also ran an BigInt.
-
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 16043083915816662841So in dieser Richtung sollte man gehen und einen deterministischen Primzahlentester bauen, der bis 2^64 sicher nicht lügt.
-
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.
-
PS: ICC würde ich definitiv mal probieren. Ist echt der Hammer was der noch manchmal rausholt.
-
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.
-
#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.
-
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`
-
Recht anschauliche Erklärung des Miller-Rabin-Tests: http://www.johannes-bauer.com/compsci/millerrabin/index.php?zahl=1001&iterationen=&calculate=1
http://www.iti.fh-flensburg.de/lang/krypto/algo/primtest.htmWie hoch ist eigentlich unsere Restfehlerwahrscheinlichkeit?
-
Ups ups ups!
mulmod_peasant
ist Quatsch. Hab' ich erst jetzt rausgefunden, als ich eine neue Reihe von Benchmarks machen wollte.Hier ist die korrigierte Version.
// Requires i < m u128 twiceBounded(u128 i, u128 m) { return i-(i > m/2)*m+i; } u128 mulmod_arcoth(u128 a, u128 b, u128 m) { u128 res = 0; for (b %= m; a != 0; a >>= 1) { if (a & 1) res += b - (b >= m-res)*m; b = twiceBounded(b, m); } return res; }
-
Ist mir leider nicht gelungen diese Variante für BigUnsigned zu übertragen. Schade, hätte sie gerne im Rennen antreten lassen.
-
Erhard Henkes schrieb:
Ist mir leider nicht gelungen diese Variante für BigUnsigned zu übertragen. Schade, hätte sie gerne im Rennen antreten lassen.
Naja also a&1 ist ein bitweise und, und in dem Fall gibt es 1 aus, wenn die Zahl ungerade ist und 0, wenn es gerade ist, also das gleiche wie a%2 == 1.
a >> ist a/2 (Verschieben der bits nach links)
-
So klappt es mit BigUnsigned:
// Requires i < m BigUnsigned twiceBounded(BigUnsigned i, BigUnsigned m) { BigUnsigned temp = (i > m / 2); return i - (temp) * m + i; } BigUnsigned mulmod_arcoth(BigUnsigned a, BigUnsigned b, BigUnsigned m) { BigUnsigned res = 0; for (b %= m; a != 0; a >>= 1) { BigUnsigned temp = (b >= m - res); if (a % 2 == 1) { res += b - (temp) * m; } b = twiceBounded(b, m); } return res; }
Es kompiliert, aber leider bittet das Programm das OS um den Exitus. Kann daher keinen Geschwindigkeitsvergleich aufzeigen. Daher verwende ich das bisherige mulmod, das sich gut bewährt, weiter:
//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; }
-
Jetzt beschleunigen wir das aber noch an anderer Stelle, wie bereits von euch vorgeschlagen:
bool IsPrime(vector<bool>& primes, BigUnsigned number, BigUnsigned primeLimit) { if (number <= primeLimit) return primes[number.toUnsignedLong()]; // lookup from bitset (in the RAM) if ((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)) { return false; } return IsPrimeMillerRabinOptimized(number); }
Vergleich bei 10^24:
Ohne die Abfrage number % ... == 0 mit 3 ... 43:
`Now Goldbach first prime number pair will be detected
1000000000000000000000000 257 time: 1887 ms
1000000000000000000000002 349 time: 2449 ms
1000000000000000000000004 307 time: 2182 ms
1000000000000000000000006 263 time: 1908 ms
1000000000000000000000008 311 time: 2172 ms
1000000000000000000000010 3 time: 33 ms
1000000000000000000000012 5 time: 61 ms
1000000000000000000000014 7 time: 92 ms
1000000000000000000000016 487 time: 3157 ms
1000000000000000000000018 11 time: 122 ms
1000000000000000000000020 13 time: 152 ms`
Mit dieser Abfrage vor dem Milli-Miller-Rabin-Test:
`Now Goldbach first prime number pair will be detected
1000000000000000000000000 257 time: 520 ms
1000000000000000000000002 349 time: 931 ms
1000000000000000000000004 307 time: 561 ms
1000000000000000000000006 263 time: 478 ms
1000000000000000000000008 311 time: 798 ms
1000000000000000000000010 3 time: 33 ms
1000000000000000000000012 5 time: 64 ms
1000000000000000000000014 7 time: 65 ms
1000000000000000000000016 487 time: 636 ms
1000000000000000000000018 11 time: 35 ms
1000000000000000000000020 13 time: 64 ms
`
Eine spürbare Beschleunigung.
Nun fehlt noch die Lügnerliste als unordered_set<BigUnsigned>.