pkglog
EN

C++ 수학 알고리즘: 에라토스테네스의 체, 유클리드 GCD, 모듈러 거듭제곱, 행렬, FFT

작성 · · 35분 읽기 · 수정 고급 실습

들어가며: “분수 계산이 오버플로우로 틀려요”

분수를 더하거나 큰 수의 거듭제곱을 모듈러로 구할 때 곱셈과 나눗셈을 그대로 쓰면 오버플로우와 시간 초과에 빠지기 쉽습니다. 핵심은 “a×b를 다 계산한 뒤 mod를 취하는 것”과 “중간 결과마다 mod를 취하고, 그 곱셈마저 넘치지 않게 하는 것”의 차이입니다.

flowchart TD
  subgraph wrong[❌ 단순 곱셈 후 mod]
    W1[a × b % mod] --> W2[a×b 오버플로우]
    W2 --> W3[잘못된 결과]
    W3 --> W4[디버깅 어려움]
  end
  subgraph right[✅ 모듈러 곱셈 활용]
    R1["(a % mod) × (b % mod) % mod"] --> R2[또는 128비트 활용]
    R2 --> R3[안전한 계산]
    R3 --> R4[암호화·코딩테스트 통과]
  end

이 글은 실제로 자주 막히는 상황부터 시작해 GCD·LCM, 소수 판별과 에라토스테네스의 체, 모듈러 연산과 빠른 거듭제곱, 행렬 거듭제곱, FFT 구현을 차례로 보고, 이어서 자주 하는 실수와 최적화 방법을 정리합니다.

오버플로우와 시간 초과가 나는 수학 문제들

분수 계산 오버플로우

a/b + c/d를 (a·d + c·b) / (b·d)로 계산하면, 분모와 분자가 10^9 수준만 돼도 곱이 10^18을 넘어 int64_t가 위험해지고, 몇 번 더하면 확실히 넘칩니다. 곱하기 전에 gcd(b, d)로 나눠 공통 인수를 없애고, 결과도 매번 기약분수로 줄여 두면 수가 불필요하게 커지지 않습니다. 그래도 범위가 모자라면 __int128로 중간값을 계산합니다.

소수 판별이 너무 느림

1부터 10^6까지의 소수 개수를 세는데 수마다 O(√n) 판별을 하면 전체가 O(n√n)이라 시간 초과가 나기 쉽습니다. 에라토스테네스의 체를 쓰면 O(n log log n)에 범위 안의 모든 소수를 한 번에 구합니다.

거듭제곱 mod 연산 시간 초과

a^b % mod에서 b가 10^18이면 b번 곱하는 루프는 끝나지 않습니다. pow(a, b) % mod는 double로 계산하므로 정밀도가 모자라고 정수 결과를 보장하지도 않습니다. 빠른 거듭제곱(분할 정복)을 쓰면 O(log b)에 구할 수 있습니다.

행렬 거듭제곱으로 피보나치

n이 10^18이면 dp[i] = dp[i-1] + dp[i-2]를 n번 반복할 수 없습니다. [fib(n+1), fib(n)]^T = M^n × [1, 0]^T 관계를 이용해 행렬 M의 거듭제곱을 O(log n)에 구하면 됩니다.

다항식 곱셈이 O(n²)로 느림

차수가 10^5인 두 다항식을 이중 루프로 곱하면 10^10번의 연산이 필요합니다. FFT(고속 푸리에 변환)를 쓰면 O(n log n)에 곱할 수 있습니다.

모듈러에서 나눗셈이 필요할 때

(a / b) % mod는 (a % mod) / (b % mod)로 계산할 수 없습니다. 모듈러 세계에는 일반적인 나눗셈이 없기 때문입니다. 대신 b의 모듈러 역원 b^(-1)을 구해 a × b^(-1) % mod로 계산합니다. mod가 소수라면 페르마 소정리로, 아니라면 확장 유클리드로 역원을 구합니다.

알고리즘 선택 가이드

문제 유형추천 기법시간 복잡도
최대공약수/최소공배수유클리드 알고리즘O(log min(a,b))
소수 판별 (단일)trial divisionO(√n)
소수 목록 (범위)에라토스테네스 체O(n log log n)
a^b % mod빠른 거듭제곱O(log b)
모듈러 역원확장 유클리드 / 페르마O(log mod)
행렬 거듭제곱분할 정복O(k³ log n)
다항식 곱셈FFTO(n log n)

수학 알고리즘의 핵심 원리와 GCD·LCM 관계

