Carmichael Zahlen (Visualisierung und lückenlose Erzeugung)



  • 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{aZggT(a,n)=1}a \in \{ a \in Z | ggT(a,n) = 1 \}
    a^{n-1} \equiv 1 \mod n

    Für a = 1 gilt dies immer, und ist damit uninteressant.
    Wir nehmen nun das kleinste a1a \neq 1 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.


Anmelden zum Antworten