2024-09-29

ゆるふわ競技プログラミングオンサイト at FORCIA #7 G - Dot Product Query (4D)

mofecoder 有志コン添字GCD畳み込み高速畳み込み計算約数系包除最大公約数クエリ処理問題エラトステネスの区間篩前処理ベズーの等式数学(整数問題) 見積り大事数列最適化問題最小回数・最小個数を求める高度典型個人的要復習 NoviSteps4D

コンテスト中最初に挑んだが解けず、その後も結局解けなかった。gcd convolution を思い出せなかった。

問題へのリンク

問題概要

長さ $N$ の正の整数からなる数列 $A_{1}, A_{2}, \dots, A_{N}$ 、 $B_{1}, B_{2}, \dots, B_{N}$ が与えられる。これらの数列に対して $Q$ 個のクエリが与えられる。

【クエリ】
正の整数 $C$ が与えられる。数列 $B$ の 0 個以上の要素を適当な整数に変更することで、

$\displaystyle \sum_{i=1}^{N} A_{i}B_{i} = C$

を満たすことができるか判定し、可能な場合は変更する要素の個数の最小値を求めよ。

制約

$1 \le N, Q, A_{i}, B_{i}, C \le 2 \times 10^{5}$

考えたこと

$D = |C - \sum_{i=1}^{N} A_{i}B_{i}|$ とする。まず、 $A$ 全体の最大公約数を $g$ として $g | D$ でない場合は -1 としてよい。結局は次の問題となる。

$A$ から最小個数の要素 $A_{i_{1}}, A_{i_{2}}, \dots, A_{i_{K}}$ を選んで、

$\mathrm{gcd}(A_{i_{1}}, A_{i_{2}}, \dots, A_{i_{K}}) | D$

となるようにせよ。

ここで、 $D$ のとり得るは $10^{10}$ オーダーの大きな値ではあるものの、その値の範囲は $2 \times 10^{5}$ と小さい ( $C \le 2 \times 10^{5}$ より) ことに注意しよう。

エラトステネスの区間篩のアナロジーから、 $v = 1, 2, \dots 2 \times 10^{5}$ について、

$\mathrm{gcd}(A_{i_{1}}, A_{i_{2}}, \dots, A_{i_{K}}) | v$

を満たすような最小の $K_{v}$ を求める問題へと帰着される。

また、 $K_{v}$ がそれほど大きくならないこともわかる。なぜならば、たとえば $K_{v} = 5$ という例を構成しようとしただけでも、最小で

$A_{1} = 3 \times 5 \times 7 \times 11$
$A_{2} = 2 \times 5 \times 7 \times 11$
$A_{3} = 2 \times 3 \times 7 \times 11$
$A_{4} = 2 \times 3 \times 5 \times 11$
$A_{5} = 2 \times 3 \times 5 \times 7$

と大きな値が必要になる。ちゃんと評価すると $K \le 6$ の範囲を調べれば十分だとわかる。

gcd convolution

ここまではコンテスト中に考察できていたが、結局 $K_{v} \le 6$ がわかっても、間に合う解法が思い浮かばなかった。解説で gcd convolution というワードを見て、すべてが氷解した。gcd convolution を用いると、次のことが $O(V \log \log V)$ の計算量でできる。

整数 $v = 1, 2, \dots, V$ に対する整数値関数 $f(v), g(v)$ があるとする。

$h(v) = \displaystyle \sum_{\mathrm{gcd}(i, j) = v} f(i)g(j)$

によって定義される整数値関数 $h(v)$ を求める (以降、これを $f \cdot g$ と書くことにする)。

具体的な方法は、たとえば次の記事に書いた。

qiita.com

今回は、

$f(v)$ = ( $A_{i} = v$ となる $i$ が存在するとき 1、そうでないとき 0)

とするといい感じになる。

$(f \cdot f)(v) = 1$ かつ $f(v) = 0$ であるような $v$ に対しては、 $K_{v} = 2$ と求められる
同様に、 $(f \cdot f \cdot f)(v) = 1$ かつ $(f \cdot f)(v) = 0$ であるような $v$ に対しては、 $K_{v} = 3$ と求められる
...

