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


Anmelden zum Antworten