Carmichael Zahlen (Visualisierung und lückenlose Erzeugung)
-
Vergleicht man die aus Fermats kleinem Satz erzeugten Zahlen mit "gesiebten" Primzahlen, so bilden sich als Unterschied die "Carmichael Zahlen" heraus. Dies soll folgendes Programm (nicht optimiert) visualisieren.
#include <iostream> #include <vector> #include <cmath> #include "MyBigUnsigned.h" // http://henkessoft.de/Sonstiges/MyBigUnsigned.rar // simplified, based on https://mattmccutchen.net/bigint/ using namespace std; int primeLimit = 3000; int main() { ////////////////////////// Generate primes lookup table /////////////////////////////////////////////////////////// vector<bool> primes(primeLimit + 1, true); // calculated primes primes.at(0) = primes.at(1) = false; // Erastothenes sieving loop int iMax = (int)sqrt((double)primeLimit) + 1; for (int i = 2; i < iMax; ++i) if (primes.at(i)) for (int j = 2 * i; j < primeLimit + 1; j += i) primes.at(j) = false; ////////////////////////// Fermat's little theorem /////////////////////////////////////////////////////////////// for (int p = 2; p <= primeLimit; ++p) { int pro = 0, contra = 0; for (BigUnsigned a = 2; a <= 5; ++a) { BigUnsigned result = a; for (BigUnsigned i = 1; i < p; ++i) result *= a; result %= p; if (result != a%p) ++contra; }//for a if (contra == 0) { cout << "\t" << p; if (primes.at(p) == false) cout << " <== Carmichael number"; } }//for p }
output: http://pastebin.com/SAkhkiGz
Hier die schnellere Variante mit mpirxx.h, mulmod und powmod (Volkard), nur Ausgabe von Carmichaelzahlen mit Zeiten:
#define NOMINMAX // due to conflict with max() #include <windows.h> // warning: no C++ standard! #include <iostream> #include <vector> #include <cmath> #include <cstdint> #include "mpirxx.h" using namespace std; int primeLimit = 1000000; mpz_class mulmod(mpz_class a, mpz_class b, const mpz_class& mod) { mpz_class x = 0; mpz_class& y = a; while (b > 0) { if ((b & 1) == 1) { x += y; x %= mod; } y <<= 1; y %= mod; b >>= 1; } return x; } mpz_class powmod_Volkard(mpz_class base, mpz_class exp, mpz_class modul) { mpz_class a1 = base, z1 = exp, x = 1, m = modul; while (z1 != 0) { while ((z1 & 1) == false) { z1 >>= 1; a1 = mulmod(a1, a1, m); } --z1; x = mulmod(x, a1, m); } return x; } int main() { uint64_t frequency; QueryPerformanceFrequency((LARGE_INTEGER*)&frequency); cout << "cpu frequency: " << frequency << " Hz" << endl; cout << "Carmichael Numbers:" << endl << endl; ////////////////////////// Generate primes lookup table /////////////////////////////////////////////////////////// vector<bool> primes(primeLimit + 1, true); // calculated primes primes.at(0) = primes.at(1) = false; int iMax = (int)sqrt((double)primeLimit) + 1; for (int i = 2; i < iMax; ++i) if (primes.at(i)) for (int j = 2 * i; j < primeLimit + 1; j += i) primes.at(j) = false; // Erastothenes sieve ////////////////////////// Fermat's little theorem /////////////////////////////////////////////////////////////// uint64_t startCount = 0; QueryPerformanceCounter((LARGE_INTEGER*)&startCount); for (int p = 2; p <= primeLimit; ++p) { int pro = 0, contra = 0; for (mpz_class a = 2; a <= 5; ++a) { mpz_class result = powmod_Volkard(a, p, p); if (result != a%p) ++contra; }//for a if (contra == 0) { if (primes.at(p) == false) { uint64_t lastCount = 0; QueryPerformanceCounter((LARGE_INTEGER*)&lastCount); uint64_t delta = lastCount - startCount; cout << "p = " << p; cout << " \ttime: " << 1000 * delta / frequency << " ms" << endl; } } }//for p }
output:
cpu frequency: 3127001 Hz Carmichael Numbers: p = 561 time: 74 ms p = 1105 time: 180 ms p = 1729 time: 320 ms p = 2465 time: 501 ms p = 2821 time: 595 ms p = 6601 time: 1714 ms p = 8911 time: 2480 ms p = 10585 time: 3057 ms p = 15841 time: 5019 ms p = 29341 time: 10517 ms p = 41041 time: 15726 ms p = 46657 time: 18320 ms p = 52633 time: 21152 ms p = 62745 time: 26114 ms p = 63973 time: 26746 ms p = 75361 time: 32327 ms p = 101101 time: 45857 ms p = 115921 time: 53933 ms p = 126217 time: 59688 ms p = 162401 time: 80152 ms p = 172081 time: 85723 ms p = 188461 time: 95459 ms p = 252601 time: 134760 ms p = 278545 time: 150967 ms p = 294409 time: 160917 ms p = 314821 time: 173788 ms p = 334153 time: 186319 ms p = 340561 time: 190402 ms p = 399001 time: 229582 ms p = 410041 time: 236820 ms p = 449065 time: 263149 ms p = 488881 time: 290746 ms p = 512461 time: 307550 ms p = 530881 time: 320544 ms p = 552721 time: 335307 ms p = 656601 time: 409225 ms p = 658801 time: 410696 ms p = 670033 time: 418498 ms p = 721801 time: 456352 ms <-- wrong p = 748657 time: 476039 ms p = 825265 time: 533049 ms p = 838201 time: 542607 ms p = 852841 time: 553756 ms p = 873181 time: 568716 ms <-- wrong p = 997633 time: 664459 ms
Interessant ist, dass die Zahlen 721801, 873181 in folgender Tabelle nicht auftauchen: https://de.wikibooks.org/wiki/Pseudoprimzahlen:_Tabelle_Carmichael-Zahlen oder http://www.mathe.tu-freiberg.de/~hebisch/skripte/zahlenth/zahlenth.pdf
Ich verwende die Ganzzahlen a = 2,3,4,5. Bei weniger Ganzzahlen (z.B. nur 2,3,4) erscheinen bereits unter 100000 zusätzliche Carmichael "Lügner". Die interessante Frage ist, wieviele Ganzzahlen muss man einsetzen, um sicher sein zu können, nur "echte" Carmichael-Zahlen zu generieren?
Ein Versuch mit 2, ..., 13 "killt" die beiden "Lügner". Man muss nun heraus finden, welche Ganzzahlen wichtig sind für diese beiden "Contras".
Versuchsergebnis:
721801 wird durch a=11 und a=13 entlarvt.
873181 wird durch a=7 und a=13 entlarvt.Die 13 hat es also in sich und gibt den Lügnern contra.
Wer mit Fermats kleinem Satz spielen will: http://codepad.org/V0UOd7oN
-
for (int p = primeStart; p <= primeLimit; ++p) { int pro = 0, contra = 0; for (BigUnsigned a = 2; a < 4; ++a) // 2,3 { BigUnsigned result = powmod(a, p, p); if (result != a%p) ++contra; } if (p > 2700) { BigUnsigned a = 5; BigUnsigned result = powmod(a, p, p); if (result != a%p) ++contra; } if (p > 721800) { BigUnsigned a = 13; BigUnsigned result = powmod(a, p, p); if (result != a%p) ++contra; } if (contra == 0) { if (primes.at(p) == false) { uint64_t lastCount = 0; QueryPerformanceCounter((LARGE_INTEGER*)&lastCount); uint64_t delta = lastCount - startCount; cout << "p = " << p; cout << " \ttime: " << 1000 * delta / frequency << " ms" << endl; } } }//for p
Damit klappt es zumindest bis 10^6. Kennt da jemand einen verfeinerten Algorithmus? Mit dem vollständigen Wissen von https://de.wikibooks.org/wiki/Pseudoprimzahlen:_Tabelle_Fermatsche_Pseudoprimzahlen könnte man sicher eine geschwindigkeits-optimierte Erzeugung (ohne liar lookup) schaffen, oder man müsste alle potentiellen Carmichael-Zahlen einem weiteren Test (?) unterziehen.
-
Bei echten Carmichael Zahlen musst du tatsächlich alle teilerfremde a kleiner deiner Zahl einsetzen. Gibt aber bestimmt abküzungen, ich überleg mir da noch mal was.
Glaube die teilerfremden Primzahlen bis zur Carmichael Zeil reiche, hab aber noch keinen Beweis.
-
Hat deine GMPlib kein powmod oder so?
-
@Bengo: Generell irgendwelche Carmichaelzahlen nach dem Muster (6m+1) * (12m+1) * (18m+1) mit allen Faktoren isPrime zu generieren, ist kein Problem ...
#define NOMINMAX // due to conflict with max() #include <windows.h> // warning: no C++ standard! #include <iostream> #include <vector> #include <cmath> #include <cstdint> #include <cassert> using namespace std; uint64_t primeLimit = 1000000; int main() { uint64_t frequency; QueryPerformanceFrequency((LARGE_INTEGER*)&frequency); cout << "cpu frequency: " << frequency << " Hz" << endl; cout << "Carmichael Numbers:" << endl << endl; ////////////////////////// Generate primes lookup table /////////////////////////////////////////////////////////// vector<bool> primes(primeLimit + 1, true); // calculated primes primes.at(0) = primes.at(1) = false; int iMax = (int)sqrt((double)primeLimit) + 1; for (int i = 2; i < iMax; ++i) if (primes.at(i)) for (int j = 2 * i; j < primeLimit + 1; j += i) primes.at(j) = false; // Erastothenes sieve ////////////////////////// Fermat's little theorem /////////////////////////////////////////////////////////////// uint64_t startCount = 0; QueryPerformanceCounter((LARGE_INTEGER*)&startCount); for (uint64_t m = 1; m<10000; ++m) { if (primes.at(m*6+1) && primes.at(m*12+1) && primes.at(m*18+1)) { uint64_t lastCount = 0; QueryPerformanceCounter((LARGE_INTEGER*)&lastCount); uint64_t delta = lastCount - startCount; cout << "m = " << m << "\tc = " << (6*m+1) * (12*m+1) * (18*m+1) << "\ttime: " << 1000 * delta / frequency << " ms" << endl; } } }
cpu frequency: 2435947 Hz Carmichael Numbers: m = 1 c = 1729 time: 0 ms m = 6 c = 294409 time: 0 ms m = 35 c = 56052361 time: 0 ms m = 45 c = 118901521 time: 1 ms m = 51 c = 172947529 time: 1 ms m = 55 c = 216821881 time: 2 ms m = 56 c = 228842209 time: 2 ms m = 100 c = 1299963601 time: 3 ms m = 121 c = 2301745249 time: 3 ms m = 195 c = 9624742921 time: 3 ms m = 206 c = 11346205609 time: 5 ms m = 216 c = 13079177569 time: 5 ms m = 255 c = 21515221081 time: 6 ms m = 276 c = 27278026129 time: 6 ms m = 370 c = 65700513721 time: 7 ms m = 380 c = 71171308081 time: 7 ms m = 426 c = 100264053529 time: 8 ms m = 506 c = 168003672409 time: 8 ms m = 510 c = 172018713961 time: 8 ms m = 511 c = 173032371289 time: 9 ms m = 710 c = 464052305161 time: 9 ms m = 741 c = 527519713969 time: 10 ms m = 800 c = 663805468801 time: 10 ms m = 825 c = 727993807201 time: 11 ms m = 871 c = 856666552249 time: 12 ms m = 930 c = 1042789205881 time: 13 ms m = 975 c = 1201586232601 time: 14 ms m = 1025 c = 1396066334401 time: 15 ms m = 1060 c = 1544001719761 time: 16 ms m = 1115 c = 1797002211241 time: 17 ms m = 1140 c = 1920595706641 time: 17 ms m = 1161 c = 2028691238689 time: 18 ms m = 1270 c = 2655343122121 time: 19 ms m = 1280 c = 2718557844481 time: 20 ms m = 1281 c = 2724933935809 time: 21 ms m = 1311 c = 2920883888089 time: 22 ms m = 1336 c = 3091175755489 time: 23 ms m = 1361 c = 3267961077889 time: 23 ms m = 1365 c = 3296857440241 time: 24 ms m = 1381 c = 3414146271409 time: 25 ms m = 1420 c = 3711619793521 time: 26 ms m = 1421 c = 3719466204049 time: 27 ms m = 1441 c = 3878725359169 time: 28 ms m = 1490 c = 4287981117241 time: 28 ms m = 1515 c = 4507445537641 time: 30 ms m = 1696 c = 6323547512449 time: 30 ms m = 1805 c = 7622722964881 time: 31 ms m = 1875 c = 8544361005001 time: 32 ms m = 1885 c = 8681793690961 time: 33 ms m = 1931 c = 9332984447209 time: 34 ms m = 2040 c = 11004252611041 time: 35 ms m = 2065 c = 11413778221441 time: 36 ms m = 2086 c = 11765530852489 time: 37 ms m = 2191 c = 13633039686169 time: 38 ms m = 2235 c = 14470947115561 time: 38 ms m = 2246 c = 14685655594249 time: 39 ms m = 2256 c = 14882678745409 time: 40 ms m = 2271 c = 15181505298649 time: 41 ms m = 2366 c = 17167430884969 time: 42 ms m = 2425 c = 18483957064801 time: 43 ms m = 2520 c = 20742413217121 time: 43 ms m = 2565 c = 21873528379441 time: 44 ms m = 2571 c = 22027380041449 time: 45 ms m = 2656 c = 24285059687809 time: 46 ms m = 2681 c = 24977268314209 time: 47 ms m = 2711 c = 25825129162489 time: 48 ms m = 2876 c = 30833142247729 time: 49 ms m = 2960 c = 33614369156161 time: 50 ms m = 3020 c = 35700127755121 time: 51 ms m = 3075 c = 37686301288201 time: 52 ms m = 3131 c = 39782913594409 time: 53 ms m = 3341 c = 48336382727569 time: 53 ms m = 3451 c = 53269464581929 time: 54 ms m = 3531 c = 57060521336809 time: 55 ms m = 3566 c = 58774132848169 time: 56 ms m = 3636 c = 62303597046289 time: 57 ms m = 3741 c = 67858397221969 time: 58 ms m = 3796 c = 70895483772049 time: 58 ms m = 3835 c = 73103085605161 time: 59 ms m = 3916 c = 77833567590769 time: 60 ms m = 3971 c = 81159260227849 time: 61 ms m = 3976 c = 81466208375329 time: 62 ms m = 4056 c = 86483161466209 time: 63 ms m = 4140 c = 91968282854641 time: 64 ms m = 4195 c = 95682503446921 time: 65 ms m = 4235 c = 98445661027561 time: 65 ms m = 4340 c = 105950928237841 time: 66 ms m = 4426 c = 112374872517529 time: 67 ms m = 4510 c = 118895125737961 time: 68 ms m = 4511 c = 118974229155289 time: 69 ms m = 4556 c = 122570307044209 time: 70 ms m = 4615 c = 127393969917241 time: 71 ms m = 4636 c = 129140929242289 time: 71 ms m = 4731 c = 137243534644009 time: 72 ms m = 5061 c = 168011973623089 time: 73 ms m = 5155 c = 177548395641481 time: 74 ms m = 5221 c = 184455452572849 time: 75 ms m = 5295 c = 192410140250521 time: 76 ms m = 5326 c = 195809339861929 time: 76 ms m = 5376 c = 201375886537729 time: 77 ms m = 5550 c = 221568419989801 time: 78 ms m = 5571 c = 224093003069449 time: 79 ms m = 5581 c = 225301895806609 time: 80 ms m = 5645 c = 233141908767121 time: 81 ms m = 5791 c = 251703127095769 time: 82 ms m = 5875 c = 262815637149001 time: 83 ms m = 6006 c = 280790932830409 time: 84 ms m = 6226 c = 312790579286329 time: 85 ms m = 6265 c = 318705390188641 time: 86 ms m = 6470 c = 351025246957321 time: 86 ms m = 6650 c = 381144706349401 time: 87 ms m = 6656 c = 382177291511809 time: 88 ms m = 6685 c = 387194417159761 time: 89 ms m = 6706 c = 390854788519609 time: 90 ms m = 6835 c = 413847154073161 time: 91 ms m = 6846 c = 415848433183849 time: 91 ms m = 6915 c = 428549255564041 time: 92 ms m = 6975 c = 439801455648601 time: 93 ms m = 6981 c = 440937387145009 time: 94 ms m = 7066 c = 457240489374169 time: 95 ms m = 7196 c = 482944146230449 time: 96 ms m = 7286 c = 501291932351689 time: 97 ms m = 7331 c = 510637565929609 time: 98 ms m = 7510 c = 548962252005961 time: 99 ms m = 7590 c = 566692953864841 time: 99 ms m = 7616 c = 572536569523969 time: 100 ms m = 7741 c = 601192212565969 time: 101 ms m = 7826 c = 621214363151929 time: 102 ms m = 7860 c = 629346067180561 time: 103 ms m = 7976 c = 657623122439329 time: 104 ms m = 8015 c = 667316922191641 time: 104 ms m = 8081 c = 683938014196609 time: 105 ms m = 8151 c = 701865606427129 time: 106 ms m = 8246 c = 726693182050249 time: 107 ms m = 8305 c = 742403294138881 time: 108 ms m = 8441 c = 779475417411169 time: 109 ms m = 8470 c = 787536877909321 time: 109 ms m = 8500 c = 795934611306001 time: 110 ms m = 8536 c = 806090432846689 time: 111 ms m = 8651 c = 839110734385129 time: 112 ms m = 8680 c = 847577589374881 time: 114 ms m = 8801 c = 883519506462529 time: 115 ms m = 8900 c = 913671191480401 time: 116 ms m = 8921 c = 920153949774049 time: 117 ms m = 8925 c = 921392227198801 time: 117 ms m = 8976 c = 937277770955329 time: 118 ms m = 8980 c = 938531360353681 time: 119 ms m = 8981 c = 938844932257009 time: 120 ms m = 9025 c = 952711345022401 time: 120 ms m = 9046 c = 959377262271049 time: 121 ms m = 9186 c = 1004612946644089 time: 122 ms m = 9355 c = 1061085945064681 time: 123 ms m = 9600 c = 1146654351705601 time: 124 ms m = 9641 c = 1161408537694369 time: 125 ms m = 9676 c = 1174103262876529 time: 125 ms m = 9761 c = 1205317701684289 time: 126 ms m = 9851 c = 1238966116844329 time: 127 ms m = 9890 c = 1253739456971641 time: 128 ms m = 9981 c = 1288666276813009 time: 129 ms
... aber lückenlos alle vollständig und zuverlässig zu finden, das erscheint mir nicht so einfach, zumindest nicht schnell.
-
selfistheman schrieb:
Hat deine GMPlib kein powmod oder so?
Doch, muss ich mal testen. Bei MPIR (GMP) sicher gut optimiert, bei BigUnsigned eher weniger, habe den algo in meinem reduzierten Paket nicht drinnen. BigUnsigned ist leichter zu testen als MPIR (für x64 kompilieren ==> h,lib, dll; Verzeichnisse setzen; lib einbinden; dll zur exe in den Ordner)
-
Ich meinte den Fermat-Test für alle Basen zu machen, die teilerfremd und kleiner als die Carmichael Zahl sind. 2,3,4,5 reicht einfach nicht. Und 4 ist glaube unnötig, da 2^2.
-
Das erfordert aber eine aufwändige Primzahlzerlegung der Carmichael-Zahl? Damit überführt man die Lügner sicher schon.
-
Erhard Henkes schrieb:
Das erfordert aber eine aufwändige Primzahlzerlegung der Carmichael-Zahl? Damit überführt man die Lügner sicher schon.
Naja kannst auch alle Zahlen testen, und wenn statt es nicht 1 bei fermat ergibt, dann guckst du ob es ein Teiler der Carmichael Zahl ist.
Ich glaube nicht, dass man an der Primzahlfaktorzerlegung ablesen kann, ob Carmicheal oder nicht.
Ups da hätte ich mich mal schlau machen müssen: http://www.mathepedia.de/Carmichael-Zahlen.aspx Das Theorem von Korselt kann nämlich genau das.
-
Ich versuch mich mal an einem Beweis:
Wir haben eine Zahl n, von der wir wissen wollen, ob sie eine Carmichael Zahl ist, oder nicht.
Die Bedinung dafür ist, dass für alle
a^{n-1} \equiv 1 \mod nFür a = 1 gilt dies immer, und ist damit uninteressant.
Wir nehmen nun das kleinste und zeigen, dass die Fermat-Bedinung wahr ist. Damit ist sie auch für alle Potenzen dieses a wahr. (Einfach 2 mal die Bedinung und dann multiplizieren).
Allgemein gilt: Wenn die Fermat-Bedinung für d und e gilt, dann gilt sie auch für d*e. Es reicht also tatsächlich aus, alle Primzahlen bis n zu testen, wenn die Bedingung aber nicht erfüllt ist, muss noch sichergestellt werden, dann dieses a kein Teiler von n ist. Wenn das für alle a gilt, dann ist n eine Carmichael Zahl.
-
Ich habe gcd und powm aus GMP umgesetzt.
Alle a ist mir etwas zu viel.
Ich nehme mal 2 bis 19, um Lügner zu unterdrücken (bei 2 bis 13 tauchen die tatsächlich auf) und achte auf den gcd(a,p) == 1 und nehme nur a, die prim sind! Nur dann wird a^p % p = a % p geprüft (a^(p-1)%p = 1 klappt bei mir nicht, habe ich wahrscheinlich falsch verstanden):
#define NOMINMAX // due to conflict with max() #include <windows.h> // warning: no C++ standard! #include <iostream> #include <vector> #include <cmath> #include <cstdint> #include "mpirxx.h" using namespace std; int primeLimit = 1000000; mpz_class powm(mpz_class base, mpz_class exp, mpz_class mod) { mpz_class result = 0; mpz_powm (result.get_mpz_t(), base.get_mpz_t(), exp.get_mpz_t(), mod.get_mpz_t()); return result; } mpz_class gcd(mpz_class bu1, mpz_class bu2) { mpz_class result = 0; mpz_gcd (result.get_mpz_t(), bu1.get_mpz_t(), bu2.get_mpz_t()); return result; } int main() { uint64_t frequency; QueryPerformanceFrequency((LARGE_INTEGER*)&frequency); cout << "cpu frequency: " << frequency << " Hz" << endl; cout << "Carmichael Numbers:" << endl << endl; ////////////////////////// Generate primes lookup table /////////////////////////////////////////////////////////// vector<bool> primes(primeLimit + 1, true); // calculated primes primes.at(0) = primes.at(1) = false; int iMax = (int)sqrt((double)primeLimit) + 1; for (int i = 2; i < iMax; ++i) if (primes.at(i)) for (int j = 2 * i; j < primeLimit + 1; j += i) primes.at(j) = false; // Erastothenes sieve ////////////////////////// Fermat's little theorem /////////////////////////////////////////////////////////////// uint64_t startCount = 0; QueryPerformanceCounter((LARGE_INTEGER*)&startCount); for (int p = 2; p <= primeLimit; ++p) { uint64_t contra = 0; for (mpz_class a = 2; a <= 19; ++a) { mpz_class d = gcd(a, p); if (d != 1 || ((primes.at(a.get_ui())) == false)) continue; if (powm(a, p - 1, p) != 1) ++contra; }//for a if (contra == 0) { if (primes.at(p) == false) { uint64_t lastCount = 0; QueryPerformanceCounter((LARGE_INTEGER*)&lastCount); uint64_t delta = lastCount - startCount; cout << "p = " << p; cout << " \ttime: " << 1000 * delta / frequency << " ms" << endl; } } }//for p }
cpu frequency: 3127001 Hz Carmichael Numbers: p = 561 time: 10 ms p = 1105 time: 22 ms p = 1729 time: 36 ms p = 2465 time: 52 ms p = 2821 time: 61 ms p = 6601 time: 144 ms p = 8911 time: 199 ms p = 10585 time: 239 ms p = 15841 time: 359 ms p = 29341 time: 706 ms p = 41041 time: 990 ms p = 46657 time: 1130 ms p = 52633 time: 1276 ms p = 62745 time: 1524 ms p = 63973 time: 1556 ms p = 75361 time: 1838 ms p = 101101 time: 2493 ms p = 115921 time: 2863 ms p = 126217 time: 3123 ms p = 162401 time: 4047 ms p = 172081 time: 4302 ms p = 188461 time: 4729 ms p = 252601 time: 6368 ms p = 278545 time: 7033 ms p = 294409 time: 7445 ms p = 314821 time: 7979 ms p = 334153 time: 8496 ms p = 340561 time: 8670 ms p = 399001 time: 10215 ms p = 410041 time: 10501 ms p = 449065 time: 11519 ms p = 488881 time: 12556 ms p = 512461 time: 13171 ms p = 530881 time: 13654 ms p = 552721 time: 14226 ms p = 656601 time: 16989 ms p = 658801 time: 17048 ms p = 670033 time: 17344 ms p = 748657 time: 19445 ms p = 825265 time: 21484 ms p = 838201 time: 21820 ms p = 852841 time: 22201 ms p = 997633 time: 25997 ms
Der Ausdruck im Vergleich mit https://de.wikibooks.org/wiki/Pseudoprimzahlen:_Tabelle_Carmichael-Zahlen passt exakt bis 10^6!
Das Programm läuft wirklich schnell durch. Echt klasse. Danke Bengo!
-
Die Vorstufe des Fermats den du nennst a^p = a mod p. Gild für alle a.
Der eigentliche Fermat a^(p-1) = 1 mod p. (Die Vorstufe durch a geteilt) geht nur, wenn a und p teilerfremd sind, also gcd(a,p) = 1.Und ja es macht nur Sinn prime a zu prüfen, da wie ich bereits schrieb, wenn es für die Basen a und b gilt, dann auch für die Base a*b.
-
Könntest du bitte die high_resolution_clock verwenden? Ich sehe keinen Grund, wieso du windows.h brauchst.
Wieso brauchst du eigentlich GMP? Du könntest alles mit ints machen, das wäre um ein Vielfachesn schneller.
Und bitte mach das at() weg, das sieht ja hässlich aus.
-
Der eigentliche Fermat a^(p-1) = 1 mod p. (Die Vorstufe durch a geteilt) geht nur, wenn a und p teilerfremd sind, also gcd(a,p) = 1.
gcd(a,p) == 1 prüfen wir in der for-loop. Also geht auch das einfachere:
if (powm(a, p-1, p) != 1)
Habe es oben ausgebessert. Bringt nochmal 1 sec bis 1 Mio.
Könntest du bitte die high_resolution_clock verwenden? Ich sehe keinen Grund, wieso du windows.h brauchst.
Ja, kann man machen. Mache ich bei der nächsten Änderung.
Wieso brauchst du eigentlich GMP? Du könntest alles mit ints machen, das wäre um ein Vielfachesn schneller.
Wir rechnen hier (19 ^ 1000000) % 19. Da wäre uint64_t schnell out_of_range.
Und bitte mach das at() weg, das sieht ja hässlich aus.
Stimmt, das ist hässliches C++, dafür wird der Bereich auf Gültigkeit geprüft. Safety first!
Sehe ich das richtig, dass ich mir so etwas wie
mpz_class powm(mpz_class base, mpz_class exp, mpz_class mod)
selbst schreiben muss? Habe dies in mpirxx.h nicht gefunden.
-
Erhard Henkes schrieb:
Wir rechnen hier (19 ^ 1000000) % 19. Da wäre uint64_t schnell out_of_range.
Quatsch. Das ist doch der Trick an PowMod
-
Zuerst mal vielen Dank an alle, die hier bei diesem staubtrockenen Thema konstruktiv mitmischen! Keine Ahnung, warum mich diese Sachen als Nicht-Mathematiker faszinieren.
-
volkard schrieb:
Erhard Henkes schrieb:
Wir rechnen hier (19 ^ 1000000) % 19. Da wäre uint64_t schnell out_of_range.
Quatsch. Das ist doch der Trick an PowMod
Ja, Du hast recht, aber das powm von GMP ist echt schnell. Kannst es ja mal mit deinem vergleichen.
-
Erhard Henkes schrieb:
volkard schrieb:
Erhard Henkes schrieb:
Wir rechnen hier (19 ^ 1000000) % 19. Da wäre uint64_t schnell out_of_range.
Quatsch. Das ist doch der Trick an PowMod
Ja, Du hast recht, aber das powm von GMP ist echt schnell. Kannst es ja mal mit deinem vergleichen.
Oder besser...
997633 passt ja auch gut tausendmal in einen uint32_t.
Falls Du noch kein schnelles mulmod64 hast, ein schnelles mulmod32 wäre return uint64(a)*b;. Könnte mir vorstellen, daß es 100-mal so schnell davon wird.
-
Mit a^1000000 mod b = a^(1000000 mod b) mod b und b<=19 passen alle Zwischenresultate in einen 64-Bit-Integer und du musst nur einmal am Schluss modulo rechnen.
-
Erhard Henkes schrieb:
Zuerst mal vielen Dank an alle, die hier bei diesem staubtrockenen Thema konstruktiv mitmischen! Keine Ahnung, warum mich diese Sachen als Nicht-Mathematiker faszinieren.
Quatsch Zahlentheorie ist gleich nach der Algebra das schönste was es gibt