[数论][莫比乌斯反演] BZOJ 4816: 数字表格

Description

求

∏i=1n∏i=1mF(i,j)

1≤n,m≤106。

Solution

推一下柿子：

∏i=1n∏i=1mF(i,j)====∏d=1nF∑ni=1∑mj=1[(i,j)=d]d∏d=1nF∑⌊nd⌋k=1μ(k)⌊ndk⌋⌊mdk⌋d∏d=1nF∑nD=1μ(Dk)⌊nD⌋⌊mD⌋d∏D=1n(∏d∣DF(d)μ(Td))⌊nD⌋⌊mD⌋

设

G(D)=∏d∣DF(d)μ(Td)

所以

∏i=1n∏i=1mF(i,j)=∏D=1nG(D)⌊nD⌋⌊mD⌋

G可以O(nlnn)时间预处理出来。F也可以直接预处理出来。
这里有一个小trick：
因为每次累乘计算F−1i的时候都要乘上一个O(logP)的时间，是不能跑过去的。所以可以先算出Fi的前缀积Ti。那么就有

F−1i=T−1iTi−1,T−1i=T−1i+1Fi+1

就可以O(n)计算出来啦。

#include <bits/stdc++.h>using namespace std;const int N = 1010101;const int MOD = 1000000007;typedef long long ll;inline char get(void) {    static char buf[100000], *S = buf, *T = buf;    if (S == T) {        T = (S = buf) + fread(buf, 1, 100000, stdin);        if (S == T) return EOF;    }    return *S++;}inline void read(int &x) {    static char c; x = 0;    for (c = get(); c < '0' || c > '9'; c = get());    for (; c >= '0' && c <= '9'; c = get()) x = x * 10 + c - '0';}int f[N], fi[N], pref[N], prefi[N];int g[N], gi[N], preg[N], pregi[N];int prime[N], vis[N], mu[N];int lim = 1000000;int n, m, Pcnt, x, y, test, pos, ans;inline void Add(int &x, int t) {    x += t; while (x >= MOD) x -= MOD;}inline int Pow(int a, ll b) {    int c = 1;    while (b) {        if (b & 1) c = (ll)c * a % MOD;        b >>= 1; a = (ll)a * a % MOD;    }    return c;}inline int Inv(int x) {    return Pow(x, MOD - 2);}int main(void) {    freopen("1.in", "r", stdin);    mu[1] = 1;    for (int i = 2; i <= lim; i++) {        if (!vis[i]) {            prime[++Pcnt] = i; mu[i] = -1;        }        for (int j = 1; j <= Pcnt && (x = i * prime[j]) <= lim; j++) {            vis[x] = 1;            if (i % prime[j]) {                mu[x] = -mu[i];            } else {                mu[x] = 0; break;            }        }    }    pref[0] = prefi[0] = 1; f[0] = fi[0] = 0;    pref[1] = prefi[1] = f[1] = fi[1] = 1;    for (int i = 2; i <= lim; i++) {        Add(f[i] = f[i - 1], f[i - 2]);        pref[i] = (ll)pref[i - 1] * f[i] % MOD;    }    prefi[lim] = Inv(pref[lim]);    for (int i = lim; i >= 2; i--) {        prefi[i - 1] = (ll)prefi[i] * f[i] % MOD;        fi[i] = (ll)prefi[i] * pref[i - 1] % MOD;    }    for (int i = 1; i <= lim; i++) g[i] = gi[i] = 1;    for (int i = 1; i <= lim; i++)        for (int j = i; j <= lim; j += i) {            if (mu[j / i] == 1) {                x = f[i]; y = fi[i];            } else if (mu[j / i] == 0) {                x = y = 1;            } else {                x = fi[i]; y = f[i];            }            g[j] = (ll)g[j] * x % MOD;            gi[j] = (ll)gi[j] * y % MOD;        }    preg[0] = pregi[0] = 1;    for (int i = 1; i <= lim; i++) {        preg[i] = (ll)preg[i - 1] * g[i] % MOD;        pregi[i] = (ll)pregi[i - 1] * gi[i] % MOD;    }    read(test);    while (test--) {        read(n); read(m); ans = 1;        if (n > m) swap(n, m);        for (int i = 1; i <= n; i = pos + 1) {            pos = min(n / (n / i), m / (m / i));            ans = (ll)ans * Pow((ll)preg[pos] * pregi[i - 1] % MOD, (ll)(n / i) * (m / i)) % MOD;        }        printf("%d\n", ans);    }    return 0;}