Original

인수분해를 이용한 이진 거듭제곱

정수 $a$, $x$, $y$와 $d \geq 3$이 주어지고 $x$가 홀수일 때, $ax^y \pmod{2^d}$를 계산하는 문제를 생각해 봅시다.

아래 알고리즘을 사용하면 $O(d)$번의 덧셈과 이진 연산, 그리고 $y$에 의한 곱셈 한 번으로 이 문제를 풀 수 있습니다.

$2^d$를 법으로 하는 곱셈군의 구조 때문에, $x \equiv 1 \pmod 4$를 만족하는 임의의 수 $x$는 다음과 같이 나타낼 수 있습니다.

$$ x \equiv b^{L(x)} \pmod{2^d}, $$

여기서 $b \equiv 5 \pmod 8$입니다. $x \mapsto -x$ 및 $a \mapsto (-1)^{y} a$로 치환하면 $x \equiv 3 \pmod 4$인 경우를 $x \equiv 1 \pmod 4$인 경우로 바꿀 수 있으므로, 일반성을 잃지 않고 $x \equiv 1 \pmod 4$라고 가정합니다. 이 표기법에서 $ax^y$는 다음과 같이 나타납니다.

$$ a x^y \equiv a b^{yL(x)} \pmod{2^d}. $$

이 알고리즘의 핵심 아이디어는 $2^d$를 법으로 계산한다는 사실을 이용하여 $L(x)$와 $b^{y L(x)}$의 계산을 단순화하는 것입니다. 이유는 나중에 분명해지겠지만, $L(x)$ 대신 $4L(x)$를 다루되 $2^{d-2}$가 아닌 $2^d$를 법으로 하여 취합니다.

이 글에서는 $32$비트 정수에 대한 구현을 다룹니다. 다음과 같이 둡시다.

• mbin_log_32(r, x)는 $r+4L(x) \pmod{2^d}$를 계산하는 함수입니다;
• mbin_exp_32(r, x)는 $r b^{\frac{x}{4}} \pmod{2^d}$를 계산하는 함수입니다;
• mbin_power_odd_32(a, x, y)는 $ax^y \pmod{2^d}$를 계산하는 함수입니다.

그러면 mbin_power_odd_32는 다음과 같이 구현됩니다.

uint32_t mbin_power_odd_32(uint32_t rem, uint32_t base, uint32_t exp) {
    if (base & 2) {
        /* divider is considered negative */
        base = -base;
        /* check if result should be negative */
        if (exp & 1) {
            rem = -rem;
        }
    }
    return (mbin_exp_32(rem, mbin_log_32(0, base) * exp));
}

x로부터 4L(x) 계산하기

$x$를 $x \equiv 1 \pmod 4$를 만족하는 홀수라고 합시다. 이는 다음과 같이 나타낼 수 있습니다.

$$ x \equiv (2^{a_1}+1)\dots(2^{a_k}+1) \pmod{2^d}, $$

여기서 $1 < a_1 < \dots < a_k < d$입니다. 각 인수는 $4$를 법으로 $1$과 같으므로, 여기서 $L(\cdot)$은 각 인수에 대해 잘 정의됩니다. 따라서 다음이 성립합니다.

$$ 4L(x) \equiv 4L(2^{a_1}+1)+\dots+4L(2^{a_k}+1) \pmod{2^{d}}. $$

그러므로 모든 $1 < k < d$에 대해 $t_k = 4L(2^n+1)$을 미리 계산해 두면, 임의의 수 $x$에 대해 $4L(x)$를 계산할 수 있습니다.

32비트 정수의 경우 다음 표를 사용할 수 있습니다.

const uint32_t mbin_log_32_table[32] = {
    0x00000000, 0x00000000, 0xd3cfd984, 0x9ee62e18,
    0xe83d9070, 0xb59e81e0, 0xa17407c0, 0xce601f80,
    0xf4807f00, 0xe701fe00, 0xbe07fc00, 0xfc1ff800,
    0xf87ff000, 0xf1ffe000, 0xe7ffc000, 0xdfff8000,
    0xffff0000, 0xfffe0000, 0xfffc0000, 0xfff80000,
    0xfff00000, 0xffe00000, 0xffc00000, 0xff800000,
    0xff000000, 0xfe000000, 0xfc000000, 0xf8000000,
    0xf0000000, 0xe0000000, 0xc0000000, 0x80000000,
};

실제로는 위에서 설명한 것과 약간 다른 방법을 사용합니다. $x$의 인수분해를 찾는 대신, $x$가 $2^d$를 법으로 $1$이 될 때까지 $x$에 $2^n+1$을 차례로 곱합니다. 이렇게 하면 $x^{-1}$의 표현을 찾게 되며, 이는 다음과 같습니다.

$$ x (2^{a_1}+1)\dots(2^{a_k}+1) \equiv 1 \pmod {2^d}. $$

이를 위해 $1 < n < d$를 만족하는 $n$에 대해 반복합니다. 현재 $x$의 $n$번째 비트가 켜져 있으면 $x$에 $2^n+1$을 곱하며, 이는 C++에서 x = x + (x << n)으로 간편하게 수행할 수 있습니다. $x$가 홀수이므로, 이 연산은 $n$보다 낮은 비트는 바꾸지 않지만 $n$번째 비트는 0으로 만듭니다.

이 모든 것을 고려하면, 함수 mbin_log_32(r, x)는 다음과 같이 구현됩니다.

