Carmichael Zahlen (Visualisierung und lückenlose Erzeugung)

Erhard Henkes

@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.

Erhard Henkes

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)

Bengo

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.

Erhard Henkes

Das erfordert aber eine aufwändige Primzahlzerlegung der Carmichael-Zahl? Damit überführt man die Lügner sicher schon.

Bengo

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.

Bengo

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 \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 $a \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.

Erhard Henkes

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!

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.

Erhard Henkes

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.

volkard

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

Erhard Henkes

Zuerst mal vielen Dank an alle, die hier bei diesem staubtrockenen Thema konstruktiv mitmischen! Keine Ahnung, warum mich diese Sachen als Nicht-Mathematiker faszinieren.

Erhard Henkes

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.

volkard

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.

Bengo

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

Erhard Henkes

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 )

Erhard Henkes

@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

Bengo

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