flowchart LR
  subgraph core[핵심 원리]
    A[GCD] --> A1["유클리드: gcd(a,b)=gcd(b,a%b)"]
    B[소수] --> B1[에라토스테네스: 배수 제거]
    C[모듈러] --> C1["(a×b)%m = ((a%m)×(b%m))%m"]
    D[행렬] --> D1["M^n = M^(n/2)×M^(n/2)"]
  end

GCD와 LCM 사이에는 gcd(a, b) × lcm(a, b) = a × b라는 관계가 있습니다. 그래서 lcm = a × b / gcd로 구할 수 있지만, a × b를 먼저 계산하면 오버플로우가 나므로 a / gcd × b 순서로 계산합니다.

// gcd(a, b) × lcm(a, b) = a × b
// 따라서 lcm(a, b) = a × b / gcd(a, b)
// 주의: a*b에서 오버플로우 가능 → lcm = a / gcd(a,b) * b

GCD·소수·모듈러 거듭제곱·행렬·FFT 구현

GCD (최대공약수) - 유클리드 알고리즘

gcd(a, b) = gcd(b, a % b)를 b가 0이 될 때까지 반복하면 남은 a가 최대공약수입니다. 나머지가 매 두 단계마다 절반 이하로 줄어들기 때문에 반복 횟수는 O(log min(a, b))입니다.

#include <cstdint>
// 반복문 버전 (권장: 스택 오버플로우 없음)
// 시간: O(log min(a,b))
int64_t gcd_iter(int64_t a, int64_t b) {
    a = std::abs(a);
    b = std::abs(b);
    while (b != 0) {
        int64_t t = b;
        b = a % b;
        a = t;
    }
    return a;
}
// 재귀 버전
int64_t gcd_rec(int64_t a, int64_t b) {
    return b == 0 ? std::abs(a) : gcd_rec(b, a % b);
}
// C++17: std::gcd, std::lcm 사용 가능
#include <numeric>
// int64_t g = std::gcd(a, b);
// int64_t l = std::lcm(a, b);

재귀 깊이도 O(log) 수준이라 재귀 버전이 실제로 스택을 넘치게 할 일은 거의 없지만, 반복문 버전이 디버깅과 인라인 최적화 면에서 조금 더 다루기 쉽습니다. C++17 이상이라면 <numeric>의 std::gcd를 쓰면 됩니다.

LCM (최소공배수)

// 오버플로우 방지: 먼저 나누고 곱하기
// lcm(a,b) = a * (b / gcd(a,b))
int64_t lcm_safe(int64_t a, int64_t b) {
    if (a == 0 || b == 0) return 0;
    int64_t g = gcd_iter(a, b);
    // a/g * b 는 g로 나눠지므로 정수. 순서 중요!
    return (a / g) * b;  // a*b/g와 값은 같지만 중간값이 작음
}
// 여러 수의 LCM
int64_t lcm_range(const std::vector<int64_t>& arr) {
    int64_t result = 1;
    for (int64_t x : arr) {
        result = lcm_safe(result, x);
        if (result == 0) break;
    }
    return result;
}

나누기를 먼저 하면 중간값이 줄어들 뿐, 최종 LCM 자체가 64비트를 넘는다면 여전히 오버플로우가 납니다. 여러 수의 LCM은 금방 커지므로, 범위를 넘을 수 있다면 아래 “에러 처리와 검증”의 검사를 함께 둡니다.

확장 유클리드 - 베주 항등식과 모듈러 역원

확장 유클리드는 ax + by = gcd(a, b)를 만족하는 정수 x, y를 함께 구합니다. gcd(a, m) = 1이면 ax + my = 1이 되고, 양변을 m으로 나눈 나머지를 보면 ax ≡ 1 (mod m)이므로 x가 곧 a의 역원입니다. mod가 소수가 아니어도 쓸 수 있다는 점이 페르마 방식과 다릅니다.

// 확장 유클리드: ax + by = gcd(a,b) 의 x, y 반환
// 반환값: gcd(a, b)
int64_t ext_gcd(int64_t a, int64_t b, int64_t& x, int64_t& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    int64_t x1, y1;
    int64_t g = ext_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}
// 모듈러 역원: a^(-1) mod m (gcd(a,m)=1 일 때)
// ax ≡ 1 (mod m) 인 x
int64_t mod_inverse(int64_t a, int64_t m) {
    int64_t x, y;
    int64_t g = ext_gcd((a % m + m) % m, m, x, y);
    if (g != 1) return -1;  // 역원 없음
    return (x % m + m) % m;
}

소수 판별 (단일 수)

#include <cmath>
// Trial division: O(√n)
bool is_prime_trial(int64_t n) {
    if (n < 2) return false;
    if (n == 2 || n == 3) return true;
    if (n % 2 == 0 || n % 3 == 0) return false;
    for (int64_t i = 5; i * i <= n; i += 6) {
        if (n % i == 0 || n % (i + 2) == 0)
            return false;
    }
    return true;
}
// 6k±1 최적화: 2,3 배수 건너뛰기