uint32_t mbin_log_32(uint32_t r, uint32_t x) {
    uint8_t n;

    for (n = 2; n < 32; n++) {
        if (x & (1 << n)) {
            x = x + (x << n);
            r -= mbin_log_32_table[n];
        }
    }

    return r;
}

$4L(x) = -4L(x^{-1})$이므로, $4L(2^n+1)$을 더하는 대신 처음에 $0$인 $r$에서 이를 뺀다는 점에 유의하십시오.

4L(x)로부터 x 계산하기

$k \geq 1$에 대해 다음이 성립합니다.

$$ (a 2^{k}+1)^2 = a^2 2^{2k} +a 2^{k+1}+1 = b2^{k+1}+1, $$

이로부터 (제곱을 반복하여) 다음을 유도할 수 있습니다.

$$ (2^a+1)^{2^b} \equiv 1 \pmod{2^{a+b}}. $$

이 결과를 $a=2^n+1$과 $b=d-k$에 적용하면 $2^n+1$의 곱셈적 위수가 $2^{d-n}$의 약수임을 유도할 수 있습니다.

이는 다시 $L(2^n+1)$이 $2^{n}$으로 나누어떨어져야 함을 의미합니다. $b$의 위수는 $2^{d-2}$이고 $b^y$의 위수는 $2^{d-2-v}$이기 때문입니다(여기서 $2^v$는 $y$를 나누는 $2$의 가장 큰 거듭제곱입니다). 따라서 다음이 필요합니다.

$$ 2^{d-k} \equiv 0 \pmod{2^{d-2-v}}, $$

따라서 $v$는 $k-2$ 이상이어야 합니다. 이는 다소 깔끔하지 않으며, 이를 완화하기 위해 처음에 $L(x)$에 $4$를 곱한다고 언급했습니다. 이제 $4L(x)$를 알고 있다면, $4L(x)$의 비트를 차례로 확인하여 이를 $4L(2^n+1)$들의 합으로 유일하게 분해할 수 있습니다. $n$번째 비트가 $1$이면, 결과에 $2^n+1$을 곱하고 현재의 $4L(x)$에서 $4L(2^n+1)$을 뺍니다.

따라서 mbin_exp_32는 다음과 같이 구현됩니다.

uint32_t mbin_exp_32(uint32_t r, uint32_t x) {
    uint8_t n;

    for (n = 2; n < 32; n++) {
        if (x & (1 << n)) {
            r = r + (r << n);
            x -= mbin_log_32_table[n];
        }
    }

    return r;
}

추가 최적화

$4L(2^{d-1}+1)=2^{d-1}$이고 $2k \geq d$에 대해 다음이 성립한다는 점에 주목하면 반복 횟수를 절반으로 줄일 수 있습니다.

$$ (2^n+1)^2 \equiv 2^{2n} + 2^{n+1}+1 \equiv 2^{n+1}+1 \pmod{2^d}, $$

이로부터 $2n \geq d$일 때 $4L(2^n+1)=2^n$임을 유도할 수 있습니다. 따라서 $\frac{d}{2}$까지만 진행한 뒤 위의 사실을 이용하여 나머지 부분을 비트 연산으로 계산하도록 알고리즘을 단순화할 수 있습니다.

uint32_t mbin_log_32(uint32_t r, uint32_t x) {
    uint8_t n;

    for (n = 2; n != 16; n++) {
        if (x & (1 << n)) {
            x = x + (x << n);
            r -= mbin_log_32_table[n];
        }
    }

    r -= (x & 0xFFFF0000);

    return r;
}

uint32_t mbin_exp_32(uint32_t r, uint32_t x) {
    uint8_t n;

    for (n = 2; n != 16; n++) {
        if (x & (1 << n)) {
            r = r + (r << n);
            x -= mbin_log_32_table[n];
        }
    }

    r *= 1 - (x & 0xFFFF0000);

    return r;
}

로그 테이블 계산

로그 테이블을 계산하려면 법이 $2$의 거듭제곱인 경우에 맞게 Pohlig–Hellman 알고리즘을 수정하면 됩니다.

여기서의 주된 과제는 $g^x \equiv y \pmod{2^d}$를 만족하는 $x$를 계산하는 것입니다. 여기서 $g=5$이고 $y$는 $2^n+1$ 꼴의 수입니다.

양변을 $k$번 제곱하면 다음을 얻습니다.

$$ g^{2^k x} \equiv y^{2^k} \pmod{2^d}. $$

$g$의 위수는 $2^{d}$보다 크지 않으므로(실제로는 $2^{d-2}$보다 크지 않지만, 편의상 $2^d$를 사용하겠습니다), $k=d-1$을 사용하면 좌변에 $g^1$ 또는 $g^0$이 남게 되며, 이를 통해 $y^{2^k}$를 $g$와 비교하여 $x$의 최하위 비트를 결정할 수 있습니다. 이제 $x=x_0 + 2^k x_1$이라고 가정합시다. 여기서 $x_0$은 이미 알고 있는 부분이고 $x_1$은 아직 알지 못하는 부분입니다. 그러면

$$ g^{x_0+2^k x_1} \equiv y \pmod{2^d}. $$

양변에 $g^{-x_0}$를 곱하면 다음을 얻습니다.

$$ g^{2^k x_1} \equiv (g^{-x_0} y) \pmod{2^d}. $$

이제 양변을 $d-k-1$번 제곱하면 $x$의 다음 비트를 얻을 수 있으며, 결국 모든 비트를 복원할 수 있습니다.

참고 문헌

• M30, Hans Petter Selasky, 2009