뤼카-레머 소수판정법

개요

뤼카-레머 소수판정법(Lucas-Lehmer primality test)은 [math]\displaystyle{ p }[/math]가 3 이상인 소수일 때, 메르센 수 [math]\displaystyle{ M_p = 2^p -1 }[/math]가 소수인지 판정할 수 있게 하는 정리이다.

이 정리는 에두아르 뤼카가 1878년에 구성하고 데릭 레머(Derrick Henry Lehmer)가 1930년에 개량하였다.

진술

수열 [math]\displaystyle{ (r_n)_{n\ge 1} }[/math]을 다음과 같이 정의하자.

[math]\displaystyle{ r_n=\begin{cases} 4,&\text{if }n=1\\ r_{n-1}^2-2,&\text{otherwise} \end{cases} }[/math]

그러면 [math]\displaystyle{ p }[/math]가 3 이상인 소수일 때, [math]\displaystyle{ M_p=2^p -1 }[/math]이 소수일 필요충분조건은 [math]\displaystyle{ r_{p-1}\equiv 0\pmod{M_p} }[/math]인 것이다.

이 수열의 처음 부분을 나열하면 {4, 14, 194, 37634, 1416317954, …}이다.

구현 예시

C

다음은 C언어를 사용하여 뤼카-레머 소수판정법을 구현한 것이다. p의 값이 64까지 작동하는 것이 확인되었다. 아래 코드에서는 p의 값이 3 이상 64 이하인지를 체크한 뒤, 소수라면 1을, 소수가 아니라면 0을 리턴한다.

#include <stdio.h>
#include <stdint.h>
#include <assert.h>

int myprimary(int p)
{
	uint64_t M = (1ULL << p) - 1;
	uint64_t r = 4;
	int i;

	assert(p >= 3);
	assert(p <= 64);
	if (p == 64) M = UINT64_MAX;

	for (i = 1; i < p - 1; i++)
	{
		r = (r * r - 2) % M;
		printf("%d: %llu\n", i, r);
	}

	return (!r);
}

R

R 코딩을 이용해 뤼카-레머 소수판정법을 적용해보자. 새 스크립트를 만들어 다음과 같이 입력한 뒤 가동하자.

myprimality=function(x){
# 초깃값 정의로 r(1)=4에 해당함.
  r=4
  i=1
# 반복구문으로 r(2), r(3), ..., r(p-1)을 M(p)로 나눈 나머지를 구하고 순차적으로 출력.
  while(i<x-1){
    r=(r^2-2)%% (2^x-1)
    i=i+1
    print(r)
  }
# r(p-1)의 값이 0인지 판단하여 메시지 출력.
  if(r==0){
    message("This is prime!")
  }
  else{
    message("This is not prime!")
  }
}

이제 콘솔에 myprimality(17)을 입력하면,

> myprimality(17)
[1] 14
[1] 194
[1] 37634
[1] 95799
[1] 119121
[1] 66179
[1] 53645
[1] 122218
[1] 126220
[1] 70490
[1] 69559
[1] 99585
[1] 78221
[1] 130559
[1] 0
This is prime!

이라는 결과를 얻는다. 즉, [math]\displaystyle{ M_{17}=2^{17}-1=131071 }[/math]은 소수이다. 한편 콘솔에 myprimality(23)을 입력하면,

> myprimality(23)
[1] 14
[1] 194
[1] 37634
[1] 7031978
[1] 7033660
[1] 1176429
[1] 7643358
[1] 3179743
[1] 2694768
[1] 763525
[1] 4182158
[1] 7004001
[1] 1531454
[1] 5888805
[1] 1140622
[1] 4321431
[1] 7041324
[1] 2756392
[1] 1280050
[1] 6563009
[1] 6107895
This is not prime!

를 얻는다. 즉, [math]\displaystyle{ M_{23}=2^{23}-1=8388607 }[/math]은 소수가 아니다.

증명

주어진 점화 관계식 [math]\displaystyle{ r_1=4,\ r_n = r_{n-1}^2-2 }[/math]에서 왜 뜬금없이 -2가 들어가 있는지 의아하게 느껴지겠지만 이것은 사실 일반화된 뤼카 수열과 관계가 있다.