2와 3보다 큰 소수는 모두 6k±1 꼴이므로, 5부터 6씩 건너뛰며 i와 i+2만 검사하면 검사 횟수가 3분의 1로 줄어듭니다. 복잡도는 여전히 O(√n)이라 n이 10^18에 가까우면 이 방식으로도 느립니다. 그럴 때는 Miller-Rabin 소수 판별을 씁니다.

에라토스테네스의 체

2부터 시작해 소수를 하나 찾을 때마다 그 배수를 모두 지워 나갑니다. 끝까지 지워지지 않은 수가 소수입니다.

#include <vector>
// 에라토스테네스의 체: O(n log log n)
std::vector<int> sieve_of_eratosthenes(int n) {
    std::vector<bool> is_prime(n + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; i * i <= n; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j <= n; j += i) {
            is_prime[j] = false;
        }
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[i]) primes.push_back(i);
    }
    return primes;
}
// 소수 개수만 세기 (목록을 만들지 않음)
int count_primes(int n) {
    if (n < 2) return 0;
    std::vector<bool> is_prime(n + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; i * i <= n; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j <= n; j += i)
            is_prime[j] = false;
    }
    int count = 0;
    for (int i = 2; i <= n; ++i)
        if (is_prime[i]) ++count;
    return count;
}

바깥 루프가 i * i <= n까지만 도는 이유는, √n보다 큰 합성수는 반드시 √n 이하의 소인수를 가져서 이미 지워졌기 때문입니다.

모듈러 연산과 빠른 거듭제곱

a^b = (a^(b/2))^2 × a^(b%2)를 이용해 b를 이진수로 보고 비트마다 제곱과 곱셈을 한 번씩 합니다.

// 빠른 거듭제곱 (분할 정복)
// a^b = (a^(b/2))^2 * a^(b%2)
// 주의: result * a, a * a가 int64_t에 들어가려면 mod < 약 3e9 이어야 함
int64_t mod_pow(int64_t a, int64_t b, int64_t mod) {
    if (mod == 1) return 0;
    int64_t result = 1;
    a %= mod;
    while (b > 0) {
        if (b & 1) {
            result = (result * a) % mod;
        }
        a = (a * a) % mod;
        b >>= 1;
    }
    return result;
}
// 모듈러 곱셈 (오버플로우 방지)
// a * b % mod 를 64비트만으로 안전하게 (mod < 2^62)
int64_t mod_mul(int64_t a, int64_t b, int64_t mod) {
    a %= mod;
    b %= mod;
    int64_t result = 0;
    while (b > 0) {
        if (b & 1) {
            result = (result + a) % mod;
        }
        a = (a * 2) % mod;
        b >>= 1;
    }
    return result;
}
// __int128 사용 (GCC/Clang)
int64_t mod_mul_128(int64_t a, int64_t b, int64_t mod) {
    return static_cast<int64_t>((__int128)a * b % mod);
}

mod_pow는 곱셈 결과가 int64_t에 들어간다는 가정 위에 있습니다. 코딩 테스트에서 흔한 1'000'000'007 정도면 문제가 없지만, mod가 10^18 수준이면 곱셈을 mod_mul이나 mod_mul_128로 바꿔야 합니다. mod_mul은 곱셈을 O(log b)번의 덧셈으로 바꾸므로 느리고, __int128을 쓸 수 있다면 그쪽이 훨씬 빠릅니다. MSVC에는 __int128이 없으므로 _umul128 같은 내장 함수나 mod_mul을 씁니다.

페르마 소정리를 이용한 모듈러 역원

p가 소수이고 a가 p의 배수가 아니면 a^(p-1) ≡ 1 (mod p)입니다. 양변에 a^(-1)을 곱하면 a^(-1) ≡ a^(p-2) (mod p)가 되므로 빠른 거듭제곱 한 번으로 역원을 구할 수 있습니다.

// mod가 소수일 때: a^(-1) = a^(p-2) mod p
int64_t mod_inverse_fermat(int64_t a, int64_t p) {
    return mod_pow(a, p - 2, p);
}
// nCr % p (p 소수): 팩토리얼 역원 전처리
// nCr = n! / (r! * (n-r)!) → n! * inv(r!) * inv((n-r)!)

a가 p의 배수이면 이 식은 0을 돌려주는데, 이때는 역원이 존재하지 않는 것이므로 호출하는 쪽에서 걸러야 합니다.

행렬 연산과 행렬 거듭제곱

