競技プログラミングにおける高度な数値最適化と整数論アルゴリズム

合同算術に基づく最短経路アプローチ

広範な値域における到達可能性の判定や組合せカウント、極値探索に有効な手法です。基準となる法 \(M\) を定め、各剰余類 \(0 \dots M-1\) をグラフの頂点としてモデル化します。状態 \(D[r]\) を、剰余類 \(r\) に属する値を実現するための最小コスト(または最小値)と定義します。状態間の遷移を重み付き有向辺で表現し、ダイクストラ法などの最短経路アルゴリズムを適用することで、各剰余クラスにおける達成可能な下限値を効率的に導出できます。このアプローチは、値が巨大でもモジュロ空間に閉じたグラフ探索に帰着させることで計算量を削減するものです。

高次元累積和(SOS DP)

ビット列で表現された集合について、部分集合または上位集合の属性値の総和を高速に計算する動的計画法です。\(N\) ビットのマスク空間において、各ビットの包含・除外関係を順次処理することで、\(O(N \cdot 2^N)\) の時間で全マスクに対する累積値を求められます。以下の実装は部分集合和を計算する典型例です。

void compute_subset_sums(std::vector<int>& dp, int n) {
    int limit = 1 << n;
    for (int bit = 0; bit < n; ++bit) {
        for (int mask = 0; mask < limit; ++mask) {
            if (mask & (1 << bit)) {
                dp[mask] += dp[mask ^ (1 << bit)];
            }
        }
    }
}

逆に、特定のマスクを包含する上位集合(スーパーセット)の累積値を求める場合は、判定条件と遷移方向を反転させます。

void compute_superset_sums(std::vector<int>& dp, int n) {
    int limit = 1 << n;
    for (int bit = 0; bit < n; ++bit) {
        for (int mask = 0; mask < limit; ++mask) {
            if (!(mask & (1 << bit))) {
                dp[mask] += dp[mask | (1 << bit)];
            }
        }
    }
}

加算演算を最小値・最大値演算に置き換えることで、偏序関係における範囲最適化問題にも直接適用可能です。

0/1 分数最適化

\(N\) 個の候補から任意の部分集合を選び、属性 \(A\) の総和と属性 \(B\) の総和の比率 \(\frac{\sum A_i}{\sum B_i}\) を最大化する問題に対して適用されます。目標の最適値 \(\lambda\) を仮定すると、条件式は \(\sum A_i - \lambda \sum B_i \geq 0\) と変形でき、さらに \(\sum (A_i - \lambda B_i) \geq 0\) となります。これにより、全体最適化問題が各要素のスコア \(C_i(\lambda) = A_i - \lambda B_i\) に基づく単一基準の選択問題に帰着します。実際の解法では、二分探索を用いて \(\lambda\) を狭め、各試行において \(C_i(\lambda)\) の上位 \(K\) 個または正の値を持つ要素を採用することで、条件を満たすか検証します。

三分探索法

単峰性(または単谷性)を持つ関数の極値を数値的に近似する手法です。探索区間 \([L, R]\) を三等分し、内部の二点 \(m_1, m_2\) における関数値 \(f(m_1), f(m_2)\) を比較します。\(f(m_1) < f(m_2)\) の場合、極小値は区間 \([m_1, R]\) 側に存在すると判断して左側を切り捨て、逆の場合は右側を切り捨てます。この操作を収束条件まで反復することで、高精度な近似解が得られます。適用前に導関数の符号変化や離散サンプルから単峰性を確認する必要があり、数学的証明が困難な場合は実験的検証に頼るケースもあります。

Pollard's Rho による素因数分解

巨大な合成数を多項式時間で素因数分解する確率的アルゴリズムです。ミラー・ラビン素数判定法と組み合わせて、再帰的に因子を分離します。実装においては、巡回検出と GCD の頻度計算をバランスさせ、メモリ使用量を抑えつつ高速化を図ります。以下に再構成した実装例を示します。

#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <algorithm>

typedef long long int64;

int64 calc_gcd(int64 a, int64 b) {
    while (b) { int64 t = a % b; a = b; b = t; }
    return a;
}

int64 pow_mod(int64 base, int64 exp, int64 mod) {
    int64 res = 1;
    base %= mod;
    while (exp) {
        if (exp & 1) res = (__int128)res * base % mod;
        base = (__int128)base * base % mod;
        exp >>= 1;
    }
    return res;
}

bool is_prime(int64 n) {
    if (n < 2) return false;
    if (n == 2 || n == 3) return true;
    if (!(n & 1)) return false;
    int64 d = n - 1, shift = 0;
    while (!(d & 1)) { d >>= 1; shift++; }
    for (int k = 0; k < 8; ++k) {
        int64 a = 2 + rand() % (n - 2);
        int64 x = pow_mod(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool ok = false;
        for (int r = 0; r < shift - 1; ++r) {
            x = (__int128)x * x % n;
            if (x == n - 1) { ok = true; break; }
        }
        if (!ok) return false;
    }
    return true;
}

int64 rho_find_factor(int64 n) {
    if (n % 2 == 0) return 2;
    int64 x = 2, y = 2, c = 1 + rand() % (n - 1);
    int64 prod = 1, step = 1;
    while (true) {
        x = (( (__int128)x * x + c) % n);
        prod = ((__int128)prod * std::abs(x - y)) % n;
        if (++step % 128 == 0) {
            int64 g = calc_gcd(prod, n);
            if (g > 1) return g;
            y = x;
        }
        if (x == y) {
            c = 1 + rand() % (n - 1);
            y = rand() % n; x = y; prod = 1; step = 1;
        }
    }
}

void decompose(int64 val, int64& max_prime) {
    if (val <= 1) return;
    if (is_prime(val)) {
        max_prime = std::max(max_prime, val);
        return;
    }
    int64 divisor = rho_find_factor(val);
    decompose(divisor, max_prime);
    decompose(val / divisor, max_prime);
}

int main() {
    srand(time(0));
    int cases;
    scanf("%d", &cases);
    while (cases--) {
        int64 target;
        scanf("%lld", &target);
        int64 answer = 0;
        decompose(target, answer);
        if (answer == target) printf("Prime\n");
        else printf("%lld\n", answer);
    }
    return 0;
}

タグ: 同余最短経路 SOS_DP 01分数計画 三分探索 Pollard_rhoアルゴリズム

5月15日 07:19 投稿