먼저 위 점화식을 푼다: [math]\displaystyle{ \displaystyle r_{n-1}=x+x^{-1} }[/math]로 치환하면 [math]\displaystyle{ \displaystyle r_{n} = (x+x^{-1})^2 -2 = x^{2}+x^{-2} }[/math]이다. 따라서 [math]\displaystyle{ \displaystyle r_1=a+a^{-1} \ (a \ge 1) }[/math]이면 [math]\displaystyle{ \displaystyle r_{n}=a^{2^{n-1}}+a^{-2^{n-1}} }[/math]임을 수학적 귀납법으로 이끌어낼 수 있다. 한편 첫 항의 조건에 따라 [math]\displaystyle{ \displaystyle 4=a+a^{-1} }[/math]이며, 이 방정식을 풀면 [math]\displaystyle{ a=2+\sqrt{3},\ a^{-1}=2-\sqrt{3} }[/math]이다. [math]\displaystyle{ a^{-1}=b }[/math]라 하면 [math]\displaystyle{ \displaystyle r_n= a^{2^{n-1}}+b^{2^{n-1}}= (2+\sqrt{3})^{2^{n-1}}+(2-\sqrt{3})^{2^{n-1}} }[/math]과 같이 표현된다.

이 수열의 [math]\displaystyle{ n }[/math]번째 항은 [math]\displaystyle{ \displaystyle \begin{cases} t_0=2,\ t_1=4 \\ t_n=4t_{n-1}-t_{n-2} \ (n \ge 2) \end{cases} }[/math]로 정의된 수열의 [math]\displaystyle{ 2^{n-1} }[/math]번째 항과 같다. 즉 [math]\displaystyle{ r_n=t_{2^{n-1}} }[/math]. 또한 [math]\displaystyle{ \{t_n\} }[/math]은 일반화된 뤼카 수열인 [math]\displaystyle{ \{V_n(4,1)\} }[/math]과 같다.

필요조건

여기서 사칙연산 및 거듭제곱 등은 환 [math]\displaystyle{ X= \{a+b\sqrt{3}|a,b \in \mathbb{Z}/M_p\mathbb{Z} \} }[/math] 내에서 이루어진다. 즉 [math]\displaystyle{ a+b\sqrt{3} \equiv c+d\sqrt{3} \pmod{M_p} }[/math]라고 써진 건 [math]\displaystyle{ a \equiv c \pmod{M_p},\ b \equiv d \pmod{M_p} }[/math]를 뜻한다고 간주한다.

[math]\displaystyle{ M_p }[/math]가 소수라고 가정하고 [math]\displaystyle{ r_{p-1} \equiv 0 \pmod{M_p} }[/math]임을 보이자. [math]\displaystyle{ p \ge 3 }[/math]일 때, [math]\displaystyle{ M_p \equiv 7 \pmod{8} }[/math] 및 [math]\displaystyle{ M_p \equiv 7 \pmod{12} }[/math]가 성립한다. (이유는 메르센 소수 문서의 성질 문단 참조) 따라서 2는 법 [math]\displaystyle{ M_p }[/math]에 대한 이차잉여이고 3은 이차잉여가 아니다.

이를 르장드르 기호 및 오일러의 규준으로 표현하면 아래와 같다.

[math]\displaystyle{ \displaystyle \begin{cases} 2^{\frac{M_p-1}{2}} \equiv \left(\frac{2}{M_p} \right) \equiv 1 \pmod{M_p} \\ 3^{\frac{M_p-1}{2}} \equiv \left(\frac{3}{M_p} \right) \equiv -1 \pmod{M_p} \end{cases} }[/math]

한편 위에서 언급한 [math]\displaystyle{ a=2+\sqrt{3},\ b=2-\sqrt{3} }[/math]를 가져온다. 식을 변형하면 [math]\displaystyle{ 2a = 4+2\sqrt{3} = (1+\sqrt{3})^2 }[/math]이다. 양 변에 [math]\displaystyle{ \displaystyle \frac{M_p+1}{2} }[/math]제곱을 하면 [math]\displaystyle{ \displaystyle 2^{\frac{M_p+1}{2}}a^{\frac{M_p+1}{2}} = (1+\sqrt{3})^{M_p+1} }[/math]이다. 각 항은 아래와 같이 식을 간단히 할 수 있다.