선형 점화식은 행렬 곱으로 표현할 수 있고, 행렬도 빠른 거듭제곱을 그대로 적용할 수 있습니다. k×k 행렬 곱셈이 O(k³)이므로 M^n은 O(k³ log n)입니다.

#include <vector>
using Matrix = std::vector<std::vector<int64_t>>;
const int64_t MOD = 1'000'000'007;
// 행렬 곱셈: A * B
Matrix mat_mul(const Matrix& A, const Matrix& B) {
    int n = static_cast<int>(A.size());
    int m = static_cast<int>(A[0].size());
    int p = static_cast<int>(B[0].size());
    Matrix C(n, std::vector<int64_t>(p, 0));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < p; ++j) {
            for (int k = 0; k < m; ++k) {
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
            }
        }
    }
    return C;
}
// 행렬 거듭제곱: M^n
Matrix mat_pow(Matrix M, int64_t n) {
    int k = static_cast<int>(M.size());
    Matrix result(k, std::vector<int64_t>(k, 0));
    for (int i = 0; i < k; ++i) result[i][i] = 1;  // 단위행렬
    while (n > 0) {
        if (n & 1) result = mat_mul(result, M);
        M = mat_mul(M, M);
        n >>= 1;
    }
    return result;
}
// 피보나치: fib(n) = M^(n-1)[0][0], M = [[1,1],[1,0]]
int64_t fib_matrix(int64_t n) {
    if (n <= 1) return n;
    Matrix M = {{1, 1}, {1, 0}};
    Matrix Mn = mat_pow(M, n - 1);
    return Mn[0][0];
}

M = [[1,1],[1,0]]의 k제곱은 [[F(k+1), F(k)], [F(k), F(k-1)]]이 되므로, M^(n-1)의 왼쪽 위 원소가 F(n)입니다.

소인수 분해

// Trial division으로 소인수 분해
std::vector<std::pair<int64_t, int>> factorize(int64_t n) {
    std::vector<std::pair<int64_t, int>> factors;
    for (int64_t i = 2; i * i <= n; ++i) {
        if (n % i != 0) continue;
        int cnt = 0;
        while (n % i == 0) {
            n /= i;
            ++cnt;
        }
        factors.emplace_back(i, cnt);
    }
    if (n > 1) factors.emplace_back(n, 1);
    return factors;
}
// 약수 개수: (e1+1)*(e2+1)*... (각 소인수 지수+1의 곱)
int count_divisors(int64_t n) {
    auto factors = factorize(n);
    int result = 1;
    for (const auto& [p, e] : factors) result *= (e + 1);
    return result;
}

루프가 끝났을 때 n이 1보다 크면, 남은 n은 √(원래 n)보다 큰 소인수 하나입니다. 같은 범위의 수를 여러 번 분해해야 한다면, 체를 돌리면서 각 수의 가장 작은 소인수를 기록해 두면 분해가 O(log n)으로 빨라집니다.

FFT (고속 푸리에 변환)

다항식을 계수 대신 n개 점에서의 값으로 바꾸면, 곱셈은 점마다 값을 곱하기만 하면 됩니다. FFT는 이 변환과 역변환을 1의 n제곱근에서 분할 정복으로 O(n log n)에 수행합니다. 아래는 Cooley-Tukey 방식의 반복 구현입니다.

// FFT는 복소수 연산과 분할 정복을 이용
// Cooley-Tukey 알고리즘이 대표적
// 실전에서는 검증된 구현(kactl, AtCoder Library 등) 사용 권장
#include <complex>
#include <vector>
using Complex = std::complex<double>;
const double PI = 3.14159265358979323846;
void fft(std::vector<Complex>& a, bool invert) {
    int n = static_cast<int>(a.size());
    for (int i = 1, j = 0; i < n; ++i) {
        int bit = n >> 1;
        for (; j >= bit; bit >>= 1) j -= bit;
        j += bit;
        if (i < j) std::swap(a[i], a[j]);
    }
    for (int len = 2; len <= n; len <<= 1) {
        double angle = 2 * PI / len * (invert ? -1 : 1);
        Complex wlen(std::cos(angle), std::sin(angle));
        for (int i = 0; i < n; i += len) {
            Complex w(1);
            for (int j = 0; j < len / 2; ++j) {
                Complex u = a[i + j], v = a[i + j + len/2] * w;
                a[i + j] = u + v;
                a[i + j + len/2] = u - v;
                w *= wlen;
            }
        }
    }
    if (invert) {
        for (int i = 0; i < n; ++i) a[i] /= n;
    }
}
std::vector<int64_t> multiply_poly(const std::vector<int64_t>& a,
                                   const std::vector<int64_t>& b) {
    int n = 1;
    while (n < static_cast<int>(a.size()) + static_cast<int>(b.size())) n <<= 1;
    std::vector<Complex> fa(n), fb(n);
    for (size_t i = 0; i < a.size(); ++i) fa[i] = a[i];
    for (size_t i = 0; i < b.size(); ++i) fb[i] = b[i];
    fft(fa, false);
    fft(fb, false);
    for (int i = 0; i < n; ++i) fa[i] *= fb[i];
    fft(fa, true);
    std::vector<int64_t> result(n);
    for (int i = 0; i < n; ++i) {
        result[i] = static_cast<int64_t>(std::round(fa[i].real()));
    }
    return result;
}

