Test Lucasa-Lehmera

Test Lucasa-Lehmera – test pierwszości dla liczb Mersenne’a. Test został ułożony przez Edwarda Lucasa w 1856, a następnie ulepszony przez niego w 1878. W 1930 test został zmodyfikowany przez Derricka Henry’ego Lehmera.

Niech $M_p$ oznacza liczbę Mersenne’a $M_p=2^p-1$ dla pewnej nieparzystej liczby pierwszej $p$ (tzn. liczby pierwszej większej od 2). Pierwszość liczby $p$ można sprawdzić za pomocą prostego algorytmu podziału, gdzie $p$ jest wykładniczo mniejsza od $M_p.$ Definiuje się następujący ciąg liczb naturalnych $(S_k):$

$S_k=\{\begin{array}{cc}4& \text{gdy }k=0,\\ S_{k-1}^2-2& \text{gdy }k\ne 0.\end{array}$

Oto kilka początkowych wyrazów tego ciągu: 4, 14, 194, 37634, ... (ciąg A003010 w OEIS).

Test Lucasa-Lehmera orzeka, że liczba $M_p$ jest pierwsza wtedy i tylko wtedy, gdy jest dzielnikiem wyrazu o numerze $(p-2)$ w tym ciągu, co krótko zapisuje się kongruencją:

$S_{p-2}\equiv 0(\mathrm{mod}{M}_p).$

Resztę z dzielenia liczby $S_{p-2}$ przez $M_p$ nazywa się residuum Lucasa-Lehmera liczby $p.$ Istotę testu można zatem streścić sformułowaniem:

liczba Mersenna $M_p$ jest pierwsza wtedy i tylko wtedy, gdy residuum Lucasa-Lehmera liczby $p$ równe jest zeru.

Za pomocą pseudokodu test można przedstawić w następujący sposób:

// Ustalamy, czy Mp = 2p − 1 jest pierwsza
Lucas-Lehmer(p)
    niech s ← 4
    niech M ← 2p − 1
    powtórz p − 2 razy:
        s ← ((s × s) − 2) mod M
    jeśli s = 0 zwróć PIERWSZA w przeciwnym razie zwróć ZŁOŻONA

Działanie mod M w każdej iteracji zapewnia, że wszystkie pośrednie wyniki są co najwyżej $p$-bitowe (w innym przypadku liczba bitów byłaby dublowana w każdej iteracji).

Złożoność czasowa

W powyższym algorytmie podczas każdej iteracji wykonywane są dwie kosztowne operacje: mnożenie s × s i operacja modulo mod M. Operacja mod M może być wykonywana szczególnie efektywnie na standardowych dwójkowych systemach poprzez dokonanie następującej obserwacji

$k\equiv (kmod{2}^n)+\lfloor k/2^n\rfloor (\mathrm{mod}{2}^n-1).$

Mówi ona, że najmniej znaczące $n$ bity spośród $k$ plus pozostałe bity $k,$ są równoważne $k$ modulo $2^n-1.$ Równoważność może być stosowana powtarzalnie do momentu co najwyżej $n$ bitów reszty. W tym postępowaniu reszta po wykonaniu dzielenia $k$ przez liczbę Mersenne’a $2^n-1$ jest obliczona bez używania dzielenia. Na przykład:

$\begin{array}{cc}916{mod2}^5-1& =1110010100_2{mod2}^5-1\\ & =11100_2+10100_2{mod2}^5-1\\ & =110000_2{mod2}^5-1\\ & =1_2+10000_2{mod2}^5-1\\ & =10001_2{mod2}^5-1\\ & =10001_2\\ & =17.\end{array}$

Ponadto, ponieważ s × s nigdy nie przekroczy $M^2<2^{2p},$ ta prosta technika zbiega w co najwyżej 1 $p$-bitowy dodatku (ewentualnie przeprowadza z $p$-tego bitu w pierwszy bit), co może być wykonane w liniowym czasie. Algorytm ten posiada wyjątkowy przypadek. Daje $2^n-1$ dla mnożenia modulo zamiast poprawnej wartości 0. Jednak ten przypadek jest łatwy do wykrycia i naprawienia.

Przykłady

Liczba Mersenne’a M₃ = 7 jest pierwsza. Test Lucasa-Lehmera sprawdza to w następujący sposób. Początkowo $s$ jest ustalona i równa 4, później jest modyfikowana 3−2 = 1 raz:

• s ← ((4 × 4) − 2) mod 7 = 0.

Ponieważ końcowa wartość $s$ wynosi 0, wnioskujemy, że M₃ jest pierwsza.

Z drugiej strony, M₁₁ = 2047 = 23 × 89 nie jest pierwsza. Powtórnie, $s$ jest ustalona i równa 4, ale teraz jest modyfikowana 11−2 = 9 razy:

• s ← ((4 × 4) − 2) mod 2047 = 14
• s ← ((14 × 14) − 2) mod 2047 = 194
• s ← ((194 × 194) − 2) mod 2047 = 788
• s ← ((788 × 788) − 2) mod 2047 = 701
• s ← ((701 × 701) − 2) mod 2047 = 119
• s ← ((119 × 119) − 2) mod 2047 = 1877
• s ← ((1877 × 1877) − 2) mod 2047 = 240
• s ← ((240 × 240) − 2) mod 2047 = 282
• s ← ((282 × 282) − 2) mod 2047 = 1736

Ponieważ końcowa wartość $s$ nie jest równa  0, wnioskujemy, że M₁₁ = 2047 nie jest pierwsza. Mimo że M₁₁ = 2047 posiada nietrywialne czynniki, test Lucasa-Lehmera nie daje żadnych wskazówek, jakie mogą one być.

Dowód

Dowód poprawności testu przedstawiany tutaj jest prostszy niż oryginalny dowód podany przez Lehmera. Przywołajmy definicję

$s_i=\{\begin{array}{cc}4& \text{gdy }i=0,\\ s_{i-1}^2-2& \text{w pozostałych przypadkach.}\end{array}$

Naszym celem jest pokazanie, że $M_p$ jest pierwsza wtedy i tylko wtedy, gdy $s_{p-2}\equiv 0(\mathrm{mod}{M}_p).$

Ciąg $⟨s_i⟩$ jest dany rekurencyjnie z rozwiązaniem w jawnej postaci. Niech $\omega =2+\sqrt{3}$ oraz $\bar{\omega }=2-\sqrt{3}.$ Poprzez indukcję dostajemy $s_i={\omega }^{2^i}+\stackrel{2^i}{\bar{\omega }}$ dla każdego $i:$

$s_0={\omega }^{2^0}+\stackrel{2^0}{\bar{\omega }}=(2+\sqrt{3})+(2-\sqrt{3})=4$

i

$\begin{array}{cc}{s}_n& =s_{n-1}^2-2\\ & ={({\omega }^{2^{n-1}}+\stackrel{2^{n-1}}{\bar{\omega }})}^2-2\\ & ={\omega }^{2^n}+\stackrel{2^n}{\bar{\omega }}+2(\omega \bar{\omega }{)}^{2^{n-1}}-2\\ & ={\omega }^{2^n}+\stackrel{2^n}{\bar{\omega }}.\end{array}$

Ostatni krok wykorzystuje $\omega \bar{\omega }=(2+\sqrt{3})(2-\sqrt{3})=1.$ Jawna forma jest używana zarówno w dowodzie dostateczności, jak i konieczności.

Dostateczność

Celem jest pokazanie, że $s_{p-2}\equiv 0(\mathrm{mod}{M}_p)$ implikuje, że $M_p$ jest pierwsza. Bezpośredni dowód wykorzystujący elementarną teorię grup daną przez J. W. Bruce, jak nawiązuje Jason Wojciechowski.

Przypuśćmy, że $s_{p-2}\equiv 0(\mathrm{mod}{M}_p).$ Wówczas

${\omega }^{2^{p-2}}+\stackrel{2^{p-2}}{\bar{\omega }}=k{M}_p$

dla pewnego całkowitego $k,$ więc

${\omega }^{2^{p-2}}=k{M}_p-\stackrel{2^{p-2}}{\bar{\omega }}.$

Mnożąc przez ${\omega }^{2^{p-2}},$ otrzymujemy

${\left({\omega }^{2^{p-2}}\right)}^2=k{M}_p{\omega }^{2^{p-2}}-(\omega \bar{\omega }{)}^{2^{p-2}}.$

Zatem

${\omega }^{2^{p-1}}=k{M}_p{\omega }^{2^{p-2}}-1.(1)$

Dla uzyskania sprzeczności, przypuśćmy, że $M_p$ jest złożona, i niech $q$ będzie najmniejszym pierwszym czynnikiem liczby $M_p.$ Liczby Mersenne’a są nieparzyste, więc $q>2.$ Niech ${\mathbb{Z}}_q$ będzie zbiorem liczb całkowitych modulo $q$ i niech $X=\{a+b\sqrt{3}\mid a,b\in {\mathbb{Z}}_q\}.$ Mnożenie w $X$ jest zdefiniowane, jak następuje

$(a+\sqrt{3}b)(c+\sqrt{3}d)=[(ac+3bd)modq]+\sqrt{3}[(ad+bc)modq].$

Oczywiście to mnożenie jest działaniem wewnętrznym, to znaczy produkt liczb z $X$ przechodzi w siebie. Rozmiar $X$ oznaczamy przez $|X|.$

Gdy $q>2,$ $\omega$ i $\bar{\omega }$ należą do $X$. Podzbiór elementów $X$ z ich odwrotnościami tworzy grupę, którą oznaczamy $X^{*},$ a jej rząd $|X^{*}|.$ Jeden z elementów $X$ nie posiada odwrotności, jest to 0, więc $|X^{*}|⩽|X|-1=q^2-1.$

Teraz $M_p\equiv 0(\mathrm{mod}q)$ i $\omega \in X,$ więc

$k{M}_p{\omega }^{2^{p-2}}=0$

w $X.$

Następnie, wykorzystując równanie (1),

${\omega }^{2^{p-1}}=-1$

w $X$ i podnosząc obie strony do kwadratu, otrzymujemy

${\omega }^{2^p}=1.$

Zatem $\omega$ należy do $X^{*}$ i posiada odwrotność ${\omega }^{2^p-1}.$ Co więcej, rząd $\omega$ dzieli $2^p$, jednak ${\omega }^{2^{p-1}}\ne 1,$ więc nie dzieli on $2^{p-1}$. Wobec tego rząd $\omega$ wynosi dokładnie $2^p.$

Rząd elementu jest nie większy niż rząd (rozmiar) grupy, więc

$2^p⩽|X^{*}|⩽q^2-1<q^2.$

Ale $q$ jest najmniejszym pierwszym czynnikiem rozkładu liczby złożonej $M_p,$ więc

$q^2⩽M_p=2^p-1.$

To prowadzi do sprzeczności $2^p<2^p-1.$ Dlatego $M_p$ jest pierwsza.

Konieczność

Z drugiej strony, celem jest pokazanie, że pierwszość $M_p$ implikuje $s_{p-2}\equiv 0(\mathrm{mod}{M}_p).$ Poniższy uproszczony dowód przedstawił Öystein J.R.

Gdy $2^p-1\equiv 7(\mathrm{mod}12)$ dla nieparzystych $p>1,$ z własności symbolu Legendre’a wynika, że $(3|M_p)=-1.$ To oznacza, że 3 jest podwójnym nieresiduum modulo $M_p.$ Dzięki kryterium Eulera, jest to równoważne

$3^{\frac{M_p-1}{2}}\equiv -1(\mathrm{mod}{M}_p).$

Natomiast 2 jest podwójnym residuum modulo $M_p,$ ponadto $2^p\equiv 1(\mathrm{mod}{M}_p)$ i stąd $2\equiv 2^{p+1}={\left(2^{\frac{p+1}{2}}\right)}^2(\mathrm{mod}{M}_p).$ Tym razem kryterium Eulera pozwala napisać

$2^{\frac{M_p-1}{2}}\equiv 1(\mathrm{mod}{M}_p).$

Połączenie tych dwóch równoważnych relacji daje

$24^{\frac{M_p-1}{2}}\equiv {\left(2^{\frac{M_p-1}{2}}\right)}^3\left(3^{\frac{M_p-1}{2}}\right)\equiv (1{)}^3(-1)\equiv -1(\mathrm{mod}{M}_p).$

Niech $\sigma =2\sqrt{3}$ i zdefiniujmy $X$ tak jak wcześniej jako pierścień $X=\{a+b\sqrt{3}\mid a,b\in {\mathbb{Z}}_{M_p}\}.$ Wówczas w pierścieniu $X$ otrzymujemy

$\begin{array}{cc}(6+\sigma {)}^{M_p}& =6^{M_p}+\left(2^{M_p}\right)\left({\sqrt{3}}^{M_p}\right)\\ & =6+2\left(3^{\frac{M_p-1}{2}}\right)\sqrt{3}\\ & =6+2(-1)\sqrt{3}\\ & =6-\sigma ,\end{array}$

gdzie pierwsza równość wykorzystuje dwumian Newtona w skończonej dziedzinie, która przedstawia się

$(x+y{)}^{M_p}\equiv x^{M_p}+y^{M_p}(\mathrm{mod}{M}_p),$

a druga równość wykorzystuje małe twierdzenie Fermata, które brzmi

$a^{M_p}\equiv a(\mathrm{mod}{M}_p)$

dla każdej liczby całkowitej $a.$ Wartość $\sigma$ została wybrana tak, aby $\omega =\frac{(6+\sigma {)}^2}{24}.$ Co za tym idzie, może być wykorzystana do wyliczenia ${\omega }^{\frac{M_p+1}{2}}$ w pierścieniu $X$ jako

$\begin{array}{cc}{\omega }^{\frac{M_p+1}{2}}& =\frac{(6+\sigma {)}^{M_p+1}}{24^{\frac{M_p+1}{2}}}\\ & =\frac{(6+\sigma )(6+\sigma {)}^{M_p}}{24\cdot 24^{\frac{M_p-1}{2}}}\\ & =\frac{(6+\sigma )(6-\sigma )}{-24}\\ & =-1.\end{array}$

To co pozostaje, to mnożenie obu stron równania przez $\stackrel{\frac{M_p+1}{4}}{\bar{\omega }}$ i wykorzystanie $\omega \bar{\omega }=1,$ co daje

$\begin{array}{cc}{\omega }^{\frac{M_p+1}{2}}\stackrel{\frac{M_p+1}{4}}{\bar{\omega }}& =-\stackrel{\frac{M_p+1}{4}}{\bar{\omega }}\\ {\omega }^{\frac{M_p+1}{4}}+\stackrel{\frac{M_p+1}{4}}{\bar{\omega }}& =0\\ {\omega }^{\frac{2^p-1+1}{4}}+\stackrel{\frac{2^p-1+1}{4}}{\bar{\omega }}& =0\\ {\omega }^{2^{p-2}}+\stackrel{2^{p-2}}{\bar{\omega }}& =0\\ s_{p-2}& =0.\end{array}$

Skoro $s_{p-2}$ jest 0 w $X,$ jest również 0 modulo $M_p.$

Zastosowanie

Test Lucasa-Lehmera jest testem pierwszości używanym przez Great Internet Mersenne Prime Search do znajdowania dużych liczb pierwszych. Te poszukiwania są bardzo owocne i przyczyniły się do znalezienia wielu największych liczb pierwszych znanych dzisiaj. Test jest uważany za wartościowy, ponieważ potrafi zbadać pierwszość olbrzymich liczb zawartych w zbiorach o dużej liczbie elementów w przeciągu przystępnego wymiaru czasu. W odróżnieniu od równie szybkiego testu Pépina dla każdej liczby Fermata, który może być używany jedynie na wiele mniejszych zbiorach tak olbrzymich liczb, zanim zakres obliczeniowy się wyczerpie.

Uwagi: test jest bardzo szybki i bardzo prosty. Przy jego użyciu znaleziono największe dotąd znane liczby pierwsze. Test wymaga jak najszybszego algorytmu mnożenia np. szybkiej transformaty Fouriera. Ponieważ mnożenie liczb zapisanych w systemie liczbowym o danej podstawie jest w pewnym sensie splotem funkcji (cyfra jako wartość w funkcji potęgi podstawy), a transformata Fouriera splotu dwóch funkcji jest iloczynem ich transformat, zatem wielkie liczby naturalne można mnożyć szybko przez siebie, robiąc szybką transformatę Fouriera z ich cyfr, a następnie szybką transformatę odwrotną iloczynu tych transformat.

Implementacje

Kod w Mathematica

Krótki kod w języku Mathematica (tutaj sprawdzający liczbę pierwszą Mersenne’a Davida Slowinskiego & Paula Gage’a z 4 stycznia, 1994 (znaną 33.) z wykładnikiem $p=859433$) w czasie nieco ponad 2 godzin na komputerze osobistym z procesorem Skylake o częstotliwości 2,9 GHz): Funkcja Mod[s, M] zwraca na końcu 0 jedynie jeśli $M$ jest liczbą pierwszą.

M = 2^859433 - 1;
s = 4;
AbsoluteTiming[
 Monitor[Do[
   s = Mod[s*s - 2, M], {n, 1, 859433 - 2}], {ProgressIndicator[
    n, {1, 859433 - 2}], n}]];
Print[Mod[s, M]]

Kod w Python

Przykładowy kod w języku Python 3 wykorzystujący bibliotekę gmpy2, kod zawiera optymalizację dzielenia modulo (patrz wyżej – Złożoność czasowa) oraz monitorowanie pozostałego czasu wykonania (estymowanego).

'''
Dla liczby p > 2, funkcja zwraca True jeżeli 2^p-1 jest pierwsza, w innym razie zwraca False
'''
from gmpy2 import *
from datetime import *

def Lucas_Lehmer(p):
    s = xmpz(4)
    M = xmpz(2 ** p) - 1
    startt = datetime.now()
    for i in range(0, p - 2):
        time_elapsed = datetime.now() - startt
        speed = i / (time_elapsed.total_seconds() + 1)
        remaining = (p - 2 - i) / (speed + 1)
        if i % int(speed + 1) == 0:
            print('Czas pozostały (hh:mm:ss.ms) {} czas wykonywania {}'.format(timedelta(seconds=remaining), time_elapsed))
        if (s.bit_length() <= p):
            # Podstawowa wersja algorytmu dla s o bitach <= p
            s = s * s - 2
            s = f_mod(s, M)
        else:
            # Optymalizacja dla s bitów > p
            s = s * s - 2
            md = s.bit_length() % 2
            right = s.bit_length() // 2 + md
            left = s.bit_length() // 2
            s = s[0:right] + f_mod(s[left:], M)

    if s == 0:
        return True
    else:
        return False

Zobacz też

• Édouard Lucas
• Derrick Henry Lehmer
• GIMPS

Linki zewnętrzne

• Szablon:MathWorld
• Szablon:Cytuj stronę
• Szablon:Cytuj stronę
• Szablon:Cytuj stronę