[数论][莫比乌斯反演] 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
, fi
, pref
, prefi
;
int g
, gi
, preg
, pregi
;
int prime
, vis
, mu
;
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;
}