double 기반 FFT는 결과 계수가 커질수록 반올림 오차가 커집니다. 결과 계수가 대략 10^14~10^15를 넘으면 std::round로도 정확한 정수를 얻지 못할 수 있으므로, 계수가 크거나 결과를 소수 mod로 구해야 한다면 NTT(정수 기반 수론 변환)나 계수를 쪼개 여러 번 FFT하는 방식을 씁니다.

모듈러 오버플로우·음수 mod·역원 부재 처리

모듈러 연산에서 오버플로우

(a * b) % mod가 틀린 값을 내는 경우입니다. a * b가 int64_t 범위를 넘으면 부호 있는 정수 오버플로우라 정의되지 않은 동작이 되고, 대개는 엉뚱한 값이 나옵니다.

// ❌ 잘못된 코드
int64_t bad = (a * b) % mod;  // a*b 오버플로우
// ✅ 해결 1: mod_mul 사용
int64_t good = mod_mul(a, b, mod);
// ✅ 해결 2: __int128 (GCC/Clang)
int64_t good2 = static_cast<int64_t>((__int128)a * b % mod);

LCM 계산 시 오버플로우

lcm(a, b) = a * b / gcd(a, b)를 그대로 쓰면 a * b가 먼저 계산되어 넘칩니다.

// ❌ 잘못된 코드
int64_t bad = a * b / gcd(a, b);  // a*b 오버플로우
// ✅ 올바른 코드: 먼저 나누고 곱하기
int64_t good = (a / gcd(a, b)) * b;

음수 mod 연산

C++에서 -5 % 3은 1이 아니라 -2입니다. 나눗셈이 0 방향으로 버림되기 때문에 나머지는 피제수의 부호를 따릅니다. 해시 테이블 인덱스나 원형 버퍼 위치처럼 0 이상이어야 하는 곳에 그대로 쓰면 음수 인덱스 버그가 됩니다.

// ❌ 잘못된 가정
int r = -5 % 3;  // r = -2 (C++)
// ✅ 올바른 모듈러 결과 (0 ~ mod-1)
int64_t mod_positive(int64_t a, int64_t m) {
    return (a % m + m) % m;
}

소수 판별에서 i*i 오버플로우

int64_t에서 i가 약 3.04×10^9를 넘으면 i * i가 범위를 넘습니다. n이 10^18에 가까우면 루프 조건 자체가 넘칠 수 있습니다.

// ❌ 잘못된 코드 (n이 10^18일 때)
for (int64_t i = 2; i * i <= n; ++i)  // i*i 오버플로우
// ✅ 올바른 코드
for (int64_t i = 2; i <= n / i; ++i)
// 또는
for (int64_t i = 2; i <= static_cast<int64_t>(std::sqrt(n)); ++i)

std::sqrt는 double로 계산하므로 n이 2^53을 넘으면 결과가 실제 제곱근보다 1 작거나 클 수 있습니다. 큰 n에서는 i <= n / i 쪽이 안전합니다.

모듈러 역원이 존재하지 않을 때

역원은 gcd(a, m) = 1일 때만 존재합니다. 위의 mod_inverse는 이 경우 -1을 돌려주는데, 반환값을 확인하지 않고 그대로 곱하면 조용히 틀린 결과가 나옵니다. 확장 유클리드든 페르마든 역원이 없다는 사실은 바뀌지 않으므로, 이 경우에는 문제를 다르게 풀어야 합니다. 예를 들어 mod를 소인수로 나눠 중국인의 나머지 정리로 합치거나, 나눗셈을 피하도록 식을 바꿉니다.

// ✅ 역원 계산 전 검증
int64_t inv = mod_inverse(a, m);
if (inv == -1) {
    // 역원 없음: 나눗셈을 쓰지 않는 방식으로 풀어야 함
    return ERROR;
}

에라토스테네스 체에서 i*i 시작점

