第3回 ドワンゴからの挑戦状 本選: B - ニワンゴくんの約数

AtCoder ドワコン ドワコン3 本戦 約数 Mo's Algorithm クエリ平方分割

https://atcoder.jp/contests/dwacon2017-honsen/tasks/dwango2017final_b

問題概要

N要素の正整数の数列Xが与えられる。Q個の区間[L, R]が与えられるので、それぞれについて「その区間に含まれるすべてのXの要素の積をとったとき、その数の約数は何個か」という問いに答えよ。

N, Q <= 105

max(Xi) <= 105

解法

約数の個数は各素因数について（出現回数＋１）の積をとることで求められる。なので問題で要求されている操作そのものはシンプルで、愚直にやるならばXの要素を素因数分解したうえでそれぞれの区間において各素数の出現回数の和を取ってかけ合わせればよい。

そう考えると各素数の出現回数テーブルを持って累積和をすればよいのではと一瞬思えるが、105までの素数は104個くらいあるのでそれでは計算量的に厳しい。

ということで別の道を探してmax(Xi)の小ささに注目してみる。max(Xi)がたかだか105なので、あるひとつのXiにおいてsqrt(105)より大きい素数はたかだか1回しか現れない（もし2回以上現れたら明らかに105を越えてしまう）。また105以下の素数は104個くらいあるがsqrt(105)以下の素数は70個くらいしかない。まとめると以下のような性質の違いがsqrt(105)を境に現れるといえる。

  1. sqrt(105)以下の素因数 → 種類は少ないがひとつのXiにおいて複数種・複数回現れる可能性がある
  2. sqrt(105)より大きい素因数 → 種類は多いがひとつのXiにおいてたかだか一種類・一回までしか現れない

以上の性質の違いからこれらを場合分けして独立に解くことにするとうまくいく。

sqrt(105)以下の素因数に対して

上述のとおりこのような素因数は70個くらいしかないので単純な累積和（全部の素数の出現回数をカウントするテーブルを作って累積和）を適用できる。

sqrt(105)より大きい素因数に対して

ひとつのXiにおいてたかだか一種類・一回までしか現れないという性質から、もしいまある区間の答えがわかっているならば、その区間の端をひとつだけずらしたときの答えというのが高速に計算できる（たとえば[L, R]の右端をひとつ伸ばして[L, R+1]を見たとき、X_(R+1)の「大きい素因数」の出現回数が1増える。もしもともとの出現回数が3回だったとするとこれによって答えは「元の答え / 3 × 4」に更新できる）。

このように区間をひとつずらしたときの答えの差分更新が高速にできるとき、Mo's algorithm が使える可能性がある。このアルゴリズムの詳細については以下の記事がとてもわかりやすかった。

具体的な話はこの記事を読んでもらうとして自分の理解をざっくりまとめておくと、Mo's algorithmは「区間クエリのいい感じの並び替え方」についてのアルゴリズムである。具体的には「区間の右端か左端をひとつずつずらして差分更新していって最終的にクエリで飛んでくる全部の区間の答えを求めたい、というとき、Mo's algorithm にしたがって区間を並べ替えるとずらしの回数が(N+Q)sqrt(N)に抑えられる」というもの（Nは数列の長さ、Qはクエリ数）。

このケースはまさにMoに適した形の構造をしている。各素因数が何回出てきたかカウントする用の105のテーブルひとつと現在の答えを保持する変数ひとつを持ったうえでMoが返してくるとおりの順番で更新を行っていけばよい。

感想

うしさんありがとう

コード (C++)

Mo の部分は上の記事のものをそのまま貼っただけ（うしさんありがとう）

#include <bits/stdc++.h>
using namespace std;
#define REP(i,n) for (int i=0;i<(n);i++)
#define REP2(i,m,n) for (int i=m;i<(n);i++)
typedef long long ll;

const ll MOD = 1000000007;
const int SQ = 350;
int N, Q;
int A[101010];
int B[101010][100];
int L[101010];
int R[101010];
ll ans[101010];
vector<int> P;
ll tmp = 1;
ll cnt[101010];
ll inv[101010];

ll powmod(ll a, ll x) {
    ll ret = 1;
    while (x) {
        if (x % 2) ret = ret * a % MOD;
        a = a * a % MOD;
        x /= 2;
    }
    return ret;
}

// from https://ei1333.hateblo.jp/entry/2017/09/11/211011
struct Mo
{
    vector< int > left, right, order;
    vector< bool > v;
    int width;
    int nl, nr, ptr;

    Mo(int n) : width((int) sqrt(n)), nl(0), nr(0), ptr(0), v(n) {}

    void insert(int l, int r) /* [l, r) */
    {
        left.push_back(l);
        right.push_back(r);
    }

    /* ソート */
    void build()
    {
        order.resize(left.size());
        iota(begin(order), end(order), 0);
        sort(begin(order), end(order), [&](int a, int b)
        {
            if(left[a] / width != left[b] / width) return left[a] < left[b];
            return right[a] < right[b];
        });
    }

    /* クエリを 1 つぶんすすめて, クエリのidを返す */
    int process()
    {
        if(ptr == order.size()) return (-1);
        const auto id = order[ptr];
        while(nl > left[id]) distribute(--nl);
        while(nr < right[id]) distribute(nr++);
        while(nl < left[id]) distribute(nl++);
        while(nr > right[id]) distribute(--nr);
        return (order[ptr++]);
    }

    inline void distribute(int idx)
    {
        v[idx].flip();
        if(v[idx]) add(idx);
        else del(idx);
    }

    void add(int idx) {
        if (A[idx] == 1) return;
        tmp = tmp * inv[cnt[A[idx]] + 1] % MOD;
        cnt[A[idx]] += 1;
        tmp = tmp * (cnt[A[idx]] + 1) % MOD;
    }

    void del(int idx) {
        if (A[idx] == 1) return;
        tmp = tmp * inv[cnt[A[idx]] + 1] % MOD;
        cnt[A[idx]] -= 1;
        tmp = tmp * (cnt[A[idx]] + 1) % MOD;
    }
};

void enum_small_primes() {
    for (int i = 2; i < SQ; ++i) {
        bool ok = true;
        for (auto p: P) if (i % p == 0) {
            ok = false;
            break;
        }
        if (ok) {
            P.push_back(i);
        }
    }
}

void solve_small() {
    REP(i, N) REP(j, (int)P.size()) while (A[i] % P[j] == 0) {
        A[i] /= P[j];
        B[i+1][j] += 1;
    }
    REP(i, N) REP(j, (int)P.size()) {
        B[i+1][j] += B[i][j];
    }
    REP(i, Q) REP(j, (int)P.size()) {
        ans[i] = ans[i] * (B[R[i]+1][j] - B[L[i]][j] + 1) % MOD;
    }
}

void solve_large() {
    Mo mo = Mo(N);
    REP(i, Q) {
        mo.insert(L[i], R[i] + 1);
    }
    mo.build();
    REP(i, Q) {
        int idx = mo.process();
        ans[idx] = ans[idx] * tmp % MOD;
    }
}

void solve() {
    cin >> N >> Q;
    REP(i, N) cin >> A[i];
    REP(i, Q) cin >> L[i] >> R[i], --L[i], --R[i];
    REP(i, Q) ans[i] = 1;
    REP(i, 101010) inv[i] = powmod(i, MOD-2);

    enum_small_primes();
    solve_small();
    solve_large();

    REP(i, Q) cout << ans[i] << "\n";
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    solve(); 
}