[math]\displaystyle{ a^{\frac{M_p+1}{2}} = a^{2^{p-1}} }[/math]

[math]\displaystyle{ 2^{\frac{M_p+1}{2}} = 2\cdot 2^{\frac{M_p-1}{2}} \equiv 2 \pmod{M_p} }[/math]

[math]\displaystyle{ (1+\sqrt{3})^{M_p+1} = (1+\sqrt{3})^{M_p}(1+\sqrt{3}) }[/math]

[math]\displaystyle{ \displaystyle (1+\sqrt{3})^{M_p} = \sum_{k=0}^{M_p} \binom{M_p}{k} (\sqrt{3})^k \equiv 1+(\sqrt{3})^{M_p} \equiv 1+3^{\frac{M_p-1}{2}}\sqrt{3} \equiv 1-\sqrt{3} \pmod{M_p} }[/math]

[math]\displaystyle{ (1+\sqrt{3})^{M_p+1} \equiv (1-\sqrt{3})(1+\sqrt{3}) \equiv -2 \pmod{M_p} }[/math]

위 식들을 대입하면 [math]\displaystyle{ \displaystyle 2a^{2^{p-1}} \equiv -2 \pmod{M_p} }[/math]이 된다. 양 변을 2로 나누고 [math]\displaystyle{ b^{2^{p-2}} }[/math]를 곱하면 [math]\displaystyle{ \displaystyle a^{2^{p-1}}b^{2^{p-2}} \equiv -b^{2^{p-2}} \pmod{M_p} }[/math]이다. [math]\displaystyle{ ab=1 }[/math]이므로 [math]\displaystyle{ \displaystyle a^{2^{p-2}} \equiv -b^{2^{p-2}} \pmod{M_p} }[/math]이다.

따라서 우변을 이항하면 [math]\displaystyle{ \displaystyle r_{p-1} = a^{2^{p-2}} + b^{2^{p-2}} \equiv 0 \pmod{M_p} }[/math]임을 알 수 있다.

충분조건

[math]\displaystyle{ r_{p-1} \equiv 0 \pmod{M_p} }[/math]가 성립한다고 가정하고 [math]\displaystyle{ M_p }[/math]가 소수임을 보이자.

만약 [math]\displaystyle{ M_p }[/math]가 합성수이고 [math]\displaystyle{ q \mid M_p,\ 3 \lt q \lt M_p }[/math]인 소인수 [math]\displaystyle{ q }[/math]가 존재한다면 모순이 발생함을 보이면 된다.

가정의 조건식은 [math]\displaystyle{ r_{p-1} \equiv 0 \pmod{q} }[/math]와 같이 쓸 수 있다. 그리고 이 식을 변형하면 [math]\displaystyle{ \displaystyle a^{2^{p-2}} + b^{2^{p-2}} \equiv 0 \pmod{q} }[/math], 즉 [math]\displaystyle{ a^{2^{p-1}} \equiv -1 \pmod{q} }[/math]이다. 양 변을 제곱하면 [math]\displaystyle{ a^{2^p} \equiv 1 \pmod{q} }[/math]이다.

[math]\displaystyle{ \lambda }[/math]를 [math]\displaystyle{ a^{\lambda} \equiv 1 \pmod{q} }[/math]를 만족하는 최소 차수라고 하면 [math]\displaystyle{ \lambda \mid 2^p,\ \lambda \nmid 2^{p-1} }[/math]이어야 하므로 [math]\displaystyle{ \lambda = 2^p }[/math]이다.

