이항계수를 구하는 알고리즘 고급편 - 페르마의 소정리-

1. 페르마의 소정리

 

p가 소수이면 모든 정수 a에 대하여 $a^p$와 a를 p로 나눈 나머지는 서로 같다.

 

만약 a가 p의 배수가 아닌 서로소라면 $a^{(p-1)}$와 1을 p로 나눈 나머지는 서로 같다

 

 

역은 성립하지 않는다.

 

즉 모든 정수 a에 대하여 $a^p \equiv a (mod p)$임에도 불구하고 p가 합성수일 수 있다

 

 

2. 응용 - 이항계수를 빠르게 구하는 방법

 

n과 r이 10만에서 100만정도 되는 매우 큰 수에 대하여 어떻게 하면 빠르게 구할 수 있을까

 

파스칼의 삼각형을 이용한 다이나믹 프로그래밍으로는 이제는 불가능하다

 

그래서 $nCr = n! / ((n-r)! * r!)$을 직접 구할 필요가 있는데 이 경우 n과 r이 매우 크면 직접 나타내기 어려우므로

 

어떤 소수 p에 대해 나머지를 구하는 방식으로 문제가 나온다

 

---------------------------------------------------------------------------------

 

$nCr = n! / ((n-r)! * r!)$에서 n!을 p로 나눈 나머지를 구하는 것은 어떻게 구하면 되겠지만 문제는 ((n-r)! * r!)으로 나눈 값의 나머지를 구해야한다는 점

 

그렇다고 해서 n!을 p로 나눈 나머지에다가 ((n-r)! * r!)을 p로 나눈 나머지를 서로 나누는 것은 모듈러 연산에서 성립하지 않는다

 

---------------------------------------------------------------------

 

2-1) 페르마의 소정리

 

하지만 페르마의 소정리를 활용하면 가능하다

 

p가 소수이고 a가 p의 배수가 아니라면, $a^{(p-1)} \equiv 1 (mod p)$

 

1)여기서 ((n-r)! * r!)이 p의 배수가 아니라는 점을 확인한다.

 

그거는 생각보다 간단한데 r!이 r, r-1, r-2,.... 1의 곱이고 (n-r)!이 n-r, n-r-1, n-r-2,....,1의 곱으로 만약 여기에 p가 없으면 배수가 아니다. p는 소수니까

 

그래서 $((n-r)!r!)^{(p-1)} \equiv 1 (mod p) $가 성립한다는 점을 확인한다면...

 

$((n-r)!r!)((n-r)!r!)^{(p-2)} \equiv 1 (mod p) $

 

2)여기서 이제 문제가 된다 ((n-r)! * r!)으로 양변을 나눠줘도 되나???

 

-------------------------------------------------------------------------------------------------------------

 

2-2) 모듈로 곱셈의 역원

 

어떤 수 a의 곱셈의 역원이란, a와 곱해서 1이 되는 수를 말한다.

 

mod p에 대하여 곱셈의 역원은 어떻게 정의할 수 있을까?

 

$a * x \equiv 1 (mod p)$를 만족하는 x이다.  그리고 이러한 x가 존재하려면 a와 p가 서로소여야한다. 

 

a와 p가 서로소일 필요충분조건이 a의 mod p에 대한 역원이 존재할 조건이다.

 

그렇다면 a와 p가 서로소라면 합동식의 양변을 a로 나눌 수가 있다

 

이에 대한 따름정리로 p가 만약 소수라면 mod p에 대한 합동식에서 p의 배수가 아닌 임의의 수로 나눌 수 있다

 

---------------------------------------------------------------------------------------------------------------

 

((n-r)! * r!)이 p의 배수가 아니라면, 

 

$((n-r)!r!)((n-r)!r!)^{(p-2)} \equiv 1 (mod p) $은 다음과 같다

 

$((n-r)!r!)^{(p-2)} \equiv ((n-r)!r!)^{-1} (mod p) $

 

양변에 n!을 곱해도 변하지 않으므로

 

