莫比乌斯反演2
「LuoGu P5221」Product
给出一个 \(n\) 。求
\[\prod_{i=1}^n\prod_{j=1}^n\dfrac{\operatorname{lcm}(i,j)}{\gcd(i,j)}
\]
先化简。
\[\prod_{i=1}^n\prod_{j=1}^n\dfrac{\frac{ij}{\gcd(i,j)}}{\gcd(i,j)}
\]
\[\prod_{i=1}^n\prod_{j=1}^n\dfrac{ij}{\gcd(i,j)^2}
\]
\[\dfrac{\prod_{i=1}^n\prod_{j=1}^nij}{(\prod_{i=1}^n\prod_{j=1}^n\gcd(i,j))^2}
\]
\[\dfrac{\prod_{i=1}^n(i^n\times n!)}{(\prod_{d=1}^nd^{\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^\frac{n}{d}[\gcd(i,j)=1]})^2}
\]
\[\dfrac{(n!)^{2n}}{(\prod_{d=1}^nd^{\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^\frac{n}{d}[\gcd(i,j)=1]})^2}
\]
辣个指数就是仪仗队啊。
我真的是个傻逼！
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL long long
#define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1 << 19, stdin), p1 == p2) ? EOF : *p1++)
char buf[1 << 19], *p1 = buf, *p2 = buf;
LL read() {
LL ret = 0, f = 1;
char ch = getchar();
while (!isdigit(ch)) {
if (ch == '-')
f = -1;
ch = getchar();
}
for (; isdigit(ch); ch = getchar()) ret = ret * 10 + ch - 48;
return ret * f;
}
LL qkpow(LL x, LL power, LL mod) {
x %= mod;
LL ans = 1;
for (; power; power >>= 1, (x *= x) %= mod)
if (power & 1)
(ans *= x) %= mod;
return ans;
}
const int MAXN = 1e6 + 5;
const int mod = 104857601;
bool vis[MAXN];
int n, len, phi[MAXN], prime[MAXN];
void seive1() {
phi[1] = 1;
for (int i = 2; i <= MAXN - 5; i++) {
if (!vis[i])
prime[++len] = i, phi[i] = i - 1;
for (int j = 1; j <= len && i * prime[j] <= MAXN - 5; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
for (int i = 2; i <= MAXN - 5; i++) phi[i] += phi[i - 1], phi[i] %= (mod - 1);
}
int main() {
seive1();
scanf("%d", &n);
int ans = 1, fact = 1;
for (int i = 2; i <= n; i++) fact = (LL)fact * (LL)i % mod;
ans = (LL)ans * (LL)qkpow(fact, n, mod) % mod;
for (int i = 1; i <= n; i++) ans = (LL)ans * (LL)qkpow(i, n, mod) % mod;
int total = 1;
for (int d = 1; d <= n; d++)
total = (LL)total * (LL)qkpow(d, (2 * phi[n / d] - 1) % (mod - 1), mod) % mod;
total = (LL)total * (LL)total % mod;
printf("%d", (LL)ans * (LL)qkpow(total, mod - 2, mod) % mod);
return 0;
}
还可以优化啊，整除分块就可以了。
这里附上优化后的代码。
\(PS\) :多次询问请提交 「LuoGu P4917」天守阁的地板 。
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL long long
#define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1 << 19, stdin), p1 == p2) ? EOF : *p1++)
char buf[1 << 19], *p1 = buf, *p2 = buf;
LL read() {
LL ret = 0, f = 1;
char ch = getchar();
while (!isdigit(ch)) {
if (ch == '-')
f = -1;
ch = getchar();
}
for (; isdigit(ch); ch = getchar()) ret = ret * 10 + ch - 48;
return ret * f;
}
LL qkpow(LL x, LL power, LL mod) {
x %= mod;
LL ans = 1;
for (; power; power >>= 1, (x *= x) %= mod)
if (power & 1)
(ans *= x) %= mod;
return ans;
}
const int MAXN = 1e6 + 5;
const int mod = 19260817;
bool vis[MAXN];
int n, len, phi[MAXN], prime[MAXN], ans[MAXN], fact[MAXN];
void seive1() {
phi[1] = 1;
for (int i = 2; i <= MAXN - 5; i++) {
if (!vis[i])
prime[++len] = i, phi[i] = i - 1;
for (int j = 1; j <= len && i * prime[j] <= MAXN - 5; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
for (int i = 2; i <= MAXN - 5; i++) phi[i] += phi[i - 1], phi[i] %= (mod - 1);
}
int calc(int n) {
int ans = 1;
for (int i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
ans = ((LL)ans * (LL)qkpow((LL)fact[j] % mod * (LL)qkpow(fact[i - 1], mod - 2, mod) % mod,
(2 * phi[n / i] - 1) % (mod - 1), mod)) %
mod;
}
return ans;
}
int main() {
int T;
seive1();
scanf("%d", &T);
fact[0] = fact[1] = 1, ans[1] = 1;
for (int i = 2; i <= MAXN - 5; i++)
fact[i] = (LL)fact[i - 1] * (LL)i % mod, ans[i] = qkpow(fact[i], 2 * i, mod);
while (T--) {
scanf("%d", &n);
int total = calc(n);
total = (LL)total * (LL)total % mod;
printf("%d\n", (LL)ans[n] * (LL)qkpow(total, mod - 2, mod) % mod);
}
return 0;
}
「LuoGu P4240」毒瘤之神的考验
\(T\) 次询问，每次询问给出 \(n,m\) 。
求
\[\sum_{i=1}^n\sum_{j=1}^m\varphi(ij)
\]
这个 \(\varphi(ij)\) 不是很好做，考虑拆开。
这里有两种处理方式：
1.直接拆成\(\varphi(i)\varphi(j)\) ，然后利用性质 \(\varphi(n)=n\prod_{d|n}[d\in\operatorname{Prime}](1-\dfrac{1}{d})\) 补上 \(\varphi(i)\) , \(\varphi(j)\) 都有的质因子。
即
\[\varphi(ij)=\varphi(i)\varphi(j)\prod_{k=1}^{\min(i,j)}[k\in\operatorname{Prime}][k|i][k|j]\dfrac{k}{k-1}
\]
2.有一个公式：
\[\varphi(ij)=\dfrac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(i,j))}
\]
我用的是第一种，代回原式可得：
\[\sum_{i=1}^n\sum_{j=1}^m\varphi(i)\varphi(j)\prod_{k=1}^{\min(i,j)}[k\in\operatorname{Prime}][k|i][k|j]\dfrac{k}{k-1}
\]
\[\sum_{i=1}^n\sum_{j=1}^m\varphi(i)\varphi(j)\prod_{k=1}^{\min(i,j)}[k\in\operatorname{Prime}][k|\gcd(i,j)]\dfrac{k}{k-1}
\]
\[\sum_{d=1}^{\min(n,m)}\prod_{k|d}[k\in\operatorname{Prime}]\dfrac{k}{k-1}\sum_{i=1}^n\sum_{j=1}^m\varphi(i)\varphi(j)[\gcd(i,j)=d]
\]
\[\sum_{d=1}^{\min(n,m)}\prod_{k|d}[k\in\operatorname{Prime}]\dfrac{k}{k-1}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\varphi(id)\varphi(jd)[\gcd(i,j)=1]
\]
直接莫比乌斯反演。
\[\sum_{d=1}^{\min(n,m)}\prod_{k|d}[k\in\operatorname{Prime}]\dfrac{k}{k-1}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\varphi(id)\varphi(jd)\sum_{D|\gcd(i,j)}\mu(D)
\]
\[\sum_{d=1}^{\min(n,m)}\prod_{k|d}[k\in\operatorname{Prime}]\dfrac{k}{k-1}\sum_{D=1}^{\min(\frac{n}{d},\frac{m}{d})}\mu(D)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\varphi(id)\varphi(jd)[D|i][D|j]
\]
\[\sum_{d=1}^{\min(n,m)}\prod_{k|d}[k\in\operatorname{Prime}]\dfrac{k}{k-1}\sum_{D=1}^{\min(\frac{n}{d},\frac{m}{d})}\mu(D)\sum_{i=1}^{\frac{n}{dD}}\sum_{j=1}^{\frac{m}{dD}}\varphi(idD)\varphi(jdD)
\]
\[\sum_{T=1}^{\min(n,m)}\sum_{d|T}((\prod_{k|d}[k\in\operatorname{Prime}]\dfrac{k}{k-1})\mu(\dfrac{T}{d}))\sum_{i=1}^{\frac{n}{T}}\varphi(iT)\sum_{j=1}^{\frac{m}{T}}\varphi(jT)
\]
\(\sum_{d|T}((\prod_{k|d}[k\in\operatorname{Prime}]\dfrac{k}{k-1})\mu(\dfrac{T}{d}))\) 这一坨直接 \(\mathcal{O(n\log n)}\) 的筛法预处理。
然后对 \(T\) 进行数论分块，右边那坨 \(\varphi\) 求和直接暴力算，至此我们就得到了单次询问 \(\mathcal{O(n+m)}\) 的算法。
怎么优化呢？显然瓶颈在后面那一坨 \(\varphi\) 求和里。这时我们就需要乱搞了。
设一个阈值 \(B\) ，函数 \(F(x),G(x,y),S(x,y,z)\) 。
\[F(x)=\sum_{d|x}((\prod_{k|d}[k\in\operatorname{Prime}]\dfrac{k}{k-1})\mu(\dfrac{x}{d}))
\]
\[G(x,y)=\sum_{i=1}^x\varphi(iy)
\]
\[S(x,y,z)=\sum_{T=1}^{z}\sum_{d|T}((\prod_{k|d}[k\in\operatorname{Prime}]\dfrac{k}{k-1})\mu(\dfrac{T}{d}))\sum_{i=1}^{y}\varphi(iT)\sum_{j=1}^{x}\varphi(jT)
\]
显然有
\[G(x,y)=G(x,y-1)+\varphi(xy)
\]
\[S(x,y,z)=S(x,y,z-1)+F(z)G(z,x)G(z,y)
\]
然后对于 \(\dfrac{m}{B}\) 的部分直接用 \(G\) 函数暴力算，后面的部分整除分块即可。
时间复杂度 \(\mathcal{O(n+n\log n+nB^2+T(\sqrt{n}+\dfrac{n}{B}))}\)
#include <cstdio>
#include <vector>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL long long
const int MAXN = 1e5 + 5;
bool vis[MAXN];
int len, mu[MAXN], prime[MAXN], B = 100;
const LL mod = 998244353;
LL qkpow(LL x, LL power, LL mod) {
x %= mod;
LL ans = 1;
for (; power; power >>= 1, (x *= x) %= mod)
if (power & 1)
(ans *= x) %= mod;
return ans;
}
LL Prime[MAXN], f[MAXN], phi[MAXN];
vector<int> G[MAXN], S[105][105];
void Getmu() {
mu[1] = 1, phi[1] = 1;
Prime[1] = 1;
for (int i = 2; i <= MAXN - 5; i++) {
Prime[i] = 1;
if (!vis[i])
prime[++len] = i, mu[i] = -1, phi[i] = i - 1;
for (int j = 1; j <= len && i * prime[j] <= MAXN - 5; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
mu[i * prime[j]] = -mu[i];
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
for (int i = 1; i <= len; i++)
for (int j = prime[i]; j <= MAXN - 5; j += prime[i])
Prime[j] = (Prime[j] * prime[i] % mod * qkpow(prime[i] - 1, mod - 2, mod) % mod) % mod;
for (int i = 1; i <= MAXN - 5; i++)
for (int j = i; j <= MAXN - 5; j += i)
f[j] = (f[j] + Prime[i] * mu[j / i]) % mod, f[j] = (f[j] + mod) % mod;
for (int i = 1; i <= MAXN - 5; i++) {
G[i].push_back(0);
for (int j = 1; j <= (MAXN - 5) / i; j++) G[i].push_back((G[i][j - 1] + phi[i * j]) % mod);
}
for (int i = 1; i <= B; i++) {
for (int j = 1; j <= B; j++) {
S[i][j].push_back(0);
for (int k = 1; k <= (MAXN - 5) / max(i, j); k++) {
S[i][j].push_back((S[i][j][k - 1] + f[k] * G[k][i] % mod * G[k][j]) % mod);
}
}
}
}
int T;
LL n, m;
LL calc(LL n, LL m) {
LL ans = 0;
for (int i = 1; i <= max(n, m) / B; i++) ans = (ans + f[i] * G[i][n / i] % mod * G[i][m / i] % mod) % mod;
for (int i = max(n, m) / B + 1, j; i <= min(n, m); i = j + 1) {
j = min(n / (n / i), m / (m / i));
ans = (ans + (S[n / i][m / i][j] - S[n / i][m / i][i - 1]) % mod) % mod;
}
return (ans + mod) % mod;
}
int main() {
Getmu();
scanf("%d", &T);
while (T--) {
scanf("%lld %lld", &n, &m);
printf("%lld\n", calc(n, m));
}
return 0;
}
「LuoGu P4449」于神之怒加强版
\(T\) 次询问与一个给定的 \(k\) ，每次询问给出 \(n,m\) 求：
\[\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k
\]
\[\sum_{d=1}^{\min(n,m)}d^k\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]
\]
\[\sum_{d=1}^{\min(n,m)}d^k\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[\gcd(i,j)=1]
\]
\[\sum_{d=1}^{\min(n,m)}d^k\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\sum_{D|\gcd(i,j)}\mu(D)
\]
\[\sum_{d=1}^{\min(n,m)}d^k\sum_{D=1}^{\min(\frac{n}{d},\frac{m}{d})}\mu(D)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[D|i][D|j]
\]
\[\sum_{d=1}^{\min(n,m)}d^k\sum_{D=1}^{\min(\frac{n}{d},\frac{m}{d})}\mu(D)\lfloor\dfrac{n}{dD}\rfloor\lfloor\dfrac{m}{dD}\rfloor
\]
\[\sum_{T=1}^{\min(n,m)}\sum_{d|T}d^k\mu(\dfrac{T}{d})\lfloor\dfrac{n}{T}\rfloor\lfloor\dfrac{m}{T}\rfloor
\]
因为 \(n,m\) 的最大值高达 \(5\times 10^6\) ， \(\mathcal{O(n\log n)}\) 的筛法会被卡的死死的，于是考虑线性筛。
设 \(f(T)=\sum_{d|T}d^k\mu(\dfrac{T}{d})\) 。
发现 \(f\) 函数是一个积性函数，可以线性筛，于是我们需要推出当 \(p\) 是素数时 \(f(p),f(p^q)(q>1)\) 的值。
\[f(p)=1^k\mu(p)+p^k\mu(1)=p^k-1
\]
\[f(p^q)=\sum_{i=0}^q(p^i)^k\mu(p^{q-i})=\sum_{i=0}^qp^{ik}\mu(p^{q-i})=p^{qk}-p^{qk-k}
\]
然后直接线性筛就好了，最后每次询问对 \(T\) 整除分块。
#include <cstdio>
#include <vector>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL long long
const int MAXN = 5e6 + 5;
bool vis[MAXN];
int len, mu[MAXN], prime[MAXN], low[MAXN], cnt[MAXN];
const LL mod = 1e9 + 7;
LL qkpow(LL x, LL power, LL mod) {
x %= mod;
LL ans = 1;
for (; power; power >>= 1, (x *= x) %= mod)
if (power & 1)
(ans *= x) %= mod;
return ans;
}
LL k, phi[MAXN], f[MAXN];
void Getmu() {
mu[1] = 1, f[1] = 1;
for (int i = 2; i <= MAXN - 5; i++) {
if (!vis[i])
prime[++len] = i, mu[i] = -1, f[i] = qkpow(i, k, mod) - 1, cnt[i] = 1, low[i] = i;
for (int j = 1; j <= len && i * prime[j] <= MAXN - 5; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
low[i * prime[j]] = low[i] * prime[j];
mu[i * prime[j]] = 0;
if (low[i] == i) {
cnt[i * prime[j]] = cnt[i] + 1;
f[i * prime[j]] = (qkpow(prime[j], cnt[i * prime[j]] * k, mod) -
qkpow(prime[j], cnt[i * prime[j]] * k - k, mod)) %
mod;
} else
f[i * prime[j]] = f[i / low[i]] * f[low[i] * prime[j]] % mod;
break;
}
mu[i * prime[j]] = -mu[i];
low[i * prime[j]] = prime[j];
f[i * prime[j]] = f[i] * f[prime[j]] % mod;
}
}
for (int i = 2; i <= MAXN - 5; i++) f[i] = (f[i] + f[i - 1]) % mod;
}
int T;
LL n, m;
LL calc(LL n, LL m) {
LL ans = 0;
for (LL i = 1, j; i <= min(n, m); i = j + 1) {
j = min(n / (n / i), m / (m / i));
ans = (ans + (f[j] - f[i - 1]) % mod * (n / i) % mod * (m / i) % mod) % mod;
}
return ans;
}
signed main() {
scanf("%lld %lld", &T, &k);
Getmu();
while (T--) {
scanf("%lld %lld", &n, &m);
printf("%lld\n", (calc(n, m) % mod + mod) % mod);
}
return 0;
}
「SDOI 2014」数表
话说为啥 \(\operatorname{SD}\) 这么爱出反演啊？
题意
有一张 \(n\times m\) 的数表，其第 \(i\) 行第 \(j\) 列 \((1 \le i \le n\) ， \(1\le j\le m)\) 的数值为能同时整除 \(i\) 和 \(j\) 的所有自然数之和。给定 \(a\) ，计算数表中不大于 \(a\) 的数之和。
\(1\le n,m\le 10^5\) ， \(1\le Q\le 2\times 10^4\) 。
神仙题。
\[\sum_{i=1}^n\sum_{j=1}^m[\sum_{x}[x|i][x|j]x\le a]\sum_{x}[x|i][x|j]x
\]
\[\sum_{i=1}^n\sum_{j=1}^m[\sum_{x}[x|\gcd(i,j)]x\le a]\sum_{x}[x|\gcd(i,j)]x
\]
\[\sum_{d=1}^{\min(n,m)}[\sum_{x}[x|d]x\le a]\sum_{x}[x|d]x\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]
\]
\[\sum_{d=1}^{\min(n,m)}[\sigma(d)\le a]\sigma(d)\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]
\]
\[\sum_{d=1}^{\min(n,m)}[\sigma(d)\le a]\sigma(d)\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]
\]
\[\sum_{d=1}^{\min(n,m)}[\sigma(d)\le a]\sigma(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[\gcd(i,j)=1]
\]
上反演。
\[\sum_{d=1}^{\min(n,m)}[\sigma(d)\le a]\sigma(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\sum_{D|\gcd(i,j)}\mu(D)
\]
\[\sum_{d=1}^{\min(n,m)}[\sigma(d)\le a]\sigma(d)\sum_{D=1}^{\min(\frac{n}{d},\frac{m}{d})}\mu(D)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[D|i][D|j]
\]
\[\sum_{d=1}^{\min(n,m)}[\sigma(d)\le a]\sigma(d)\sum_{D=1}^{\min(\frac{n}{d},\frac{m}{d})}\mu(D)\lfloor\dfrac{n}{dD}\rfloor\lfloor\dfrac{m}{dD}\rfloor
\]
\[\sum_{T=1}^{\min(n,m)}(\sum_{d|T}[\sigma(d)\le a]\sigma(d)\mu(\dfrac{T}{d}))\lfloor\dfrac{n}{T}\rfloor\lfloor\dfrac{m}{T}\rfloor
\]
重点是中间那一坨很烦， \(a\) 不是固定的。
考虑把询问离线下来，按照 \(a\) 从小到大排序，在把 \(1\) 到 \(100000\) 每个数的约数和都预处理出来并排序，看能不能加进去，能加就加，最后对于每个询问整除分块求解，因为要用 \((\sum_{d|T}[\sigma(d)\le a]\sigma(d)\mu(\dfrac{T}{d}))\) 的前缀和，所以用树状数组维护。
请注意你的常数因子大小！
#include <cstdio>
#include <vector>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL unsigned int
const int MAXN = 1e5 + 5;
bool vis[MAXN];
int len, mu[MAXN], prime[MAXN], low[MAXN], cnt[MAXN];
LL mod = 2147483648ll;
#define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1 << 19, stdin), p1 == p2) ? EOF : *p1++)
char buf[1 << 19], *p1 = buf, *p2 = buf;
LL read() {
LL ret = 0, f = 1;
char ch = getchar();
while (!isdigit(ch)) {
if (ch == '-')
f = -1;
ch = getchar();
}
for (; isdigit(ch); ch = getchar()) ret = ret * 10 + ch - 48;
return ret * f;
}
struct BIT {
LL tree[MAXN];
void MemSet() { fill(tree + 1, tree + 100000 + 1, 0); }
void Add(int i, LL x) {
while (i <= 100000) {
tree[i] += x, i += i & -i;
tree[i] %= mod;
}
}
LL query(int i) {
LL s = 0;
while (i > 0) {
s += tree[i];
s %= mod;
i -= i & -i;
}
return s % mod;
}
} bit;
struct node {
LL n, m, a, Id;
} problem[20005];
struct node1 {
LL sigma, Id;
} fuck[MAXN];
bool cmp(node a, node b) { return a.a < b.a; }
bool Cmp(node1 a, node1 b) { return a.sigma < b.sigma; }
void Getmu() {
mu[1] = 1;
for (int i = 2; i <= MAXN - 5; i++) {
if (!vis[i])
prime[++len] = i, mu[i] = -1;
for (int j = 1; j <= len && i * prime[j] <= MAXN - 5; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
}
mu[i * prime[j]] = -mu[i];
}
}
for (int i = 1; i <= MAXN - 5; i++) fuck[i].Id = i;
for (int i = 1; i <= MAXN - 5; i++)
for (int j = i; j <= MAXN - 5; j += i) fuck[j].sigma += i;
}
int T;
LL n, m;
LL calc(LL n, LL m) {
LL ans = 0;
for (LL i = 1, j; i <= min(n, m); i = j + 1) {
LL mmp = n / i, mmp1 = m / i;
j = min(n / mmp, m / mmp1);
ans = (ans + ((bit.query(j) - bit.query(i - 1)) % mod) * (mmp * mmp1 % mod)) % mod;
}
return ans;
}
LL ans[MAXN];
void Add(int x, LL y) {
for (int i = 1; i * x <= 100000; i++) bit.Add(i * x, y * mu[i] % mod);
}
signed main() {
Getmu();
T = read();
for (int i = 1; i <= T; i++)
problem[i].Id = i, problem[i].n = read(), problem[i].m = read(), problem[i].a = read();
sort(problem + 1, problem + T + 1, cmp);
int l = 1;
sort(fuck + 1, fuck + 100001, Cmp);
for (int i = 1; i <= T; i++) {
while (fuck[l].sigma <= problem[i].a && l <= 100000) Add(fuck[l].Id, fuck[l].sigma % mod), l++;
ans[problem[i].Id] = calc(problem[i].n, problem[i].m);
}
for (int i = 1; i <= T; i++) printf("%lld\n", (ans[i] % mod + mod) % mod);
return 0;
}
「LuoGu P3768」简单的数学题
给出 \(p,n\) 。
求
\[(\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j))\bmod p
\]
模数先不管，放一边。发现这题跟 \(Crash\) 的数字表格 很像，只是变成了乘以 \(\gcd(i,j)\) 。
\[(\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j))\bmod p
\]
\[\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^nij[\gcd(i,j)=d]
\]
\[\sum_{d=1}^nd\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ijd^2[\gcd(i,j)=1]
\]
\[\sum_{d=1}^nd^3\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ij[\gcd(i,j)=1]
\]
\[\sum_{d=1}^nd^3\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ij\sum_{D|\gcd(i,j)}\mu(D)
\]
\[\sum_{d=1}^nd^3\sum_{D=1}^{\frac{n}{d}}\mu(D)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ij[D|i][D|j]
\]
\[\sum_{d=1}^nd^3\sum_{D=1}^{\frac{n}{d}}\mu(D)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ij[D|i][D|j]
\]
后面这坨跟 \(Crash\) 的数字表格 的处理方式一模一样。
\[\sum_{d=1}^nd^3\sum_{D=1}^{\frac{n}{d}}\mu(D)D^2(\dfrac{(1+\lfloor\frac{n}{dD}\rfloor)\lfloor\frac{n}{dD}\rfloor}{2})^2
\]
\[\sum_{T=1}^n\sum_{d|T}d\mu(\dfrac{T}{d})T^2(\dfrac{(1+\lfloor\frac{n}{T}\rfloor)\lfloor\frac{n}{T}\rfloor}{2})^2
\]
\[\sum_{T=1}^n\varphi(T)T^2(\dfrac{(1+\lfloor\frac{n}{T}\rfloor)\lfloor\frac{n}{T}\rfloor}{2})^2
\]
因为 \(n\le 10^{10}\) ，所以要杜教筛。
构造
\[f=\varphi\cdot\operatorname{id}\cdot\operatorname{id}
\]
\[g=\operatorname{id}\cdot\operatorname{id}
\]
\[(f*g)(n)=\sum_{d|n}(\varphi(d)\cdot d^2)\dfrac{n^2}{d^2}=\sum_{d|n}\varphi(d)n^2=n^3
\]
然后一个问题就是求 \(n^3\) 的前缀和。
\[\sum_{i=1}^ni^3=(\dfrac{n(n+1)}{2})^2
\]
套上整除分块即可。
真 \(TM\) 恶心！
#include <cstdio>
#include <map>
#include <iostream>
#include <algorithm>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>
using namespace std;
using namespace __gnu_pbds;
#define LL long long
LL qkpow(LL x, LL power, LL mod) {
x %= mod;
LL ans = 1;
for (; power; power >>= 1, (x *= x) %= mod)
if (power & 1)
(ans *= x) %= mod;
return ans;
}
const int MAXN = 1e7 + 5;
bool vis[MAXN];
LL n, mod, len, phi[MAXN], prime[MAXN], mu[MAXN], inv2, inv6;
gp_hash_table<LL, LL> Phi;
void seive1() {
phi[1] = 1, mu[1] = 1;
for (LL i = 2; i <= MAXN - 5; i++) {
if (!vis[i])
prime[++len] = i, phi[i] = i - 1, mu[i] = -1;
for (LL j = 1; j <= len && i * prime[j] <= MAXN - 5; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
mu[i * prime[j]] = -mu[i];
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
for (LL i = 1; i <= MAXN - 5; i++)
phi[i] = phi[i] * (i * i % mod) % mod, phi[i] += phi[i - 1], phi[i] %= mod;
}
LL sum2(LL x) { return x % mod * ((x + 1) % mod) % mod * ((2 * x + 1) % mod) % mod * inv6 % mod; }
LL calcphi(LL n) {
if (n <= MAXN - 5)
return phi[n];
if (Phi[n])
return Phi[n];
LL res = qkpow(n % mod * ((n + 1) % mod) % mod * inv2 % mod, 2, mod);
for (LL i = 2, j; i <= n; i = j + 1) {
j = n / (n / i);
res = res - ((sum2(j) - sum2(i - 1)) % mod * calcphi(n / i) % mod);
res %= mod;
}
return Phi[n] = (res + mod) % mod;
}
int T;
LL calc(LL n) {
LL ans = 0;
for (LL i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
ans = (ans + (((calcphi(j) - calcphi(i - 1)) % mod + mod) % mod *
qkpow(((1 + n / i) % mod) * (n / i % mod) % mod * inv2 % mod, 2, mod))) %
mod;
}
return ans;
}
int main() {
scanf("%lld %lld", &mod, &n);
inv2 = qkpow(2, mod - 2, mod);
inv6 = qkpow(6, mod - 2, mod);
seive1();
printf("%lld", (calc(n) + mod) % mod);
return 0;
}
「CQOI 2015」选数
给出$ N,K,L,H$ 。求
\[\sum_{i_1=L}^H\sum_{i_2=L}^H...\sum_{i_N=L}^H[\gcd(i_1,i_2...,i_N)=K]
\]
\[\sum_{i_1=\frac{L-1}{K}}^{\frac{H}{K}}\sum_{i_2=\frac{L-1}{K}}^{\frac{H}{K}}...\sum_{i_N=\frac{L-1}{K}}^{\frac{H}{K}}[\gcd(i_1,i_2...,i_N)=1]
\]
\[\sum_{i_1=\frac{L-1}{K}}^{\frac{H}{K}}\sum_{i_2=\frac{L-1}{K}}^{\frac{H}{K}}...\sum_{i_N=\frac{L-1}{K}}^{\frac{H}{K}}\sum_{D|\gcd(i_1,i_2...,i_N)}\mu(D)
\]
\[\sum_{D=1}^{\frac{H}{K}}\mu(D)\sum_{i_1=\frac{L-1}{K}}^{\frac{H}{K}}\sum_{i_2=\frac{L-1}{K}}^{\frac{H}{K}}...\sum_{i_N=\frac{L-1}{K}}^{\frac{H}{K}}[D|i_1][D|i_2]...[D|i_N]
\]
\[\sum_{D=1}^{\frac{H}{K}}\mu(D)(\lfloor\dfrac{H}{KD}\rfloor-\lfloor\dfrac{L-1}{KD}\rfloor)^N
\]
整除分块+杜教筛即可。
#include <cstdio>
#include <map>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL __int128
__int128 read() {
__int128 x = 0, f = 1;
char ch = getchar();
while (ch < '0' || ch > '9') {
if (ch == '-')
f = -1;
ch = getchar();
}
while (ch >= '0' && ch <= '9') {
x = x * 10 + ch - '0';
ch = getchar();
}
return x * f;
}
void write(__int128 x) {
if (x < 0) {
putchar('-');
x = -x;
}
if (x > 9)
write(x / 10);
putchar(x % 10 + '0');
}
const int MAXN = 1e6 + 5;
bool vis[MAXN];
const LL mod=1e9+7;
LL n, m, len, phi[MAXN], prime[MAXN], mu[MAXN];
map<LL,LL>Mu;
void seive1() {
phi[1] = 1, mu[1] = 1;
for (int i = 2; i <= MAXN-5; i++) {
if (!vis[i])
prime[++len] = i, phi[i] = i - 1, mu[i] = -1;
for (int j = 1; j <= len && i * prime[j] <= MAXN-5; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
mu[i * prime[j]] = - mu[i];
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
for (int i = 1; i <= MAXN-5; i++) phi[i] += phi[i - 1], mu[i]+=mu[i-1],mu[i]%=mod,mu[i]=(mu[i]+mod)%mod;
}
LL qkpow(LL x, LL power, LL mod) {
x %= mod;
LL ans = 1;
for (; power; power >>= 1, (x *= x) %= mod)
if (power & 1)
(ans *= x) %= mod;
return ans;
}
LL calcmu(LL n){
if(n<=MAXN-5)return mu[n];
if(Mu[n])return Mu[n];
LL res=1ll;
for(LL i=2,j;i<=n;i=j+1){
j=n/(n/i);
res=res-((j-i+1)%mod*calcmu(n/i)%mod);
res%=mod;
}
return Mu[n]=(res%mod+mod)%mod;
}
LL N,K,L,H;
LL calc(LL n){
LL ans=0;
for(LL i=1,j;i<=n;i=j+1){
j=min(n/(n/i),((L-1)/K)/i?((L-1)/K)/((((L-1)/K))/i):n/(n/i));
ans=(ans+(((calcmu(j)-calcmu(i-1))%mod+mod)%mod*qkpow((H/K)/i-((L-1)/K)/i,N,mod)%mod))%mod;
}
return ans%mod;
}
int main() {
seive1();
N=read(),K=read(),L=read(),H=read();
write((calc(H/K)+mod)%mod);
return 0;
}
