빠른 거듭제곱

거듭제곱

a의 b승을 c로 나눈 나머지를 구해야 할 때가 있습니다.
참고로 거듭제곱을 영어로는 power라고 합니다.

일단 단순히 구현하면 이렇게 할 수 있습니다.

// integer power function
int ipow(int a, int b, int c){
    int ans = 1;
    for(int i = 0; i < b; ++i) {
        ans *= a;
        ans %= c;
    }
    return ans;
}

한 가지 조심해야 하는 점은, 값의 범위에 따라서는 중간에 ans * a의 값이 int 범위를 넘어갈 수 있습니다.
이 문제를 피하기 위해선 중간에 잠깐 long long으로 캐스팅하면 됩니다.

...
    ans = ((long long) ans * a) % c;
    /* alternatively:
    ans = (ans * 1LL * a) % c;
    */
...

아니면 캐스팅을 하지 말고, 임시 저장 변수를 long long으로 선언할 수도 있습니다.

int ipow(int a, int b, int c){
    long long ans = 1;
    for(int i = 0; i < b; ++i) {
        ans *= a;
        ans %= c;
    }
    return ans;
}

빠른 거듭제곱

그런데 지수 b가 무지막지하게 클 때가 있습니다.
예를 들어 b가 10억 정도라면 위의 코드에서 반복문이 정확히 b번 돌기 때문에 시간이 너무 오래 걸립니다.

이때 b의 홀짝성을 가지고 계산량을 줄이는 방법이 있겠습니다.

$$ a^{2s} = \left(a^s\right)^2 \ a^{2s+1} = \left(a^s\right)^2 \cdot a$$

계산량이 어떻게 줄었는지 감이 오시나요?
a를 2s번 곱할 땐, s번만 곱하고 나면 나머지 절반은 똑같기 때문에, $a^s$만 계산한 후 그 값을 제곱하면 됩니다.
b가 홀수(=2s+1)일 때도 비슷하게 가능하지요.

즉 지수가 $\left \lfloor b/2 \right \rfloor$일 때의 결과를 구한 후, 제곱하고, b의 홀짝성에 따라 a를 곱해주거나 말거나 하면 되겠습니다.
지수가 $\left \lfloor b/2 \right \rfloor$일 때의 그 결과는 어떻게 구할까요?
이것도 그 횟수만큼 돌 필요 없이 절반으로 나누어 처리하면 됩니다. 즉, 재귀적으로 처리하면 됩니다.

결국 재귀 함수를 사용해 이런 코드를 짤 수 있습니다.

int ipow(int a, int b, int c){
    if(b == 0) return 1;

    long long tmp = ipow(a, b/2, c);
    long long ans = (tmp * tmp) % c;
    if(b % 2 == 1) {
        ans *= a;
        ans %= c;
    }

    return ans;
}

b가 10억이었을 때의 호출과정을 생각해봅시다.

ipow(a, 1000000000, c)        -> ipow(a, 500000000, c)      
-> ipow(a, 250000000, c)      -> ipow(a, 125000000, c)      
-> ipow(a, 62500000, c)       -> ipow(a, 31250000, c)       
-> ipow(a, 15625000, c)       -> ipow(a, 7812500, c)        
-> ipow(a, 3906250, c)        -> ipow(a, 1953125, c)        
-> ipow(a, 976562, c)         -> ipow(a, 488281, c)         
-> ipow(a, 244140, c)         -> ipow(a, 122070, c)         
-> ipow(a, 61035, c)          -> ipow(a, 30517, c)          
-> ipow(a, 15258, c)          -> ipow(a, 7629, c)           
-> ipow(a, 3814, c)           -> ipow(a, 1907, c)           
-> ipow(a, 953, c)            -> ipow(a, 476, c)            
-> ipow(a, 238, c)            -> ipow(a, 119, c)            
-> ipow(a, 59, c)             -> ipow(a, 29, c)             
-> ipow(a, 14, c)             -> ipow(a, 7, c)              
-> ipow(a, 3, c)              -> ipow(a, 1, c)              
-> ipow(a, 0, c)

굉장히 복잡하고 계산량이 많아 보이지만, 총 31번의 ipow() 함수 호출 각각에는 반복문이 하나도 없기 때문에 시간이 거의 걸리지 않습니다.
매 호출마다 b가 대강 절반이 되기 때문에 $\log_2 b$번의 호출이 일어납니다. 시간복잡도는 $\Theta (\log_2 b)$입니다. 처음의 코드는 $\Theta (b)$였죠.

그런데 재귀 호출은 콜 스택(call stack)을 관리해야 하기 때문에 조금 느리다는 느낌이 들기도 합니다. 반복문으로 바꿀 수 있을까요?
일단 중간중간 ipow() 함수의 반환값 사이에 의존 관계가 있고, 그를 바탕으로 반복문을 세울 수는 있습니다.
그런데 이 원리를 이용하면 조금 어려운 것 같아요. 다른 방법을 소개해 볼게요.

---

b를 절반으로 쪼개는 아이디어는 동일하지만, 수식을 이렇게 묶어볼 수 있습니다.

$$ a^{2s} = \left(a^2\right)^s \ a^{2s+1} = \left(a^2\right)^s \cdot a$$

아까랑 똑같은 식 같아 보인다면 다시 한 번 읽어보셔요.
여기서는 a²를 s승하고 있고, 위에서는 a의 s승을 제곱했었습니다.

이걸 구현한 코드는 이렇게 될 거예요.

int ipow(int a, int b, int c){
    if (b == 0) return 1;
    long long ans = ipow((long long)a * a % c, b/2, c);
    if (b % 2 == 1) {
        ans *= a;
        ans %= c;
    }
    return ans;
}

가장 마지막의 ipow() 호출에서 반환값은 1이죠?
그리고 한 ipow() 함수가 다른 ipow() 함수의 값을 받아온 다음, a를 곱하거나 말거나 하고, 다시 넘겨주게 됩니다.
즉 하나의 반환값을 계속 서로 넘겨주며, 각자 거기에 자신의 a를 곱하거나 말거나 하는 셈입니다.
그럼 어떤 값이 곱해지는지를 고민해 보면 반복문으로 고치는 힌트를 얻을 수 있습니다.

int ipow(int a, int b, int c){
    long long ans = 1;
    while (b) {
        if (b % 2 == 1) {
            ans *= a;
            ans %= c;
        }
        a = (long long)a * a % c;
        b /= 2;
    }
    return ans;
}

이 글에서 사용한 "절반으로 쪼개는 원리"는 이밖에도 쓰임새가 많습니다.