안쪽 루프를 j = i * 2부터 시작해도 결과는 맞지만, 2i, 3i, ..., (i-1)i는 이미 더 작은 소수의 배수로 지워졌으므로 헛걸음입니다. i * i부터 시작하면 됩니다.

// ❌ 비효율 (동작은 함)
for (int j = i * 2; j <= n; j += i)
// ✅ 올바른 최적화: i*i부터 시작
for (int j = i * i; j <= n; j += i)

빠른 거듭제곱에서 mod 1

모든 정수를 1로 나눈 나머지는 0이므로 a^b mod 1은 항상 0이어야 합니다. 그런데 b가 0이면 루프를 한 번도 돌지 않아 초기값 1이 그대로 반환됩니다. 그래서 mod == 1을 먼저 처리합니다.

// ✅ mod 1 처리
int64_t mod_pow(int64_t a, int64_t b, int64_t mod) {
    if (mod == 1) return 0;  // 모든 수 mod 1 = 0
    // ...
}

GCD에 0 입력

gcd_iter(0, 0)은 루프를 돌지 않고 0을 반환하며, std::gcd(0, 0)도 0입니다. 문제는 그 뒤입니다. 반환된 0으로 나누면(기약분수화, LCM 등) 0으로 나누기가 되어 프로그램이 죽습니다. 0이 들어올 수 있는 경로라면 나누기 전에 확인합니다.

// ✅ 0 처리
int64_t gcd_safe(int64_t a, int64_t b) {
    if (a == 0 && b == 0) return 0;  // 호출자가 0으로 나누지 않도록 확인할 것
    return gcd_iter(a, b);
}

행렬 곱셈 차원 불일치

행렬 곱셈은 A의 열 수와 B의 행 수가 같아야 합니다. 위의 mat_mul은 이를 확인하지 않으므로, 크기가 맞지 않으면 범위 밖 접근이 일어납니다.

// ✅ 차원 검증
Matrix mat_mul_safe(const Matrix& A, const Matrix& B) {
    if (A.empty() || B.empty() || A[0].size() != B.size())
        throw std::invalid_argument("Invalid matrix dimensions");
    // ...
}

체 최적화·비트 압축·캐시 친화적 행렬 곱셈

에라토스테네스 체 최적화

짝수는 2를 빼면 모두 합성수이므로 홀수만 검사합니다. 홀수 소수 i의 배수 중 짝수 배수는 이미 지워졌으므로 2 * i씩 건너뛰면 됩니다.

// 짝수 최적화: 2만 따로 처리, 홀수만 체크
std::vector<int> sieve_optimized(int n) {
    if (n < 2) return {};
    std::vector<bool> is_prime(n + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 4; i <= n; i += 2) is_prime[i] = false;
    for (int i = 3; i * i <= n; i += 2) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j <= n; j += 2 * i)  // 홀수 배수만
            is_prime[j] = false;
    }
    std::vector<int> primes = {2};
    for (int i = 3; i <= n; i += 2)
        if (is_prime[i]) primes.push_back(i);
    return primes;
}

메모리: 비트 단위 압축

std::vector<bool>은 표준이 허용하는 특수화로 이미 원소 하나를 1비트로 저장합니다. 그래서 vector<char>나 bool 배열에 비해 메모리가 8분의 1이고, 아래처럼 직접 비트를 다루는 코드와 메모리 사용량은 같습니다. 직접 구현하는 이유는 메모리보다는 워드 단위 연산이나 홀수만 저장하는 압축 같은 추가 최적화를 넣기 위해서입니다.

// uint8_t 하나에 8개 수의 소수 여부를 저장
std::vector<int> sieve_bits(int n) {
    std::vector<uint8_t> is_prime((n + 8) / 8, 0xFF);
    is_prime[0] &= ~3;  // 0, 1 제거
    for (int i = 2; i * i <= n; ++i) {
        if (!(is_prime[i / 8] & (1 << (i % 8)))) continue;
        for (int j = i * i; j <= n; j += i)
            is_prime[j / 8] &= ~(1 << (j % 8));
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i)
        if (is_prime[i / 8] & (1 << (i % 8)))
            primes.push_back(i);
    return primes;
}

n이 수억 이상이라 캐시에 다 들어가지 않는다면, 범위를 캐시 크기 블록으로 나눠 처리하는 구간 체(segmented sieve)가 효과가 큽니다.

행렬 곱셈 캐시 친화적 순서

vector<vector<>>는 행 단위로 메모리가 이어져 있습니다. i-j-k 순서에서 가장 안쪽 루프가 k를 돌면 B[k][j]는 행을 건너뛰며 열 방향으로 접근하게 되어 캐시 효율이 나쁩니다. i-k-j 순서로 바꾸면 안쪽 루프에서 B[k][j]와 C[i][j] 모두 같은 행을 차례로 읽고 씁니다.