$n!((n-r)!r!)^{(p-2)} \equiv n!((n-r)!r!)^{-1} (mod p) $

 

우변은 nCr과 같으므로 $n!((n-r)!r!)^{(p-2)} \equiv nCr (mod p) $

 

그래서 nCr을 p로 나눈 나머지는 $n!((n-r)!r!)^{p-2}$을 p로 나눈 나머지와 같다

 

------------------------------------------------------------------------------------------------------------------

 

2-3) 분할정복을 이용한 거듭제곱

 

다음 문제는 p가 매우 큰 소수라면 $((n-r)!r!)^{p-2}$을 구하는데도 시간이 너무 오래걸린다는 점이다

 

이는 분할정복을 이용한 거듭제곱 알고리즘으로 해결할 수 있다

 

이런 경우 O(nlogp)만에 nCr을 p로 나눈 나머지를 구할 수 있다

 

마지막 문제는 이항계수를 많이 구해야하는 경우 팩토리얼을 어떻게 빨리 계산할 수 있을까인데,

 

미리 1!부터 n!까지 계산해두어 배열에 저장해둔다면, nCr을 계산하는데 시간을 줄일 수 있다

 

 

3. 예시 문제

 

https://www.acmicpc.net/problem/11401

 

3-1) N,K가 400만 정도 되는 매우 큰 수에 대한 이항 계수를 매우 큰 소수로 나눈 나머지를 구하는 문제

 

 

11401번: 이항 계수 3

 

모든 알고리즘을 적용하더라도 마지막으로 신경써야할 점은 거듭제곱이나 팩토리얼을 구하는 과정에서

 

모듈로 연산의 거듭제곱과 곱셈 성질을 이용하여 P로 나눈 나머지를 계속 그 수로 저장시켜줘야한다는 점

 

from sys import stdin

def power(a,n):
    
    result = 1

    while n:
        
        if n % 2 == 1:
            
            result *= a
        
        a *= a
        
        a = a % p
        
        n = n//2
    
    return result

def factorial(n):
    
    result = 1
    
    for i in range(1,n+1):
        
        result *= i
        
        result = result % p
    
    return result

n,r = map(int,stdin.readline().split())

p = 1000000007

print((factorial(n)*power(factorial(n-r)*factorial(r),p-2)) % p)

 

거듭제곱 연산 중간에 a = a%p와 result = result % p를 이용해서 계속 숫자를 줄여나가야한다는 점이 중요하다

 

 

3-2) m개의 수가 한번에 들어올 때 제한시간안에 이항계수를 큰 소수의 나머지를 구하는 문제

 

https://www.acmicpc.net/problem/13977

13977번: 이항 계수와 쿼리

 

m개의 수가 한번에 들어올때는 최댓값에 대한 factorial을 먼저 구하는데.. 이 과정에서도 factorial 점화식이

 

factorial(n-1) * n = factorial(n) 이라는 점을 이용하면 0번, 1번만 1로 저장해놓고 400만까지 for문을 돌면

 

O(400만)만에 모든 factorial을 저장해놓을 수 있다

 

당연히 이 과정에서도 p로 나눈 나머지를 저장해놔야 수를 줄여서 메모리도 아낄 수 있고 시간도 아낄 수 있다 

 

 

from sys import stdin

def power(a,n):
    
    result = 1
    
    while n:
        
        if n % 2 == 1:
            
            result *= a
        
        a *= a
        
        a = a % p
        
        n = n // 2
    
    return result

fact_list = [0] * 4000001

fact_list[0] = 1
fact_list[1] = 1

p = 1000000007

for i in range(2,4000001):
    
    fact_list[i] = fact_list[i-1] * i
    
    fact_list[i] = fact_list[i] % p

    
m = int(stdin.readline())

for _ in range(m):
    
    n,k = map(int,stdin.readline().split())
    
    print((fact_list[n]*power(fact_list[n-k]*fact_list[k],p-2)) % p)

 

참조

 

https://dimenchoi.tistory.com/50

정수론 (4) - 합동식에서의 나눗셈