한편 [math]\displaystyle{ \displaystyle a^q=(2+\sqrt{3})^q \equiv 2^q+3^\frac{q-1}{2}\sqrt{3} \pmod{q} }[/math]에서 페르마의 소정리 에 의해 [math]\displaystyle{ 2^q \equiv 2 \pmod{q} }[/math]이고, [math]\displaystyle{ q\gt 3 }[/math]이므로 [math]\displaystyle{ \displaystyle 3^{\frac{q-1}{2}} \equiv \left(\frac{3}{q}\right) \equiv \pm 1 \pmod{q} }[/math]이다. 즉 [math]\displaystyle{ \displaystyle a^q \equiv 2 \pm \sqrt{3} \equiv a^{\pm 1} \pmod{q} }[/math]이다. 다시 말해 [math]\displaystyle{ a^{q+1} \equiv 1 \pmod{q} }[/math] 또는 [math]\displaystyle{ a^{q-1} \equiv 1 \pmod{q} }[/math]이므로 [math]\displaystyle{ \lambda \mid q+1 }[/math] 또는 [math]\displaystyle{ \lambda \mid q-1 }[/math]이다. 그러므로 [math]\displaystyle{ \lambda = 2^p \le q+1 }[/math]이고, [math]\displaystyle{ q \ge 2^p-1 = M_p }[/math]이다.

하지만 이는 앞서 가정했던 [math]\displaystyle{ q\lt M_p }[/math]와 모순이다. 따라서 [math]\displaystyle{ M_p }[/math]는 소수이다.

야코비 오류 검사

큰 메르센 소수에 대해 뤼카-레머 테스트를 적용할 때, 컴퓨터로 반복 연산을 구현한다. 그런데 지수가 클수록 변수 하나의 메모리도 커지고, 연산 횟수도 늘어나는데 이 과정에서 하드웨어 오류가 하나도 없어야 한다. 가령 [math]\displaystyle{ p=57885161 }[/math]일 때, 위에서 구현한 코드에서 변수의 길이는 57885161비트이고, 반복문의 횟수는 [math]\displaystyle{ p-2 }[/math]회, 즉 57885159회이다. 여기서 컴퓨터가 (오버클럭이나 전원 단절과 같은 문제로) 한 번이라도 계산 실수를 한다면 모든 계산 과정은 수포로 돌아간다.

따라서 실제로는 반복문 중간 지점마다 적절히 세이브 파일을 저장해 두고, 아래 방법으로 검사를 했더니 오류가 발견되면 해당 체크 포인트로 돌아가는 방식을 쓴다.

검사 방법

이 검사 방법은 야코비 기호를 이용하여 셈한다고 해서 야코비 오류 검사(Jacobi error checking)라고 부른다. 실제로 GIMPS 및 Prime95에서 이 검산을 사용하고 있다.

위에서 정의한 수열[math]\displaystyle{ \{r_n\} }[/math]의 특정 항과 메르센 수 [math]\displaystyle{ M_p (p \ge 3) }[/math]를 임의로 가져온다. [math]\displaystyle{ R_k }[/math]는 [math]\displaystyle{ r_k }[/math]를 [math]\displaystyle{ M_p }[/math]로 나눈 나머지로, 이 값이 메모리에 실제로 저장되어 있는 값이다.

임의의 번호 [math]\displaystyle{ k \ge 2 }[/math]에 대해 [math]\displaystyle{ R_k }[/math]가 오류 없이 제대로 계산되었다면, 아래 두 식을 만족해야 한다. 여기서 괄호 표시는 야코비 기호이다.

[math]\displaystyle{ \displaystyle \left(\frac{R_k+2}{M_p}\right) = 1,\ \left(\frac{R_k-2}{M_p}\right) = -1 }[/math]

반대로 말하면 이 등식을 만족하지 못한다면 계산 중간에 오류가 있었다는 뜻이 된다. 물론 이 검산만으로는 오류 여부를 온전히 밝힐 수는 없지만 상당 부분 계산 오류를 잡아낼 수 있다.

증명

먼저 [math]\displaystyle{ k=2, R_k \equiv 14 \pmod{M_p} }[/math]인 경우부터 알아본다.

[math]\displaystyle{ \displaystyle \left(\frac{R_k+2}{M_p}\right) = \left(\frac{16}{M_p}\right) = \left(\frac{4}{M_p}\right)^{2} = 1 }[/math]

