Description

Tweetuzki 有一個袋子，袋子中有 \(N\) 個無差別的球。Tweetuzki 每次隨機取出一個球后放回。求取遍所有球的期望次數。

取遍是指，袋子中所有球都被取出來過至少一次。

Input

一行一個整數 \(N\)

Output

一行一個整數，表示期望次數對 \(200403013\) 取模的結果

Hint

\(1~\leq~N~\leq~10^7\)

Solution

期望嘛，就是你xjb亂推好幾個看起來很靠譜事實上最多有一個是對的甚至可能一個也不對的式子，然後用蒙特卡洛方法拍一拍，挑一個最靠譜的交上去

感謝@_rqy提供指導

先給出結論：

\[E(n)~=~n~\times~(\sum_{i=1}^{n}~\frac{1}{i})\]

證明如下：

引理：對於一次試驗，有 \(p~(0~\leq~p~\leq~1)\) 的概率成功，有 \(1-p\) 的概率不成功，如果一次試驗後不成功則重複進行試驗，否則停止。則期望在 \(\frac{1}{p}\) 步後停止試驗

證明：

設期望 \(x\) 步後停止試驗，根據已知以及期望的線性相加性，可得方程：

\[x~=~1~+~p~\times~0~+~(1~-~p)~\times~x\]

上述方程的直觀解釋是當前進行完一步後有 \(p\) 的概率再進行 \(0\) 次（停止），\(1~-~p\) 的概率繼續進行，由於下次試驗和這次試驗完全相同，所以有 \(1~-~p\) 的概率期望再進行 \(x\)

解方程，得

\[x~=~\frac{1}{p}\]

引理證畢。

於是考慮摸球，首先有 \(n\) 個球沒有被摸到的時候，摸一次摸到新球的概率為 \(1\)，根據引理，期望摸 \(1\) 次。

然後剩下 \(n~-~1\) 個球沒有被摸到，摸一次摸到新球的概率為 \(\frac{n-1}{n}\)，根據引理，期望摸 \(\frac{n}{n-1}\) 次。

然後剩下 \(n~-~2\) 個球沒有摸到，摸一次摸到新球的概率為 \(\frac{n-2}{n}\) ，根據引理，期望摸 \(\frac{n}{n-2}\) 次。

……（一位讀者及時砸壞了復讀機，停止了復讀）

按照球的個數進行數學歸納，根據期望的可加性（非線性相加性），可以得到摸遍 \(n\)

\[E(n)~=~n~+~\frac{n}{n-1}~+~\frac{n}{n-1}~+~\dots~+~\frac{n}{1}\]

提取公因數 \(n\)，可得

\[E(n)~=~n~\times~(\sum_{i=1}^{n}~\frac{1}{i})\]

證畢。

然鵝最後一個式子只是為了讓結論變得優美一點，事實上用倒數第二個式子就可以做題了233

於是先篩一遍逆元，然後求答案即可。

Code

#include <cstdio>
#ifdef ONLINE_JUDGE
#define freopen(a, b, c)
#endif
#define rg register
#define ci const int
#define cl const long long

typedef long long int ll;

namespace IPT {
    const int L = 1000000;
    char buf[L], *front=buf, *end=buf;
    char GetChar() {
        if (front == end) {
            end = buf + fread(front = buf, 1, L, stdin);
            if (front == end) return -1;
        }
        return *(front++);
    }
}

template <typename T>
inline void qr(T &x) {
    rg char ch = IPT::GetChar(), lst = ' ';
    while ((ch > '9') || (ch < '0')) lst = ch, ch=IPT::GetChar();
    while ((ch >= '0') && (ch <= '9')) x = (x << 1) + (x << 3) + (ch ^ 48), ch = IPT::GetChar();
    if (lst == '-') x = -x;
}

template <typename T>
inline void ReadDb(T &x) {
    rg char ch = IPT::GetChar(), lst = ' ';
    while ((ch > '9') || (ch < '0')) lst = ch, ch = IPT::GetChar();
    while ((ch >= '0') && (ch <= '9')) x = x * 10 + (ch ^ 48), ch = IPT::GetChar();
    if (ch == '.') {
        ch = IPT::GetChar();
        double base = 1;
        while ((ch >= '0') && (ch <= '9')) x += (ch ^ 48) * ((base *= 0.1)), ch = IPT::GetChar();
    }
    if (lst == '-') x = -x;
}

namespace OPT {
    char buf[120];
}

template <typename T>
inline void qw(T x, const char aft, const bool pt) {
    if (x < 0) {x = -x, putchar('-');}
    rg int top=0;
    do {OPT::buf[++top] = x % 10 + '0';} while (x /= 10);
    while (top) putchar(OPT::buf[top--]);
    if (pt) putchar(aft);
}

const int maxn = 10000010;
const int MOD = 20040313;

int n, ans;
int inv[maxn];

void GetInv();

int main() {
    freopen("1.in", "r", stdin);
    qr(n);
    GetInv();
    for (rg int i = 0; i < n; ++i) {
        ans = (ans + 1ll * n * inv[n - i]) % MOD;
    }
    qw(ans, '\n', true);
}

void GetInv() {
    inv[1] = 1;
    for(int i = 2; i <= n; ++i) inv[i] = 1ll * (MOD - MOD / i) * inv[MOD % i] % MOD;
}

Summary

對於一次試驗，有 \(p~(0~\leq~p~\leq~1)\) 的概率成功，有 \(1-p\) 的概率不成功，如果一次試驗後不成功則重複進行試驗，否則停止。則期望在 \(\frac{1}{p}\) 步後停止試驗