// i-k-j 순서: 안쪽 루프가 B와 C의 같은 행을 연속 접근
for (int i = 0; i < n; ++i) {
    for (int k = 0; k < m; ++k) {
        int64_t aik = A[i][k];
        for (int j = 0; j < p; ++j) {
            C[i][j] = (C[i][j] + aik * B[k][j]) % MOD;
        }
    }
}

MOD가 10^9 수준이면 원소 곱 하나가 약 10^18이라, 곱을 여러 개 더한 뒤 한 번에 mod를 취하면 int64_t를 넘습니다. 누적할 때마다 mod를 취하거나, 부호 없는 64비트로 몇 개(대략 16개 이하)를 더할 때마다 줄이는 식으로 범위를 관리해야 합니다.

컴파일 타임 상수 활용

mod가 컴파일 타임 상수이면 컴파일러가 나눗셈을 곱셈과 시프트로 바꾸는 최적화를 적용할 수 있어, 런타임 변수로 나누는 것보다 빠릅니다. 템플릿 인자로 mod를 받는 래퍼가 흔히 쓰이는 이유입니다.

// mod가 컴파일 타임 상수면 컴파일러 최적화 유리
template <int64_t Mod>
struct ModInt {
    int64_t v;
    ModInt(int64_t x = 0) : v((x % Mod + Mod) % Mod) {}
    ModInt operator*(ModInt o) const {
        return ModInt(static_cast<int64_t>((__int128)v * o.v % Mod));
    }
    // ...
};

성능 비교 요약

연산나이브최적화비고
GCDO(min(a,b))O(log min(a,b))유클리드
소수 n개O(n√n)O(n log log n)에라토스테네스
a^b modO(b)O(log b)빠른 거듭제곱
행렬 M^nO(k³ n)O(k³ log n)분할 정복
다항식 곱O(n²)O(n log n)FFT

모듈러 래퍼·팩토리얼 전처리와 nCr·RSA 예시

모듈러 연산 래퍼 클래스

모듈러 연산을 여러 곳에서 쓰면 mod를 빼먹거나 음수 보정을 잊는 실수가 반드시 나옵니다. 연산자를 오버로딩한 래퍼로 감싸 두면 이런 실수를 타입 수준에서 막을 수 있습니다.

template <int64_t Mod>
struct ModInt {
    int64_t v;
    ModInt(int64_t x = 0) : v((x % Mod + Mod) % Mod) {}
    ModInt& operator+=(ModInt o) { v = (v + o.v) % Mod; return *this; }
    ModInt& operator-=(ModInt o) { v = (v - o.v + Mod) % Mod; return *this; }
    ModInt& operator*=(ModInt o) {
        v = static_cast<int64_t>((__int128)v * o.v % Mod);
        return *this;
    }
    ModInt operator+(ModInt o) const { return ModInt(*this) += o; }
    ModInt operator-(ModInt o) const { return ModInt(*this) -= o; }
    ModInt operator*(ModInt o) const { return ModInt(*this) *= o; }
    ModInt pow(int64_t b) const {
        ModInt r(1), a(*this);
        for (; b; b >>= 1, a *= a) if (b & 1) r *= a;
        return r;
    }
    ModInt inv() const { return pow(Mod - 2); }  // 페르마 (Mod 소수일 때)
};
using Mint = ModInt<1'000'000'007>;

팩토리얼 전처리와 nCr % p

n, r이 10^6 수준인 nCr을 소수 p로 나눈 나머지를 여러 번 구해야 한다면, 팩토리얼과 그 역원을 O(n)에 미리 계산해 두고 쿼리마다 O(1)에 답합니다. 역원을 매번 구하지 않고 inv_fact[n] 하나만 페르마로 구한 뒤, inv_fact[i-1] = inv_fact[i] × i 관계로 거꾸로 채우는 것이 요령입니다.

// nCr, nPr 등에 필요한 팩토리얼·역원 전처리
struct FactorialPrecompute {
    std::vector<Mint> fact, inv_fact;
    FactorialPrecompute(int n) : fact(n + 1), inv_fact(n + 1) {
        fact[0] = 1;
        for (int i = 1; i <= n; ++i) fact[i] = fact[i - 1] * i;
        inv_fact[n] = fact[n].inv();
        for (int i = n; i >= 1; --i) inv_fact[i - 1] = inv_fact[i] * i;
    }
    Mint nCr(int n, int r) {
        if (r < 0 || r > n) return 0;
        return fact[n] * inv_fact[r] * inv_fact[n - r];
    }
};

