Carmichael Zahlen (Visualisierung und lückenlose Erzeugung)
-
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
-
Ich habe es mit diesem Programm (auf Wunsch mit <chrono> ) bis 10^7 probiert: http://codepad.org/oK7tLmv6
... und schon wieder zwei neue "Lügner": 5968873 und 9699690
a = 2 ... 19 reicht da also schon wieder nicht. Für diesen Zahlenbereich muss man den range von a erweitern. Bei a = 2 ... 29 verschwindet 9699690, aber die 5968873 bleibt. Ich denke, die hat man auf https://de.wikibooks.org/wiki/Pseudoprimzahlen:_Tabelle_Carmichael-Zahlen schlicht vergessen (43*127*1093 = 5968873). https://de.numberworld.info/5968873
Hier steht sie dabei: http://oeis.org/A097061/internal
Nobody is perfect. Ich ergänze die dort mal in diesem wikibook.
Endlich mal ein nützliches Resultat.9699690 = 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 <== finde ich auch bewundernswert. hat aber keinen Namen, soweit ich das weiß.
EDIT: "The p-primorial is the product of all primes less than or equal to p. It is sometimes denoted by p#. ... First ten: 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870, 6469693230" ( http://www.numbergossip.com/list )
-
@Bengo: Habe ich das so richtig verstanden bei Dir?
a=2 ... p ?
for (mpz_class a = 2; a <= p; ++a)
Da dauert der Programmlauf ewig. Völlig indiskutabel.
Carmichael Numbers: 561 time: 122 ms 1105 time: 474 ms 1729 time: 1162 ms 2465 time: 2359 ms 2821 time: 3100 ms 6601 time: 16911 ms 8911 time: 30719 ms 10585 time: 43505 ms 15841 time: 97253 ms
-
Die Primzahlen von 2 bis zur Carmichael-Zahl reichen, wie bereits gesagt (und wenn nicht 1 rauskommt, dann prüfen ob die Zahl Teiler von Carmichael-Zahl ist, die Sache mit dem ggT kann man sich dann auch komplett sparen). Wenn du es nicht machst, kannst du Lügner aber nie ausschließen
-
Irgendetwas machst du grundsätzlich falsch...
https://ideone.com/eZX8Kr
-
@selfistheman: Danke für den Nachweis.
-
int gcd(int a, int b) { return !b ? a : gcd(b, a%b); }
Sowas versaut einem doch den ganzen Tag.
-
Du brauchst nicht mal die gdc-Funktion (Ich hab es schon 2 bis 3 mal hier geschrieben, du musst nicht vorher auf teilerfremd prüfen, es reicht wenn nachher nicht 1 rauskommt, zu prüfen ob die Zahl teilerfremd ist. Wenn du nur Primzahlen testest, ist das äquivalent dazu, dass die Primzahl ein Teiler ist), und wenn du ihn doch brauchen wilst, solltest du den euklidischen Algo auch richtig implementieren.
Ich geb dir mal ein Zahlenbeispiel:
117 113:
117 = 1 * 113 + 4
113 = 28 * 4 + 1
4 = 4 * 1 + 0
Er endet wenn 0 rauskommt, und der Rest der vorletzen Zeile ist der ggT. Du bruachst also nicht nur den Rest, sondern auch den ganzahligen Anteil der Zahl.