515. 프로젝트 오일러 #515 : 불협화음 숫자들

문제 : 

Let d(p,n,0) be the multiplicative inverse of n modulo prime p, defined as n × d(p,n,0) = 1 mod p.

Let  for k ≥ 1.

Let for all primes a ≤ p < a + b.

You are given:

• D(101,1,10) = 45
• D(10³,10²,10²) = 8334
• D(10⁶,10³,10³) = 38162302

Find D(10⁹,10⁵,10⁵).

문제를 보면, 정수론에 입각한 문제입니다.  현재 신상 문제로, 아직 푼 사람들 수가 100명이 채 안 되네요.

정수론 문제가 나오면, 많이들 어려워하시기 마련입니다.  (어떤 면에서 저는 정수론 문제가 더 잘 와닿는 것 같네요.)

수학적 베이스를 보자면,

1) 모듈라 필드에서 곱셈에 대한 역원을 찾는 것은 쉽지 않다.

2) 모듈라하는 숫자가 소수라면, 0을 제외한 모든 수는 역원을 가지고 있다.

예를 들면 7-모듈라라고 하면, 우리는 다음과 같은 형태로 역원을 찾을 수 있습니다.

1x1 = 1 (mod 7)

2x4 = 1 (mod 7)

3x5 = 1 (mod 7)

4x2 = 1 (mod 7)

5x3 = 1 (mod 7)

6x6 = 1 (mod 7)

어떤 소수 p이든, 1과 p-1은 그 자신이 역원이 됩니다.  이것은 정수에서 1과 -1이 자기 자신이 곱셈의 역원이 되는 것과 갈다고 생각하셔도 됩니다.

그런데 이 문제에서 정의된 d(p, n, k)가 가장 중요한 부분입니다.  그리고 문제에서는 언급하지 않았지만, n = p-1 이라는 것이죠.

형태로 자기호출 함수의 형태를 띄고 있습니다.  k가 줄어들어서 0이 될 때까지 엄청난 호출을 할텐데요.  당연히 이것을 곧이 곧대로 자기호출 함수를 이용해서 프로그램 짜시면, 아주 작은 수의 범위(100 이내)에서는 동작할지 몰라도, 조금만 숫자가 커지면, 대책이 없습니다.

우리가 관심가져야할 부분은 d(p, n, 0)이 과연 몇번이나 나올 것인가이죠.  위의 수식을 풀어서 표시하면 다음과 같습니다.

이 된다는 것입니다.

이것을 이용해서 k 에 대해서 배열을 선언하면,

k = 1

1 1 1 1 1 1 1 1 ...... 1

k = 2

1 2 3 4 5 6 7 8 ...... n

k = 3

1 3 6 10 15 21 28 36 ....... ??

과 같이 자신 이전의 값과 바로 전에 계산된 값을 합치면, d(p, n, 0)이 얼마나 많이 나타났는지를 알 수 있습니다.

이것을 미리 배열로 저장하면 편하겠지만, 문제는 이 숫자가 엄청 크기 때문에 64bit 정수형으로 표시할 수 없다는 것입니다.  또한 문제에서 10의 9승의 범위 숫자만큼 위 배열을 할당하는 것은 어렵습니다.  결국 어느 한가지도 해결이 안 된다는 것이죠.

사실 위의 것을 잘 유도하면 하나의 식으로 만들 수 있습니다.  (그렇다고 사실 나아지는 것을 전혀 없습니다.  매번 k수만큼 계산해야 하니까요.)

조합을 계산하기 위해서는 k-1번 곱셈과 나눗셈을 해주어야 합니다.  이것도 k가 큰 수라면 오래 걸립니다.

이것을 좀 더 빨리 하기 위해서, 문서에서 검색을 할 때 사용하는 기법을 사용할 수 있습니다.  이렇게 하면, 매 계산마다, 한번의 곱셈과 한번의 나눗셈이면 되죠.  이 방법에 대해서는 제가 쓴 옛날 글을 참조하세요.  (http://sdev.tistory.com/146)

여기까지 되셨다면, 

이제는 모듈라 p 에서의 역원을 찾는 것입니다.  역원을 찾는 방법은 모든 수를 다 대입해보는 방법도 있겠지만, 이렇게 하면, 10의 9승인 소수 p에 대해서 무려 p-1번 루프를 돌아야 하니 비효율적입니다.

이 계산을 위해서는 유클리드 호젯법과 중국인의 나머지 정리를 사용했습니다.

uint64_t findInverse(uint64_t p, uint64_t m)
{
    static int64_t rec[10000];
    uint64_t a = p;
    int count = 0;

    while( m ) { uint64_t t = m; rec[count++] = a/m; m = a%m; a = t; }
    int64_t s = 0, t = 1;
    count--;
    while( --count >= 0 )
    {
        int64_t r = t;
        t = -t*rec[count]+s;
        s = r;
    }
    t = (t>0)?t:p+t;
    return t;
}

사실 인터넷에 뒤져보면, 아마 역원 구하는 알고리즘이 나와있을 것 같은데, 저는 알고리즘을 수학적 방법으로 제가 알고 있는 그대로를 구현했기 때문에 성능은 보장하지 못 합니다.  그리고 오버플로우가 발생할 수도 있고요.

이렇게 하면, 역원을 구하는 것과 d(p, n, 0)의 갯수를 세는 것이 다 되었습니다.

그런데 이 계산을 하면, 10의 9승 오더의 소수와 10의 5승 범위 때문에 엄청나게 많은 시간을 소모합니다.  이 시간을 현재 이 알고리즘으로는 줄일 수 있는 방법이 없습니다.

알고리즘 최적화는 앞서 말한 것 정도가 최선이 아닐까 생각합니다.

좀 더 나아가고자 한다면, 다음을 고려해보셔야 합니다.

1) d(p, n, k)에서 d(p, i, 0)의 갯수가 모듈라 p에 의해서 달라진다.

2) d(p, i, 0)의 갯수를 모두 다 더하면, 모듈라 p의 값은?

3) 실제 계산해야할 i의 범위는 1 ~ p-1 까지인가?  그보다 더 적어질 수는 없을까?