という感じだ。これを高々 6 回繰り返せば良い。なお、この 6 という部分はちゃんと計算量評価すると、 $A, C$ の値の最大値を $M$ として $O(\log M)$ と書ける。

以上の解法の計算量を評価しよう。区間篩パートも考慮すると $O(N + Q + M \log M)$ と評価できる。

コード

#include <bits/stdc++.h>
using namespace std;

// f(k) = sum_{GCD(i, j) = k} g(i)h(j)
// F(k) = sum_{k | i} f(i)
// F(k) = G(k)H(k)
namespace FastGCDConvolution {
    vector<bool> eratos(int N) {
        vector<bool> isprime(N, true);
        isprime[0] = isprime[1] = false;
        for (int p = 2; p < N; ++p) {
            if (!isprime[p]) continue;
            for (int i = p*2; i < N; i += p)
                isprime[i] = false;
        }
        return isprime;
    }

    template<class T> void zeta(vector<T> &v, const vector<bool> &isprime) {
        int N = (int)v.size();
        for (int p = 2; p < N; ++p) {
            if (!isprime[p]) continue;
            for (int i = (N-1)/p; i >= 1; --i)
                v[i] += v[i*p];
        }
    }

    template<class T> void mebius(vector<T> &v, const vector<bool> &isprime) {
        int N = (int)v.size();
        for (int p = 2; p < N; ++p) {
            if (!isprime[p]) continue;
            for (int i = 1; i*p < N; ++i)
                v[i] -= v[i*p];
        }
    }

    template<class T> vector<T> mul(const vector<T> &a, const vector<T> &b) {
        int N = max((int)a.size(), (int)b.size());
        const auto &isprime = eratos(N);
        vector<T> A(N, 0), B(N, 0);
        for (int i = 0; i < a.size(); ++i) A[i] = a[i];
        for (int i = 0; i < b.size(); ++i) B[i] = b[i];
        zeta(A, isprime), zeta(B, isprime);
        vector<T> C(N);
        for (int i = 1; i < N; ++i) C[i] = A[i] * B[i];
        mebius(C, isprime);
        return C;
    }
};

const int MAX = 210000;
int main() {
    long long N, Q, C, base = 0, g = 0;
    cin >> N;
    vector<long long> A(N), B(N), va(MAX, 0), tmp(MAX, 0), num(MAX, -1);
    for (int i = 0; i < N; i++) {
        cin >> A[i]; 
        g = gcd(g, A[i]);
        va[A[i]] = tmp[A[i]] = num[A[i]] = 1;
    }
    for (int i = 0; i < N; i++) cin >> B[i], base += A[i] * B[i];

    for (int iter = 2; iter < 10; iter++) {
        tmp = FastGCDConvolution::mul(tmp, va);
        for (int v = 0; v < MAX; v++) {
            if (tmp[v] > 0 && num[v] == -1) num[v] = iter;
        }
    }

    if (base <= MAX)  {
        vector<long long> res(MAX, -1);
        res[0] = 0;
        for (int v = 1; v < MAX; v++) {
            if (num[v] < 1) continue;
            for (int v2 = v; v2 < MAX; v2 += v) {
                if (res[v2] == -1 || res[v2] > num[v]) res[v2] = num[v];
            }
        }
        cin >> Q;
        while (Q--) {
            cin >> C;
            cout << res[abs(base - C)] << endl;
        }
    } else {
        // res[i] := base - MAX + i の答え
        vector<long long> res(MAX, -1);
        for (long long v = 1; v < MAX; ++v) {
            if (num[v] < 1) continue;
            for (long long v2 = (base - MAX + v - 1) / v * v; v2 < base; v2 += v) {
                if (res[v2 - base + MAX] == -1 || res[v2 - base + MAX] > num[v]) {
                    res[v2 - base + MAX] = num[v];
                }
            }
        }
        cin >> Q;
        while (Q--) {
            cin >> C;
            cout << res[MAX - C] << endl;
        }
    }
}

けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

ゆるふわ競技プログラミングオンサイト at FORCIA #7 G - Dot Product Query (4D)

問題概要

制約

考えたこと

gcd convolution

コード