Perfekte, superperfekte Zahlen und Mersenne Primzahlen ermitteln (Miller-Rabin, Lucas-Lehmer)



  • Beim Durchgang durch die natürlichen Zahlen als kleines Mathe-Programm in C++ die Ermittlung der perfekten (vollkommenen), der superperfekten Zahlen, incl. p und bei den superperfekten der Mersenne-Primzahl.

    Hintergrund:
    https://de.wikipedia.org/wiki/Vollkommene_Zahl
    https://de.wikipedia.org/wiki/Superperfekte_Zahl

    #define NOMINMAX     // due to conflict with max()
    #include <windows.h> // warning: no C++ standard!
    
    #include <iostream>
    #include <cstdint>
    #include <chrono>
    
    using namespace std;
    
    void setConsole() // windows.h
     {
         _CONSOLE_SCREEN_BUFFER_INFO info;
         HANDLE handle = GetStdHandle(STD_OUTPUT_HANDLE);
         GetConsoleScreenBufferInfo(handle, &info);
         COORD c;
         c.X = std::max<SHORT>(info.dwSize.X, 150);
         c.Y = std::max<SHORT>(info.dwSize.Y, 1000);
         SetConsoleScreenBufferSize(handle, c);
    
         SMALL_RECT RECT;
         RECT.Left = 0;
         RECT.Top = 0;
         RECT.Bottom = std::max(info.srWindow.Bottom - info.srWindow.Top, 40 - 1);
         RECT.Right = std::max(c.X - 1, 100 - 1);
         SetConsoleWindowInfo(handle, TRUE, &RECT);
     }
    
     void textcolor(unsigned short color = 15)  // windows.h
     {
         SetConsoleTextAttribute(::GetStdHandle(STD_OUTPUT_HANDLE), color);
     }
    
    uint64_t getSumOfDivisors(uint64_t n)
    {
        uint64_t sum = 0;
    
        for (uint64_t i=1; i<=n; ++i)
        {
            if (n%i == 0)
            {
                sum += i;
            }
        }
    
        return sum;
    }
    
    //
    uint64_t getP(uint64_t n)
    {
        uint64_t p=1;
    
        while (n>1)
        {
            n>>=1;
            ++p;
        }
    
        return p;
    }
    
    main()
    {
        setConsole();    
        chrono::time_point<chrono::system_clock> start, last;
        start = chrono::system_clock::now();
    
        for (uint64_t n=2; n<=10000000; ++n)
        {
            uint64_t sum1 = getSumOfDivisors(n);
            uint64_t sum2 = getSumOfDivisors(sum1);
    
            if (sum1 == 2 * n)
            {
                textcolor(14);
                cout << "perfect number:       " << n << " \t                ";
                uint64_t p = getP(n)/2 + 1;
                cout << "\tp: " << p;
    
                last = chrono::system_clock::now();
                uint64_t elapsedTime = chrono::duration_cast<chrono::seconds>(last-start).count();
                cout << "\ttime: " << elapsedTime << " s" <<  endl;
            }
    
            if (sum2 == 2 * n)
            {
                textcolor(10);
                cout << "super perfect number: " << n << " \tmersenne prime: " << sum2-1 << " ";
                uint64_t p = getP(n);
                cout << "\tp: " << p;
    
                last = chrono::system_clock::now();
                uint64_t elapsedTime = chrono::duration_cast<chrono::seconds>(last-start).count();
                cout << "\ttime: " << elapsedTime << " s" <<  endl;
            }
        }
    }
    

    Ausgabe:

    super perfect number: 2         mersenne prime: 3       p: 2    time: 0 s
    super perfect number: 4         mersenne prime: 7       p: 3    time: 0 s
    perfect number:       6                                 p: 2    time: 0 s
    super perfect number: 16        mersenne prime: 31      p: 5    time: 0 s
    perfect number:       28                                p: 3    time: 0 s
    super perfect number: 64        mersenne prime: 127     p: 7    time: 0 s
    perfect number:       496                               p: 5    time: 0 s
    super perfect number: 4096      mersenne prime: 8191    p: 13   time: 0 s
    perfect number:       8128                              p: 7    time: 0 s
    super perfect number: 65536     mersenne prime: 131071  p: 17   time: 50 s
    super perfect number: 262144    mersenne prime: 524287  p: 19   time: 805 s
    

    Das Programm geht natürlich schnell in die Knie. Frage an die math geeks: Welche Daten kann man hier sinnvoll in lookup tables schieben, damit man rasch an die Big Unsigned heran kann? 😉

    Man könnte das Prg auch allgemein auf m-superperfect numbers erweitern.



  • Erhard Henkes schrieb:

    Welche Daten kann man hier sinnvoll in lookup tables schieben, damit man rasch an die Big Unsigned heran kann? 😉

    Weiß ich nicht, aber es hilft schon mal viel, bei getSumOfDivisors nur bis zur Wurzel zu laufen und die Teiler als Paare zu summieren:

    uint64_t getSumOfDivisors(uint64_t n)
    {
        //assert(n>=2);
    
        uint64_t sum = 1+n; 
        uint64_t i;
    
        for (i=2; i*i<n; ++i) 
            if (n%i == 0) 
                sum += i + n/i;
    
        if (i*i==n)
            sum += i;
    
        return sum;
    }
    


  • Ja, das beschleunigt signifikant! 🙂

    super perfect number: 2         mersenne prime: 3       p: 2    time: 0 s
    super perfect number: 4         mersenne prime: 7       p: 3    time: 0 s
    perfect number:       6                                 p: 2    time: 0 s
    super perfect number: 16        mersenne prime: 31      p: 5    time: 0 s
    perfect number:       28                                p: 3    time: 0 s
    super perfect number: 64        mersenne prime: 127     p: 7    time: 0 s
    perfect number:       496                               p: 5    time: 0 s
    super perfect number: 4096      mersenne prime: 8191    p: 13   time: 0 s
    perfect number:       8128                              p: 7    time: 0 s
    super perfect number: 65536     mersenne prime: 131071  p: 17   time: 0 s
    super perfect number: 262144    mersenne prime: 524287  p: 19   time: 2 s
    

    Echt klasse, aber das reicht noch lange nicht. Wir brauchen den super Turbo.
    "There are no odd superperfect numbers below 7 * 10^24"
    Multithreading haben wir noch in der Hinterhand.



  • Dein Hauptproblem ist ja die Funktion getSumOfDivisors möglichst schnell hinzukriegen. Ich vermute der beste Ansatz ist über die Primfaktorzerlegung, denn aus dieser lässt sich die gesuchte Summe einfach berechnen (siehe hier). Jetzt steht man allerdings vor dem Problem die Primfaktorzerlegung zu finden, was auch alles andere als einfach ist. Allerdings bin ich der Meinung, dass dies einfacher ist als alle Teiler einer Zahl zu suchen. Besonders wenn die zu prüfende Zahl viele kleine Primfaktoren besitzt, sind diese schnell zu finden und man kann durch diese Teilen und arbeitet mit einer deutlich kleineren Zahl weiter. Außerdem düfte Google so einiges hergeben zu faktorisierungs Algorithmen.



  • Off Topic, aber wenn es schon entfernt um Optimierungen geht:

    divide schrieb:

    ...
        for (i=2; i*i<n; ++i)
    ...
    

    Ich weiss dass sqrt() ein wenig verpönt ist, aber könnte eine Wurzel zu ziehen bei hinreichend großen n nicht effizienter sein, als bei jedem Schleifendurchlauf eine Integer-Multiplikation zu machen?



  • n/i > i
    hast ja unten eh n%i

    denke, dass das der compiler rafft



  • Ich fange mal mit der ersten klugen Optimierung an: Primzahlzerlegung.

    uint64_t gesodi(uint64_t n) {
      uint64_t res = 1;
      while (n != 1) {
        for (uint64_t i = 2; i*i <= n; ++i)
          if (n%i == 0) {
            uint64_t f = 1 + i;
            n /= i;
            while (n%i == 0) {
              n /= i;
              f = 1 + i*f;
            }
            res *= f;
            goto next_loop;
          }
        res *= n + 1;
        break;
      next_loop:;
      }
      return res;
    }
    

    Mit den Standard-Primzahlsiebtechniken lässt sich das natürlich noch weiter optimieren.



  • Die getSumOfDivisors(n) geht so:
    Primfaktorenpotenzen auflisten.
    Dazu bei kleinen Zahlen bis 2^64 gerne PollardRho benutzen, später wasauchimmer die BigNum-Bibliothek hat.
    Mit einer rekursiven Funktion den Teilerverband aufspannen und summieren.
    Fertig.



  • Durch die gesteigerte Geschwindigkeit überwindet man die große Lücke im unteren Millionenbereich und findet eine neue "perfect number":

    calculation process at: 33500000
    perfect number:       33550336                          p: 13   time: 284 s
    
    calculation process at: 34000000
    

    Vgl.: https://en.wikipedia.org/wiki/List_of_perfect_numbers
    Ab dann muss man allerdings bis in den Milliardenbereich ausharren.

    Ich habe versucht, das i*i<=n durch ein i < iMax mit wurzelgezogener Konstante zu ersetzen. Zu meinem Erstaunen war das jedoch langsamer (vlt. arbeitet hier der Compiler intelligent).



  • Erhard Henkes schrieb:

    Welche Daten kann man hier sinnvoll in lookup tables schieben, damit man rasch an die Big Unsigned heran kann? 😉

    Ich überlege wiedermal für Suche bis 2^128, vorerst nicht weiter, weils eh nicht erschöpfend geht.
    Da es wie immer darum geht, Bereiche abzuscannen, nimm Dir einen Bereich von zu untersuchenden Zahlen, a<=n<b, und hau ganz gemein fast wie im Eratosthenes drauf. Also nimmst alle Primzahlen von 2 bis 2^32 (die kriegste in ner Minute mit dem kleinen Erathosthenes aus dem anderen Thread) und "streichst" die Zahlen im hohen Bereich fast wie üblich. Brauchst fürs Streichen nur eine Division um den Startwert zu bestimmen und dann gehts ja mit streichindex<b;streichindex+=smallprime) herauf.
    Nur haste kein array<bool>, sondern ein array<UInt32> als Streichfeld, und Du schreibst nicht true, sondern smallprime da rein.
    Und Schwupps, haste für alle Zahlen a<=n<b ein nettes Vorwissen: Ist array[n]!=0, dann ist array[n] der größte Primteiler unter 2^32 von n. Den kannste erstmal abspalten und dann normal weitermachen.



  • Du kannst auch eine schwache abschätzung der Wurzel nach oben verwenden. Wenn man weniger als die Hälfte der Bits wegnimmt, sollte das größer als die Wurzel sein.



  • Echt jetzt? Keiner bringt die offensichtliche Optimierung zu meiner Funktion?

    uint64_t gesodi(uint64_t n) {
      uint64_t res=1;
      auto update = [&](int p){
        uint64_t f = 1 + p;
        n /= p;
        while (n%p == 0)
          n /= p, f = 1 + p*f;
        res *= f;
      };
      if (n%2==0) update(2);
      if (n%3==0) update(3);
      if (n%5==0) update(5);
      for (uint64_t i=7; i*i<=n;) {
        if (n%i == 0) update(i);
        i += 4;
        if (n%i == 0) update(i);
        i += 2;
      }
      if (n != 1)
        res *= n+1;
      return res;
    }
    


  • Erhard Henkes schrieb:

    Ich habe versucht, das i*i<=n durch ein i < iMax mit wurzelgezogener Konstante zu ersetzen. Zu meinem Erstaunen war das jedoch langsamer (vlt. arbeitet hier der Compiler intelligent).

    Nee.
    Aber EINMAL sqrt(n) ist schon ein wenig teurer als EINMAL ii. Und vielleicht auch als ZWEIMAL oder DREIMAL. Natürlich viel viel billiger als n-MAL.
    Nur isse so, daß 50% der Zahlen schon bei i=2 rausfliegen und 16% bei i=3 und desdawegen kommt verdammt oft EINMAL oder ZWEIMAL vor.
    Aslo hast schon recht, sqrt auszulagern ist das logischte von der Welt. Nur zufällig isses seit ich es beobachte (ca 15 Jahre) immer ein wenig langsamer als i
    i.
    Gut wäre es vielleicht, mit i*i die kleinen Teiler zu machen und so ab 17 oder wo(??) umzuschwenken auf sqrt, weils dann eh vermutlich eine Primzahl ist.
    Optimös wäre es, eh, die Teilersuche mit i=2, i=3, i=5… zu hardcoden, denn da braucht der Compiler nur eine Plutimikation dafür statt einer teuren Division.
    Also raus aus der dicken Funktion damit und findFirstDivisor(n) und dadrin kannste lecker if-return, if-return, if-return, for… machen.


  • Mod

    Hat jemand was von isqrt erwähnt? 🤡 Folgendes scheint nach etwas Probiererei zu funktionieren:

    using u128 = unsigned __int128;
    
    unsigned clz(u128 u) {
    	auto a = __builtin_clzll(u), b = __builtin_clzll(u >> 64);
    	if (b == 64)
    		return a + 64;
    	return b;
    }
    
    u128 isqrt( u128 n )
    {
    	u128 x = n >> (128-clz(n))/2, y;
    	for(;;)
    	{
    		y = x;
    		x += n/x;
    		auto rest = x%2;
    		x /= 2;
    		std::cout << (int)x << ' ' << (int)y << '\n';
    		if( x == y || (y == x+1 && rest) )
    			break;
    	}
    	return y*y == n? y : 0;
    }
    


  • @gesodi: echt stark!

    calculation process at: 33500000
    perfect number:       33550336                          p: 13   time: 87 s
    


  • Aktueller Stand:

    #define NOMINMAX     // due to conflict with max()
    #include <windows.h> // warning: no C++ standard!
    
    #include <iostream>
    #include <cstdint>
    #include <cstdlib>
    #include <algorithm>
    #include <chrono>
    
    using namespace std;
    
    void wait()
    {
    	cout << "Press any key to continue." << endl;
    	cin.clear();
    	cin.ignore(numeric_limits<streamsize>::max(), '\n');
    	cin.get();
    }
    
    void setConsole() // windows.h
    {
    	_CONSOLE_SCREEN_BUFFER_INFO info;
    	HANDLE handle = GetStdHandle(STD_OUTPUT_HANDLE);
    	GetConsoleScreenBufferInfo(handle, &info);
    	COORD c;
    	c.X = max<SHORT>(info.dwSize.X, 150);
    	c.Y = max<SHORT>(info.dwSize.Y, 1000);
    	SetConsoleScreenBufferSize(handle, c);
    
    	SMALL_RECT RECT;
    	RECT.Left = 0;
    	RECT.Top = 0;
    	RECT.Bottom = max(info.srWindow.Bottom - info.srWindow.Top, 40 - 1);
    	RECT.Right = max(c.X - 1, 100 - 1);
    	SetConsoleWindowInfo(handle, TRUE, &RECT);
    }
    
    void textcolor(unsigned short color = 15)  // windows.h
    {
    	SetConsoleTextAttribute(::GetStdHandle(STD_OUTPUT_HANDLE), color);
    }
    
    uint64_t getSumOfDivisors(uint64_t n)
    {
    	uint64_t res = 1;
    
    	auto update = [&](uint64_t p) 
    	{
    		uint64_t f = 1 + p;
    		n /= p;
    		while (n%p == 0)
    			n /= p, f = 1 + p*f;
    		res *= f;
    	};
    
    	if (n % 2 == 0) update(2);
    	if (n % 3 == 0) update(3);
    	if (n % 5 == 0) update(5);
    
    	for (uint64_t i = 7; i*i <= n;) 
    	{
    		if (n%i == 0) update(i);
    		i += 4;
    		if (n%i == 0) update(i);
    		i += 2;
    	}
    
    	if (n != 1)
    		res *= n + 1;
    
    	return res;
    }
    
    uint64_t getP(uint64_t n)
    {
    	uint64_t p = 1;
    
    	while (n>1)
    	{
    		n >>= 1;
    		++p;
    	}
    
    	return p;
    }
    
    bool checkForPerfectNumber(uint64_t n, uint64_t sum, chrono::time_point<chrono::system_clock>& start)
    {
    	if (sum == 2 * n)
    	{
    		textcolor(14);
    		cout << "perfect number:       " << n << " \t                ";
    		uint64_t p = getP(n) / 2 + 1;
    		cout << "\tp: " << p;
    		uint64_t elapsedTime = chrono::duration_cast<chrono::seconds>(chrono::system_clock::now() - start).count();
    		cout << "\ttime: " << elapsedTime << " s" << endl;
    		return true;
    	}
    	return false;
    }
    
    bool checkForSuperPerfectNumber(uint64_t n, uint64_t sum, chrono::time_point<chrono::system_clock>& start)
    {
    	if (sum == 2 * n)
    	{
    		textcolor(10);
    		cout << "super perfect number: " << n << " \tmersenne prime: " << sum - 1 << " ";
    		uint64_t p = getP(n);
    		cout << "\tp: " << p;
    		uint64_t elapsedTime = chrono::duration_cast<chrono::seconds>(chrono::system_clock::now() - start).count();
    		cout << "\ttime: " << elapsedTime << " s" << endl;
    		return true;
    	}
    	return false;
    }
    
    void checkNumber(uint64_t n, chrono::time_point<chrono::system_clock>& start)
    {	
    	uint64_t sum1 = getSumOfDivisors(n);
    	uint64_t sum2 = getSumOfDivisors(sum1);
    
    	checkForPerfectNumber(n, sum1, std::ref(start));
    	checkForSuperPerfectNumber(n, sum2, std::ref(start));	
    }
    
    int main()
    {
    	setConsole();
    	chrono::time_point<chrono::system_clock> start, last;
    	start = chrono::system_clock::now();
    	chrono::time_point<chrono::system_clock>& refStart = start;
    
    	const uint64_t lastNumber = 10000000000;
    
    	for (uint64_t n = 2; n <= lastNumber; ++n)
    	{		
    		if (n % 1000000 == 0)
    		{
    			textcolor(15);
    			cout << "calc. at: " << n << endl;			
    		}
    
    		checkNumber(n, std::ref(refStart));				
    	}
    
    	wait();
    }
    


  • Hat man eine Mersenne-Primzahl, also 2^n - 1 (Mersenne-Zahl) und Primzahl, so kann man daraus mittels n perfekte und superperfekte Zahlen ableiten:
    Perfekte Zahl: 2(n-1)*(2n -1) ==> mit Summe der Teiler(n) ist gleich n
    Superperfekte Zahl: 2^(n-1) ==> Summe der Teiler (Summe der Teiler(n)) ist gleich 2*n

    Daher habe ich versucht ausgehend von den Mersenne-Zahlen 2^n - 1 mittels Divisions-Precheck (mit Primzahlen von 3 bis 3 Mio.) und anschließendem iterativem Miller-Rabin mit verschiedenen zufälligen Basen als "Zeuge" die Mersenne-Primzahlen zu bestimmen. Das gelingt im Bereich n: 2 ... 2300 zuverlässig (mathematisch und bez. multithreading noch nicht optimiert), ist aber für hohe Potenzen n einfach noch nicht schnell genug:

    http://pastebin.com/ENNCpcWV <== Sourcecode

    Resultat:

    p: 2      time: 9 ms
    p: 3      time: 0 ms
    p: 5      time: 0 ms
    p: 7      time: 0 ms
    p: 13     time: 1 ms
    p: 17     time: 0 ms
    p: 19     time: 0 ms
    p: 31     time: 3 ms
    p: 61     time: 16 ms
    p: 89     time: 49 ms
    p: 107    time: 66 ms
    p: 127    time: 93 ms
    p: 521    time: 26615 ms
    p: 607    time: 16048 ms
    p: 1279   time: 599282 ms
    p: 2203   time: 3938991 ms
    p: 2281   time: 825088 ms
    total time: 6296721 ms
    

    Komplett: http://pastebin.com/7tNMn36k

    Eine Mersenne-Primzahl mit p=2281 ergibt mit 2^p - 1 = 4,460875571837584295711517064021 * 10^686

    Das sind gewaltige Zahlenbereiche. Welche Methodik verwendet man hier in der Top-Liga, und welche Chancen hat der Hobbyist mit seinem PC?



  • Beispiel: http://www.math.rwth-aachen.de/~Gabriele.Nebe/Vorl/SEMCA/Mersenne.pdf
    https://de.wikipedia.org/wiki/Lucas-Lehmer-Test
    https://de.wikiversity.org/wiki/Benutzer:Tanik/Lucas-Lehmer-Test_für_Mersennezahlen/Beispiel1

    Mit dem (verbesserten) Lucas-Lehmer-Test:

    #define NOMINMAX     // due to conflict with max()
    #include <windows.h> // warning: no C++ standard!
    
    #include "MyBigUnsigned.h"
    
    #include <iostream>
    #include <cstdint>
    
    using namespace std;
    
    void textcolor(unsigned short color = 15)
    {
    	SetConsoleTextAttribute(::GetStdHandle(STD_OUTPUT_HANDLE), color);
    }
    
    BigUnsigned calculateMersenneNumber(BigUnsigned p)
    {
    	BigUnsigned x = 2;
    
    	for (BigUnsigned i = 1; i < p; ++i)
    	{
    		x <<= 1;
    	}
    	return (x - 1);
    }
    
    bool LucasLehmerTest(BigUnsigned p)
    {
    	BigUnsigned m = calculateMersenneNumber(p), s = 4;
    
    	for (BigUnsigned i = 2; i < p; ++i)
    		s = (s*s - 2) % m;
    
    	return (s == 0);
    }
    
    int main()
    {
    	uint64_t startCount = 0, lastCount = 0, delta = 0, frequency;
    
    	QueryPerformanceFrequency((LARGE_INTEGER*)&frequency);
    	cout << "cpu frequency: " << frequency << " Hz" << endl << endl;
    
    	BigUnsigned p;
    
    	cout << "Lucas-Lehmer-Test (Mersenne-Zahlen)" << endl; // https://de.wikipedia.org/wiki/Lucas-Lehmer-Test
    
    	QueryPerformanceCounter((LARGE_INTEGER*)&startCount);
    
    	for (BigUnsigned p = 0; p < 10000; ++p)
    	{
    		if (LucasLehmerTest(p))
    		{			
    			QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
    			uint64_t delta = lastCount - startCount;
    
    			textcolor(15);
    			cout << "p: " << p;
    			textcolor(2);
    			cout << "\ttime: " << 1000 * delta / frequency << " ms" << endl;
    		}
    		else
    		{
    			//cout << " ist keine Mersenne-Primzahl" << endl;
    		}
    	}	
    }
    
    cpu frequency: 3127021 Hz
    
    Lucas-Lehmer-Test (Mersenne-Zahlen)
    p: 3    time: 0 ms
    p: 5    time: 0 ms
    p: 7    time: 0 ms
    p: 13   time: 1 ms
    p: 17   time: 1 ms
    p: 19   time: 2 ms
    p: 31   time: 2 ms
    p: 61   time: 6 ms
    p: 89   time: 14 ms
    p: 107  time: 23 ms
    p: 127  time: 39 ms
    p: 521  time: 6140 ms
    p: 607  time: 10858 ms
    p: 1279 time: 191493 ms
    p: 2203 time: 1617345 ms
    p: 2281 time: 1855752 ms
    p: 3217 time: 7159214 ms
    p: 4253 time: 21421841 ms
    p: 4423 time: 25019798 ms
    

    Die Ergebnisse sind korrekt: http://www.mersenne.org/primes/

    Der Lucas-Lehmer-Test ist zeitlich besser als der iterierende Miller-Rabin-Test mit vorgelagertem Divisions-Precheck, aber für sehr hohe Zahlenbereiche, wie bei Mersenne-Primzahlen aktuell notwendig, leider immer noch lahm.

    Probleme: BigUnsigned, Multiplikation (==> Multiplikation mit schneller Fourier-Transformation, Schöhage und Strassen, 1971), noch single threaded

    Eine deutlich spürbare Beschleunigung erzielt man durch einen Divisions-Precheck 3 ... 293 (dort liegt in etwa das Optimum) ab p = 1000:

    #define NOMINMAX     // due to conflict with max()
    #include <windows.h> // warning: no C++ standard!
    
    #include "MyBigUnsigned.h"
    
    #include <iostream>
    #include <cstdint>
    #include <vector>
    #include <cassert>
    
    using namespace std;
    
    const uint32_t DIVISOR_PRECHECK_MAX = 293;
    
    void textcolor(unsigned short color = 15)
    {
    	SetConsoleTextAttribute(::GetStdHandle(STD_OUTPUT_HANDLE), color);
    }
    
    uint64_t convertBigToU64(BigUnsigned num)
    {
    	BigUnsigned ZweiHochVierUndSechzig = stringToBigUnsigned("18446744073709551616");
    	assert(num < ZweiHochVierUndSechzig);
    	uint64_t jHi = num.getBlock(1);
    	uint64_t jLo = num.getBlock(0);
    	return (jHi << 32 | jLo);
    }
    
    void convertU64ToBig(uint64_t num, BigUnsigned& b)
    {
    	uint32_t numLo = (uint32_t)(num & 0xFFFFFFFF);
    	uint32_t numHi = (uint32_t)((num >> 32) & 0xFFFFFFFF);
    	b.setBlock(0, numLo);
    	b.setBlock(1, numHi);
    }
    
    BigUnsigned calculateMersenneNumber(BigUnsigned p)
    {
    	BigUnsigned x = 2;
    
    	for (BigUnsigned i = 1; i < p; ++i)
    	{
    		x <<= 1;
    	}
    	return (x - 1);
    }
    
    bool LucasLehmerTest(BigUnsigned p)
    {
    	BigUnsigned m = calculateMersenneNumber(p), s = 4;
    
    	for (BigUnsigned i = 2; i < p; ++i)
    		s = (s*s - 2) % m;
    
    	return (s == 0);
    }
    
    // division test with first prime numbers
    bool divisionPrecheck(vector<bool>& primes, BigUnsigned p)
    {
    	BigUnsigned m = calculateMersenneNumber(p);
    
    	for (BigUnsigned divisor = 3; divisor <= DIVISOR_PRECHECK_MAX; ++divisor)
    	{		
    		if (primes[convertBigToU64(divisor)])
    		{
    			if (m % divisor == 0)
    			{
    				return false;
    			}
    		}
    	}
    
    	return true;
    }
    
    int main()
    {
    	uint64_t startCount = 0, lastCount = 0, delta = 0, frequency;
    	QueryPerformanceFrequency((LARGE_INTEGER*)&frequency);
    	cout << "cpu frequency: " << frequency << " Hz" << endl << endl;
    
    	////////////////////////// Generate primes lookup table /////////////////////////////////////////////////////////////
    
    	uint64_t const primeLimit = DIVISOR_PRECHECK_MAX;
    	cout << "We generate vector<bool>(primes) up to: " << primeLimit << endl;
    	vector<bool> primes(primeLimit + 1); // calculated primes (low values)
    	vector<bool>& refPrimes = primes;
    
    	primes[0] = primes[1] = false;
    	for (uint64_t i = 2; i < primeLimit + 1; ++i) { primes[i] = true; }
    
    	// Erastothenes sieving loop
    	uint64_t iMax = (uint64_t)sqrt((double)primeLimit) + 1;
    	for (uint64_t i = 2; i < iMax; ++i)
    	{
    		if (primes[i])
    		{
    			for (uint64_t j = 2 * i; j < primeLimit + 1; j += i)
    			{
    				primes[j] = false;
    			}
    		}
    	}
    
    	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    
    	BigUnsigned p;
    
    	cout << "Lucas-Lehmer-Test (Mersenne-Zahlen)" << endl; // https://de.wikipedia.org/wiki/Lucas-Lehmer-Test
    
    	QueryPerformanceCounter((LARGE_INTEGER*)&startCount);
    
    	for (BigUnsigned p = 0; p < 10000; ++p)
    	{
    		if (p > 1000)
    		{
    			if (divisionPrecheck(primes, p) == false)
    			{
    				continue;
    			}			
    		}
    
    		if (LucasLehmerTest(p))
    		{			
    			QueryPerformanceCounter((LARGE_INTEGER*)&lastCount);
    			uint64_t delta = lastCount - startCount;
    
    			textcolor(15);
    			cout << "p: " << p;
    			textcolor(2);
    			cout << "\ttime: " << 1000 * delta / frequency << " ms" << endl;
    		}
    		else
    		{
    			//cout << " ist keine Mersenne-Primzahl" << endl;
    		}
    	}	
    }
    
    cpu frequency: 3127021 Hz
    
    We generate vector<bool>(primes) up to: 293
    Lucas-Lehmer-Test (Mersenne-Zahlen)
    p: 3    time: 0 ms
    p: 5    time: 0 ms
    p: 7    time: 0 ms
    p: 13   time: 1 ms
    p: 17   time: 1 ms
    p: 19   time: 1 ms
    p: 31   time: 2 ms
    p: 61   time: 6 ms
    p: 89   time: 14 ms
    p: 107  time: 23 ms
    p: 127  time: 39 ms
    p: 521  time: 5895 ms
    p: 607  time: 10473 ms
    p: 1279 time: 90122 ms
    p: 2203 time: 326176 ms
    p: 2281 time: 369233 ms
    p: 3217 time: 1261672 ms
    p: 4253 time: 3680716 ms
    p: 4423 time: 4260730 ms
    

    Die Zeiten in ms sind die gesamte vergangene Zeit ab Start des Tests. Der optimale max. Primzahl-Teiler für das jeweilige p müsste hier erst noch untersucht werden.


Log in to reply