[math]\displaystyle{ \displaystyle \left(\frac{R_k-2}{M_p}\right) = \left(\frac{12}{M_p}\right) = \left(\frac{2}{M_p}\right)^{2} \left(\frac{3}{M_p}\right) = \left(\frac{3}{M_p}\right) = A }[/math] 여기서 [math]\displaystyle{ M_p \equiv 3 \pmod{4}, M_p \equiv 1 \pmod{3} }[/math]이므로 [math]\displaystyle{ A = -\left(\frac{M_p}{3}\right) = -\left(\frac{1}{3}\right) = -1 }[/math]

다음으로 [math]\displaystyle{ R_k }[/math]가 명제를 만족할 때, [math]\displaystyle{ R_{k+1} \equiv R_k^2-2 \pmod{M_p} }[/math]도 마찬가지로 참임을 보인다.

주어진 조건: [math]\displaystyle{ \displaystyle \left(\frac{R_k+2}{M_p}\right) = 1,\ \left(\frac{R_k-2}{M_p}\right) = -1 }[/math]

[math]\displaystyle{ \displaystyle \left(\frac{R_{k+1}+2}{M_p}\right) = \left(\frac{R_k^2}{M_p}\right) = \left(\frac{R_k}{M_p}\right)^{2} = 1 }[/math]

[math]\displaystyle{ \displaystyle \left(\frac{R_{k+1}-2}{M_p}\right) = \left(\frac{R_k^2-4}{M_p}\right) = \left(\frac{R_k+2}{M_p}\right) \left(\frac{R_k-2}{M_p}\right) = 1 \cdot -1 = -1 }[/math]

수학적 귀납법에 의해 모든 [math]\displaystyle{ k \ge 2 }[/math]와 임의의 [math]\displaystyle{ M_p (p \ge 3) }[/math]에 대해 주어진 관계식은 성립한다.

활용

13번째 메르센 소수인 [math]\displaystyle{ M_{521}=2^{521}-1 }[/math]과 그 이후의 소수는 모두 뤼카-레머 소수판정법으로 소수 여부가 밝혀졌으며, GIMPS에서 제공하는 프로그램인 Prime95는 메르센 수의 소수 여부를 판정할 때 뤼카-레머 소수판정법으로 교차검증을 한다. (서로 독립된 컴퓨터에서 각각 돌려서 결과 일치를 확인하는 방식)

소요 시간

앞서 소개된 구현 예시에서 반복문의 횟수와 변수의 길이(이진법 자릿수)는 각각 지수 [math]\displaystyle{ p }[/math]에 비례해서 커진다. 한 번 제곱(r*r)하고 나머지 연산(%)을 하는 데는 일반적인 곱셈법으로 지수의 제곱에 비례하며, 현재까지 알려진 가장 빠른 방법으로는 [math]\displaystyle{ p \log p }[/math]에 비례한다. 따라서 해당 알고리즘을 완주하는데 걸리는 시간은 지수의 제곱 이상으로 비례해서 커진다.

Prime95에서 뤼카-레머 테스트를 일반 컴퓨터로 돌릴 때, [math]\displaystyle{ p\lt 10^6 }[/math]일 때는 테스트 시간이 5분도 채 안 걸리지만 [math]\displaystyle{ p \gt 10^7 }[/math] 크기에서는 한나절 이상 잡아먹고, [math]\displaystyle{ p \gt 10^8 }[/math]의 경우 전원을 끄지 않고 4주 이상 연속으로 가동해야 한다.

하지만 이 정리는 필요조건 및 충분조건 둘 다 증명이 되어 있고, 현재로서는 다른 소수판정법보다는 압도적으로 빠른 알고리즘이다. 실행 도중 하드웨어 오류만 없다면 어떤 소수 지수에 대해 메르센 수의 소수 여부를 사람이 감당할 수 있는 시간 내에 정확하게 가려낼 수 있는 것이다. 알려진 가장 큰 소수의 지위를 메르센 소수가 장기간 차지하고 있는 것도 이 때문이다.

사실 2의 1억 제곱 크기의 거대한 수를 수개월 내로 검사할 수 있다는 것은 기적 수준으로 빠른 셈이다. 더욱이 위에서 말한 소요시간은 가정용 컴퓨터 한 대 기준이고, 슈퍼컴퓨터에 버금가는 엔진 여러 대를 마련해 둔다면 시간을 더 단축할 수 있다.