問題概要

整数 $M$ が与えられる。空の vector があって

$1$ 以上 $M$ 以下の整数の中から一様ランダムに 1 つ選び、
それを vector に push する
このとき vector 内の要素の最大公約数が 1 ならば処理を終える

というのを行う。処理を終えたときの vector のサイズの期待値を求めよ ( $p/q$ の形で表して、 $pq^{-1} ({\rm mod}. 1000000007)$ を求めよ)。

制約

$1 \le M \le 10^{5}$

解法

この記事に書いた解法を実践する

drken1215.hatenablog.com

今回の問題を言い換えると

毎回の試行で $1$ 以上 $M$ 以下の整数を挙げていく。このとき $M$ 以下の素数 $p_1, p_2, \dots, p_n$ について

ミッション 1: $p_1$ の倍数でないものが登場したらクリア
ミッション 2: $p_2$ の倍数でないものが登場したらクリア
...
ミッション n: $p_n$ の倍数でないものが登場したらクリア

というゲームのすべてのミッションをクリアするまでに必要な試行回数の期待値を求めよ

ということになる。このままでは $n$ が膨大なので一見 $2^{n}$ 通りの包除原理をやるのがつらそうなのだが、よく考えると $p$ の積が $M$ を超えるものについては考える必要が無い。

$1$ 以上 $M$ 以下の整数のうち、メビウス関数値が 0 でないものについて包除原理を行えば良い。

#include <iostream>
#include <vector>
using namespace std;


// エラトステネスやりながら「メビウス関数」を求める
const int MAX = 101010;
bool is_prime[MAX];
int min_factor[MAX];
int meb[MAX];
vector<int> preprocess(int n = MAX) {
    vector<int> res;
    for (int i = 0; i < n; ++i) is_prime[i] = true, min_factor[i] = -1,  meb[i] = 1;
    is_prime[0] = false; is_prime[1] = false;
    min_factor[0] = 0, min_factor[1] = 1;
    meb[0] = 0; meb[1] = 1;
    for (int i = 2; i < n; ++i) {
        if (!is_prime[i]) continue;
        res.push_back(i);
        min_factor[i] = i;
        meb[i] = -1;
        for (int j = i*2; j < n; j += i) {
            is_prime[j] = false;
            meb[j] *= -1;
            if (min_factor[j] == -1) min_factor[j] = i;
            if ((j/i) % i == 0) meb[j] = 0;
        }
    }
    return res;
}

long long modinv(long long a, long long mod) {
    long long b = mod, u = 1, v = 0;
    while (b) {
        long long t = a/b;
        a -= t*b; swap(a, b);
        u -= t*v; swap(u, v);
    }
    u %= mod;
    if (u < 0) u += mod;
    return u;
}

int main() {
    const long long MOD = 1000000007;
    long long M; cin >> M;
    preprocess();

    long long res = 1;
    for (long long v = 2; v <= M; ++v) {
        if (meb[v] == 0) continue;
        long long r = (M/v) * modinv(M, MOD) % MOD;
        long long add = r * modinv((MOD + 1 - r) % MOD, MOD) % MOD;
        if (meb[v] == -1) res = (res + add) % MOD;
        else res = (res - add + MOD) % MOD;
    }
    cout << res << endl;
}

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

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

Codeforces #548 Div. 2 D - Steps to One (R2300)

問題概要

制約

解法