이 방식은 n이 p보다 작을 때만 성립합니다. n ≥ p이면 n!이 p의 배수가 되어 역원이 없으므로 뤼카 정리 같은 다른 방법이 필요합니다.

에러 처리와 검증

// GCD 기반 검증
int64_t safe_lcm(int64_t a, int64_t b) {
    if (a == 0 || b == 0) return 0;
    int64_t g = gcd_iter(a, b);
    if (std::abs(a / g) > INT64_MAX / std::abs(b)) {
        throw std::overflow_error("LCM overflow");
    }
    return std::abs((a / g) * b);
}

곱하기 전에 INT64_MAX / |b|와 비교하면 곱셈을 하지 않고도 오버플로우 여부를 알 수 있습니다. GCC와 Clang에서는 __builtin_mul_overflow로 같은 검사를 더 간단히 할 수 있습니다.

단위 테스트용 테스트 케이스

// 검증용 예시
void test_math_algorithms() {
    assert(gcd_iter(48, 18) == 6);
    assert(gcd_iter(0, 5) == 5);
    assert(lcm_safe(12, 18) == 36);
    assert(is_prime_trial(97) == true);
    assert(is_prime_trial(100) == false);
    assert(mod_pow(2, 10, 1000) == 24);
    assert(mod_inverse_fermat(3, 7) == 5);  // 3*5=15≡1 (mod 7)
    assert(fib_matrix(10) == 55);
}

알고리즘 선택 플로우

flowchart TD
    A[수학 문제] --> B{문제 유형}
    B -->|최대공약수/배수| C[유클리드]
    B -->|소수 범위| D[에라토스테네스]
    B -->|소수 단일| E[Trial division]
    B -->|a^b mod| F[빠른 거듭제곱]
    B -->|나눗셈 mod| G[모듈러 역원]
    B -->|선형 점화식| H[행렬 거듭제곱]
    B -->|다항식 곱| I[FFT]

실전 활용 예시: 기약분수화

std::pair<int64_t, int64_t> to_irreducible(int64_t a, int64_t b) {
    if (b == 0) throw std::invalid_argument("denominator is zero");
    int64_t g = gcd_iter(std::abs(a), std::abs(b));
    a /= g;
    b /= g;
    if (b < 0) { a = -a; b = -b; }  // 분모를 양수로
    return {a, b};
}

b가 0이 아니면 g도 0이 아니므로 나눗셈이 안전합니다. 분모의 부호를 항상 양수로 맞춰 두면 같은 분수가 한 가지 표현만 가지게 되어 비교와 해싱이 쉬워집니다.

실전 활용 예시: RSA 스타일 거듭제곱

RSA 암호화는 m^e mod n 계산입니다. 실제 RSA의 n은 2048비트 이상이라 큰 정수 라이브러리가 필요하고, 여기서는 64비트 범위에서 원리만 보여 줍니다. n이 커도 곱셈이 넘치지 않도록 __int128로 곱합니다.

// 64비트 모듈러에서 안전한 거듭제곱 (원리 확인용, 실제 RSA에는 쓰지 말 것)
int64_t rsa_pow(int64_t m, int64_t e, int64_t n) {
    int64_t result = 1 % n;
    m %= n;
    while (e > 0) {
        if (e & 1) result = static_cast<int64_t>((__int128)result * m % n);
        m = static_cast<int64_t>((__int128)m * m % n);
        e >>= 1;
    }
    return result;
}

실제 암호 구현에서는 입력에 따라 실행 시간이 달라지지 않는 상수 시간 연산과 패딩 처리가 필수이므로, 직접 구현하지 말고 OpenSSL 같은 검증된 라이브러리를 씁니다.

자주 묻는 질문 (FAQ)

Q. C++에서 음수에 % 연산을 하면 왜 음수가 나오나요?

C++11부터 정수 나눗셈은 0 방향으로 버림하도록 정의되어 있어서, -7 % 3의 결과는 2가 아니라 -1입니다. 모듈러 연산 결과를 배열 인덱스나 해시 값으로 쓰면 음수 인덱스로 인한 버그가 생길 수 있습니다. 항상 0 이상의 값이 필요하다면 ((a % m) + m) % m처럼 m을 한 번 더해 보정해야 합니다.

Q. 오버플로우를 피하는 가장 좋은 방법은?

GCC/Clang에서는 __int128로 곱셈을 하거나, mod_mul처럼 곱셈을 덧셈으로 분해하거나, ModInt 래퍼로 연산을 캡슐화하는 방법이 있습니다. 여러 곳에서 모듈러 연산을 쓴다면 래퍼 클래스로 한곳에 모아 두는 편이 실수가 적습니다.

같이 보면 좋은 글