「例题」莫比乌斯反演
Luogu P4450 双亲数
Description
给定 \(n,m,d\)，qiu
\[\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]
\]
\(1\le n,m\le 10^6,1\le d\le min(n,m)\)
Solution
\[\begin{align*}
\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]
&=\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac m d}[\gcd(i,j)=1] \\
&=\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac m d}\sum_{k|gcd(i,j)}\mu(k) \\
&=\sum_{k=1}^n\sum_{i=1}^{\frac n {kd}}\sum_{j=1}^{\frac m {kd}}\mu(k) \\
&=\sum_{k=1}^n\left\lfloor \frac n {kd} \right\rfloor \left\lfloor \frac m {kd} \right\rfloor \mu(k)
\end{align*}
\]
然后就可以整除分块了
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 1e6 + 5;
int p[N], tot, mu[N];
bool flag[N];
void init(int n)
{
mu[1] = 1;
for(int i = 2; i <= n; i++)
{
if(!flag[i]) p[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * p[j] <= n; j++)
{
flag[i * p[j]] = 1;
if(i % p[j] == 0)
{
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for(int i = 1; i <= n; i++) mu[i] += mu[i - 1];
return;
}
ll solve(int n, int m, int d)
{
n /= d, m /= d;
if(n > m) n ^= m ^= n ^= m;
ll res = 0;
for(int l = 1, r; l <= n; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
res += 1ll * (mu[r] - mu[l - 1]) * (n / l) * (m / l);
}
return res;
}
int main()
{
init(1e6);
int n, m, d;;
read(n), read(m), read(d);
write(solve(n, m, d)), pc('\n');
return 0;
}
// A.S.
还有一道双倍经验 Luogu P3455 [POI2007]ZAP-Queries
Luogu P2522 Problem b
Description
\(n\) 组询问，每次询问给定 \(a,b,c,d,k\)，求
\[\sum_{i=a}^b\sum_{j=c}^d[\gcd(i,j)=k]
\]
\(1\le n,k \le 5\times 10^4,1\le a\le b\le 5\times 10^4,1\le c\le d\le 5\times 10^4\)
Solution
不难想到容斥，于是题目转化为求
\[\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k]
\]
就跟上一题一样了
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 5e4 + 5;
int p[N], tot, mu[N], sum[N];
bool flag[N];
void init(int n)
{
mu[1] = 1;
for(int i = 2; i <= n; i++)
{
if(!flag[i]) p[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * p[j] <= n; j++)
{
flag[i * p[j]] = 1;
if(i % p[j] == 0)
{
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for(int i = 1; i <= n; i++) sum[i] = sum[i - 1] + mu[i];
return;
}
int k;
ll calc(int n, int m)
{
n /= k, m /= k;
int d = min(n, m);
ll res = 0;
for(int l = 1, r; l <= d; l = r + 1)
{
r = min((n / (n / l)), (m / (m / l)));
res += 1ll * (sum[r] - sum[l - 1]) * (n / l) * (m / r);
}
return res;
}
int main()
{
init(N - 5);
int T; read(T);
while(T--)
{
int a, b, c, d;
read(a), read(b), read(c), read(d), read(k);
write(calc(b, d) - calc(a - 1, d) - calc(b, c - 1) + calc(a - 1, c - 1)), pc('\n');
}
return 0;
}
// A.S.
Luogu P3911 最小公倍数之和
Description
给定 \(a_1,a_2,\dots a_n\)，求
\[\sum_{i=1}^n\sum_{j=1}^n\text{lcm}(a_i,a_j)
\]
\(1\le n \le 5\times 10^4,1\le a_i \le 5 \times 10^4\)
Solution
设 \(m=max(a_i)\)，\(c_i\) 表示 \(\sum a\) 中有多少个 \(i\)。
\[\begin{align*}
\sum_{i=1}^n\sum_{j=1}^n\text{lcm}(a_i,a_j)
&=\sum_{i=1}^n\sum_{j=1}^n\dfrac{a_i·a_j}{\gcd(a_i,a_j)} \\
&=\sum_{i=1}^m\sum_{j=1}^m\dfrac{i·j}{\gcd(i,j)}·c_i·c_j \\
&=\sum_{d=1}^m\sum_{i=1}^m\sum_{j=1}^m[\gcd(i,j)=d]\dfrac{i·j}{d}·c_i·c_j \\
&=\sum_{d=1}^m\sum_{i=1}^{\frac m d}\sum_{j=1}^{\frac m d}[\gcd(i,j)=1]i·j·d·c_{id}·c_{jd} \\
&=\sum_{d=1}^md\sum_{i=1}^{\frac m d}\sum_{j=1}^{\frac m d}\sum_{k|\gcd(i,j)}\mu(k)·i·j·c_{idk}·c_{jdk} \\
&=\sum_{d=1}^md\sum_{k=1}^{\frac m d}\mu(k)·k^2\sum_{i=1}^{\frac m{dk}}\sum_{j=1}^{\frac m{dk}}i·j·c_{idk}·c_{jdk} \\
&=\sum_{dk=1}^mdk\sum_{k|dk}\mu(k)·k\sum_{i=1}^{\frac m{dk}}(i·c_{idk})^2 \\
\end{align*}
\]
令 \(T=dk\)
\[\sum_{i=1}^n\sum_{j=1}^n\text{lcm}(a_i,a_j)=
\sum_{T=1}^mT\sum_{k|T}\mu(k)·k\sum_{i=1}^{\frac m T}(i·c_{iT})^2 \\
\]
然后第二、三个求和可以预处理，复杂度是调和级数，最后再 \(O(m)\) 统计答案即可。
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 5e4 + 5;
int n, m, a[N];
int p[N], tot, mu[N];
bool flag[N];
ll c[N], sum[N], smu[N];
void init()
{
mu[1] = 1;
for(int i = 2; i <= m; i++)
{
if(!flag[i]) p[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * p[j] <= m; j++)
{
flag[i * p[j]] = 1;
if(i % p[j] == 0)
{
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for(int i = 1; i <= m; i++)
for(int j = i; j <= m; j += i)
smu[j] = smu[j] + mu[i] * i;
for(int i = 1; i <= n; i++) c[a[i]]++;
for(int i = 1; i <= m; i++)
for(int j = 1; j <= m / i; j++)
sum[i] = sum[i] + j * c[i * j];
for(int i = 1; i <= m; i++) sum[i] = sum[i] * sum[i];
return;
}
signed main()
{
read(n);
for(int i = 1; i <= n; i++) read(a[i]), m = max(m, a[i]);
init();
ll ans = 0;
for(int i = 1; i <= m; i++)
ans += i * smu[i] * sum[i];
write(ans), pc('\n');
return 0;
}
// A.S.
AT5200 LCMs
Description
给定 \(a_1,a_2,\dots a_n\)，求
\[\sum_{i=1}^n\sum_{j=i+1}^n\text{lcm}(a_i,a_j)
\]
对 \(998244353\) 取模
\(1\le n \le 2\times 10^5,1\le a_i \le 10^5\)
Solution
\[\sum_{i=1}^n\sum_{j=i+1}^n\text{lcm}(a_i,a_j)
=\dfrac{\sum_{i=1}^n\sum_{j=1}^n\text{lcm}(a_i,a_j)-sum_{i=1}^na_i}{2}
\]
剩下的就跟上一题一样了。
注意有 \(\mu\) 可能会有负数，所以还是老老实实取模吧，不然看不出来要调好久 qaq
COde
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 2e5 + 5;
const int M = 1e6 + 5;
const int mod = 998244353;
ll add(ll x) {return x < mod ? x : x - mod;}
ll sub(ll x) {return x < 0 ? x : x + mod;}
ll qpow(ll a, int b)
{
ll res = 1;
while(b)
{
if(b & 1) res = res * a % mod;
a = a * a % mod, b >>= 1;
}
return res;
}
int n, m, a[N];
int p[M], tot, mu[M];
bool flag[M];
ll c[M], sum[M], smu[M];
void init()
{
mu[1] = 1;
for(int i = 2; i <= m; i++)
{
if(!flag[i]) p[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * p[j] <= m; j++)
{
flag[i * p[j]] = 1;
if(i % p[j] == 0)
{
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for(int i = 1; i <= m; i++)
for(int j = i; j <= m; j += i)
smu[j] = add((smu[j] + mu[i] * i % mod) % mod + mod);
for(int i = 1; i <= n; i++) c[a[i]]++;
for(int i = 1; i <= m; i++)
for(int j = 1; j <= m / i; j++)
sum[i] = add(sum[i] + j * c[i * j] % mod);
for(int i = 1; i <= m; i++) sum[i] = sum[i] * sum[i] % mod;
return;
}
signed main()
{
read(n);
for(int i = 1; i <= n; i++) read(a[i]), m = max(m, a[i]);
init();
ll ans = 0;
for(int i = 1; i <= m; i++)
ans = ((ans + i * smu[i] % mod * sum[i] % mod) % mod + mod) % mod;
for(int i = 1; i <= n; i++) ans = ((ans - a[i]) % mod + mod) % mod;
ans = ans * qpow(2, mod - 2) % mod;
write(ans), pc('\n');
return 0;
}
// A.S.
Luogu P5221 Product
Description
给定 \(n\)，求
\[\prod_{i=1}^n\prod_{j=1}^n\dfrac{\text{lcm}(i,j)}{\gcd(i,j)}\quad(\bmod 104857601)
\]
\(1\le n \le 10^6\)
Solution
\[\begin{align*}
\prod_{i=1}^n\prod_{j=1}^n\dfrac{\text{lcm}(i,j)}{\gcd(i,j)}
&=\prod_{i=1}^n\prod_{j=1}^n\dfrac{i*j}{\gcd(i,j)^2} \\
&=\dfrac{(n!)^2}{\prod_{i=1}^n\prod_{j=1}^n\gcd(i,j)^2} \\
\prod_{i=1}^n\prod_{j=1}^n\gcd(i,j)^2
&=(\prod_{i=1}^n\prod_{j=1}^n\gcd(i,j))^2 \\
\prod_{i=1}^n\prod_{j=1}^n\gcd(i,j)
&=\prod_{d=1}^n\prod_{i=1}^n\prod_{j=1}^n[\gcd(i,j)=d]d \\
&=\prod_{d=1}^nd^{\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d]} \\
\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d]
&=\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac n d}[\gcd(i,j)=1] \\
&=1+2\sum_{i=1}^n\varphi(i)
\end{align*}
\]
预处理一下 \(\varphi\) 的前缀和即可，设 \(sum(i)=\sum_{j=1}^i\varphi(j)\)。
\[\prod_{i=1}^n\prod_{j=1}^n\dfrac{\text{lcm}(i,j)}{\gcd(i,j)}
=\dfrac{(n!)^2}{(\prod_{d=1}^n d^{2*sum(\frac n d)+1})^2}
\]
根据欧拉定理，可以将 \(2*sum(\frac n d)+1\bmod (p-1)\)
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 1e6 + 5;
const int p = 104857601;
ll qpow(ll a, int b)
{
ll res = 1;
while(b)
{
if(b & 1) res = res * a % p;
a = a * a % p, b >>= 1;
}
return res;
}
int pri[N], tot, phi[N];
bool flag[N];
void getphi(int n)
{
phi[1] = 1;
for(int i = 2; i <= n; i++)
{
if(!flag[i]) pri[++tot] = i, phi[i] = i - 1;
for(int j = 1; j <= tot && i * pri[j] <= n; j++)
{
flag[i * pri[j]] = 1;
if(i % pri[j] == 0)
{
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
phi[i * pri[j]] = phi[i] * (pri[j] - 1);
}
}
for(int i = 1; i <= n; i++) phi[i] = (phi[i - 1] + phi[i]) % (p - 1);
return;
}
int main()
{
int n; read(n);
getphi(n);
int ans = 1;
for(int i = 1; i <= n; i++) ans = 1ll * ans * i % p;
ans = qpow(ans, 2 * n);
int tmp = 1;
for(int i = 1; i <= n; i++)
tmp = 1ll * tmp * qpow(i, (2 * phi[n / i] - 1) % (p - 1)) % p;
ans = 1ll * ans * qpow(1ll * tmp * tmp % p, p - 2) % p;
write(ans), pc('\n');
return 0;
}
// A.S.
Luogu P4917 天守阁的地板
Description
\(T\) 组数据，每组数据给定一个 \(n\)，求
\[\prod_{i=1}^n\prod_{j=1}^n\dfrac{\text{lcm}(i,j)}{\gcd(i,j)}\quad(\bmod 19260817)
\]
\(1\le T \le 1000,1\le n \le 10^6\)
有没有发现跟上一题一样（强烈谴责题面
Solution
本题与上一题唯一的区别就是最后计算答案时需要整除分块，因为 \(O(Tn)\) 过不去，而我上一题偷懒没写（
Code
int calc(int n)
{
int ans = qpow(fac[n], n << 1), res = 1;
for(int l = 1, r; l <= n; l = r + 1)
{
r = n / (n / l);
res = 1ll * res * qpow(1ll * fac[r] * qpow(fac[l - 1], p - 2) % p, 2 * phi[n / l] - 1) % p;
}
return 1ll * ans * qpow(1ll * res * res % p, p - 2) % p;
}
Luogu P3327 约数个数和
Description
设 \(d(x)\) 为 \(x\) 的约数个数，给定 \(n,m\)，求
\[\sum_{i=1}^n\sum_{j=1}^md(i*j)
\]
\(T\) 组数据
\(1\le T,n,m\le 50000\)
Solution
首先有个性质
\[d(x*y)=\sum_{i|x}\sum_{j|y}[\gcd(i,j)=1]
\]
所以原式为
\[\sum_{x=1}^n\sum_{y=1}^m\sum_{i|x}\sum_{j|y}[\gcd(i,j)=1]
=\sum_{i=1}^n\sum_{j=1}^m\left\lfloor\dfrac n i\right\rfloor \left\lfloor\dfrac m j\right\rfloor[\gcd(i,j)=1]
\]
设
\[f(x)=\sum_{i=1}^n\sum_{j=1}^m\left\lfloor\dfrac n i\right\rfloor\left\lfloor\dfrac m j\right\rfloor[\gcd(i,j)=x] \\
g(x)=\sum_{x|d}f(d)
\]
那么
\[\begin{align*}
g(x)&=\sum_{i=1}^n\sum_{j=1}^m\left\lfloor\dfrac n i\right\rfloor\left\lfloor\dfrac m j\right\rfloor[x|\gcd(i,j)] \\
&=\sum_{i=1}^{\frac n x}\sum_{j=1}^{\frac m x}\left\lfloor\dfrac n {ix}\right\rfloor\left\lfloor\dfrac m {jx}\right\rfloor
\end{align*}
\]
因为
\[f(n)=\sum_{n|d}\mu(\frac n d)g(d)
\]
所以
\[\begin{align*}
f(1)&=\sum_{1|d}\mu(d)g(d) \\
&=\sum_{i=1}^n\mu(i)g(i)
\end{align*}
\]
考虑如何求 \(g(x)\)，设 \(s(x)=\sum_{i=1}^n\lfloor\frac n i\rfloor\)
\[g(x)=s(\lfloor\frac n i\rfloor)*s(\lfloor\frac m j\rfloor)
\]
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 5e4 + 5;
int p[N], tot;
ll mu[N], sum[N];
bool flag[N];
void init(int n)
{
mu[1] = 1;
for(int i = 2; i <= n; i++)
{
if(!flag[i]) p[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * p[j] <= n; j++)
{
flag[i * p[j]] = 1;
if(i % p[j] == 0)
{
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for(int i = 1; i <= n; i++) mu[i] += mu[i - 1];
for(int i = 1; i <= n; i++)
{
for(int l = 1, r; l <= i; l = r + 1)
{
r = i / (i / l);
sum[i] += 1ll * (r - l + 1) * (i / l);
}
}
return;
}
ll calc(int n, int m)
{
ll res = 0;
int mn = min(n, m);
for(int l = 1, r; l <= mn; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
res += (mu[r] - mu[l - 1]) * sum[n / l] * sum[m / l];
}
return res;
}
int main()
{
init(5e4);
int T; read(T);
while(T--)
{
int n, m;
read(n), read(m);
write(calc(n, m)), pc('\n');
}
return 0;
}
// A.S.
Luogu P1829 Crash的数字表格 / JZPTAB
Description
给定 \(n,m\)，求
\[\sum_{i=1}^n\sum_{j=1}^m\text{lcm}(i,j)\quad(\bmod 20101009)
\]
\(1\le n,m\le 10^7\)
Solution
默认 \(n\le m\)，除法向下取整
\[\begin{align*}
\sum_{i=1}^n\sum_{j=1}^m\text{lcm}(i,j)
&=\sum_{i=1}^n\sum_{j=1}^m\dfrac{i*j}{\gcd(i,j)} \\
&=\sum_{i=1}^n\sum_{j=1}^m\sum_{d=1}^{n}[\gcd(i,j)=d]\dfrac{i*j}{d}\\
&=\sum_{d=1}^nd\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac m d}[\gcd(i,j)=1]i*j \\
\end{align*}
\]
设
\[\begin{align*}
sum(n,m)&=\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=1]i*j \\
&=\sum_{i=1}^n\sum_{j=1}^m\sum_{d|\gcd(i,j)}\mu(d)*i*j \\
&=\sum_{d=1}^n\mu(d)*d^2\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac m d}i*j \\
\end{align*}
\]
前面预处理前缀和即可，后面就是两个等差数列的和的乘积
记
\[getsum(n,m)=\dfrac{n(n+1)} 2\times \dfrac{m(m+1)}2
\]
则
\[sum(n,m)=\sum_{d=1}^n\mu(d)*d^2*getsum(\frac n d,\frac m d)
\]
可以整除分块，然后原式为
\[\sum_{d=1}^nd*sum(\frac n d,\frac m d)
\]
发现又可以整除分块
Code
#include <bits/stdc++.h>
#define int long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 1e7 + 5;
const int mod = 20101009;
int p[N], tot, mu[N];
int sum[N];
bool flag[N];
void init(int n, int m)
{
mu[1] = 1;
for(int i = 2; i <= n; i++)
{
if(!flag[i]) p[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * p[j] <= n; j++)
{
flag[i * p[j]] = 1;
if(i % p[j] == 0)
{
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for(int i = 1; i <= n; i++) sum[i] = (sum[i - 1] + i * i % mod * (mu[i] + mod) % mod) % mod;
return;
}
int getsum(int n, int m)
{
return (n * (n + 1) / 2 % mod) * (m * (m + 1) / 2 % mod) % mod;
}
int calc(int n, int m)
{
int res = 0;
for(int l = 1, r; l <= n; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
res = (res + (sum[r] - sum[l - 1] + mod) * getsum(n / l, m / l) % mod) % mod;
}
return res;
}
int solve(int n, int m)
{
int res = 0;
for(int l = 1, r; l <= n; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
res = (res + ((l + r) * (r - l + 1) / 2 % mod) * calc(n / l, m / l) % mod) % mod;
}
return res;
}
signed main()
{
int n, m;
read(n), read(m);
if(n > m) n ^= m ^= n ^= m;
init(n, m);
write(solve(n, m)), pc('\n');
return 0;
}
// A.S.
Luogu P3768 简单的数学题
Description
给定 \(n,p\)，求
\[\left( \sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j) \right)\bmod p
\]
\(n\le 10^{10},5\times 10^8\le p \le 1.1\times 10^9\)
Solution
\[\begin{align*}
\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)
&=\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d]ij \\
&=\sum_{d=1}^nd^3\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac n d}[\gcd(i,j)=1]ij \\
&=\sum_{d=1}^nd^3\sum_{k=1}^{\frac n d}\mu(k)k^2(\sum_{i=1}^{\frac n {dk}}i)^2
\end{align*}
\]
设 \(sum(n)=\sum_{i=1}^ni\)
令 \(T=dk\)
\[\sum_{T=1}^nsum(\frac n T)^2\sum_{d|T}d^3(\frac T d)^2\mu(\frac T d) \\
\sum_{T=1}^nsum(\frac n T)^2T^2\sum_{d|T}d\mu(\frac T d)
\]
现在只需要筛出来后面的那一部分
\[id*\mu=\varphi \\
\sum_{d|T}d\mu(\frac T d)=\varphi(T) \\
T^2\sum_{d|T}d\mu(\frac T d)=T^2\varphi(T)
\]
设
\[f(i)=i^2\varphi(i),S(n)=\sum_{i=1}^nf(i)
\]
看一下杜教筛套路式子
\[g(1)S(n)=\sum_{i=1}^n(g*f)(i)-\sum_{i=2}^ng(i)S(\frac n i)
\]
考虑 \(g\) 应该是什么
因为 \(\varphi*1=id\)，令 \(g(i)=i^2\)
\[S(n)=\sum_{i=1}^n(g*f)(i)-\sum_{i=2}^ng(i)S(\frac n i) \\
\begin{align*}
(g*f)(i)&=\sum_{d|i}f(d)g(\frac i d) \\
&=\sum_{d|i}d^2\varphi(d)(\frac i d)^2 \\
&=i^2\sum_{d|i}\varphi(d) \\
&=i^3
\end{align*}
\]
所以
\[S(n)=\sum_{i=1}^ni^3-\sum_{i=2}^ni^2S(\frac n i)
\]
其中
\[\sum_{i=1}^ni^3=(\sum_{i=1}^ni)^2
\]
因此原式
\[\sum_{T=1}^nsum(\frac n T)^2f(T)
\]
就可以整除分块加杜教筛解决了
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 4e6 + 5;
ll p, n;
ll pri[N], tot, phi[N];
ll sum[N], inv6;
bool flag[N];
map <ll, ll> mp;
ll add(ll x) {return x < p ? x : x - p;}
ll sub(ll x) {return x < 0 ? x + p : x;}
ll s2(ll x) {x %= p; return x * (x + 1) % p * (2 * x + 1) % p * inv6 % p;}
ll s3(ll x) {x %= p; return (x * (x + 1) / 2 % p) * (x * (x + 1) / 2 % p) % p;}
ll qpow(ll a, int b)
{
ll res = 1;
while(b)
{
if(b & 1) res = res * a % p;
a = a * a % p, b >>= 1;
}
return res;
}
void init(int n)
{
phi[1] = 1;
for(int i = 2; i < n; i++)
{
if(!flag[i]) pri[++tot] = i, phi[i] = i - 1;
for(int j = 1; j <= tot && i * pri[j] < n; j++)
{
flag[i * pri[j]] = 1;
if(i % pri[j] == 0)
{
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
phi[i * pri[j]] = phi[i] * (pri[j] - 1);
}
}
for(int i = 1; i < n; i++) sum[i] = add(sum[i - 1] + 1ll * i * i % p * phi[i] % p);
inv6 = qpow(6, p - 2);
return;
}
ll calc(ll n)
{
if(n < N) return sum[n];
if(mp[n]) return mp[n];
ll res = s3(n);
for(ll l = 2, r; l <= n; l = r + 1)
{
r = n / (n / l);
res = sub(res - sub(s2(r) - s2(l - 1)) * calc(n / l) % p);
}
return mp[n] = res;
}
ll solve(ll n)
{
ll res = 0;
for(ll l = 1, r; l <= n; l = r + 1)
{
r = n / (n / l);
res = add(res + sub(calc(r) - calc(l - 1)) * s3(n / l) % p);
}
return res;
}
int main()
{
read(p), read(n);
init(N);
write(solve(n)), pc('\n');
return 0;
}
// A.S.
Luogu P3704 数字表格
Description
\(T\) 组数据，每组数据给定 \(n,m\)，求
\[\prod_{i=1}^n\prod_{j=1}^mf_{\gcd(i,j)}\quad (\bmod 10^9+7)
\]
其中 \(f_i\) 表示斐波那契数列的第 \(i\) 项。
\(1\le T \le 10^3,1\le n,m\le 10^6\)
Solution
\[\begin{align*}
\prod_{i=1}^n\prod_{j=1}^mf_{\gcd(i,j)}
&=\prod_{i=1}^n\prod_{j=1}^m\prod_{d=1}^n[\gcd(i,j)=d]f_{d} \\
&=\prod_{d=1}^nf_d^{\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]} \\
\end{align*}
\]
发现指数部分非常经典
\[\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]=\sum_{k=1}^{\frac n d}\mu(k)\left\lfloor \frac n{dk} \right\rfloor\left\lfloor \frac m{dk} \right\rfloor
\]
带回原式
\[\prod_{d=1}^nf_d^{\sum_{k=1}^{\frac n d}\mu(k)\left\lfloor \frac n{dk} \right\rfloor\left\lfloor \frac m{dk} \right\rfloor}
=\prod_{d=1}^n\prod_{k=1}^{\frac n d}f_d^{\mu(k)\left\lfloor \frac n{dk} \right\rfloor\left\lfloor \frac m{dk} \right\rfloor}
\]
令 \(T=dk\)
\[\prod_{T=1}^n\prod_{d|T}f_d^{\mu(\frac T d)\left\lfloor \frac n{dk} \right\rfloor\left\lfloor \frac m{dk} \right\rfloor}=\prod_{T=1}^n\left(\prod_{d|T}f_d^{\mu(\frac T d)}\right)^{\left\lfloor \frac n{dk} \right\rfloor\left\lfloor \frac m{dk} \right\rfloor}
\]
括号外面的可以整除分块，里面的可以直接预处理，复杂度是调和级数
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 1e6 + 5;
const int mod = 1e9 + 7;
int p[N], tot, mu[N];
bool flag[N];
ll f[N], inv[N], pro[N];
ll add(ll x) {return x < mod ? x : x - mod;}
ll qpow(ll a, int b)
{
ll res = 1;
a %= mod;
while(b)
{
if(b & 1) res = res * a % mod;
a = a * a % mod, b >>= 1;
}
return res;
}
void init(int n)
{
mu[1] = 1;
f[1] = 1ll, inv[1] = 1ll, pro[0] = pro[1] = 1ll;
for(int i = 2; i <= n; i++)
{
if(!flag[i]) p[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * p[j] <= n; j++)
{
flag[i * p[j]] = 1;
if(i % p[j] == 0)
{
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
f[i] = add(f[i - 1] + f[i - 2]);
inv[i] = qpow(f[i], mod - 2);
pro[i] = 1ll;
}
for(int i = 1; i <= n; i++)
for(int j = i; j <= n; j += i)
if(mu[j / i]) pro[j] = pro[j] * (mu[j / i] == 1 ? f[i] : inv[i]) % mod;
for(int i = 2; i <= n; i++) pro[i] = pro[i] * pro[i - 1] % mod;
return;
}
ll solve(int n, int m)
{
ll res = 1;
for(int l = 1, r; l <= n; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
res = res * qpow(pro[r] * qpow(pro[l - 1], mod - 2) % mod, 1ll * (n / l) * (m / l) % (mod - 1)) % mod;
}
return res;
}
int main()
{
init(1e6);
int T; read(T);
while(T--)
{
int n, m; read(n), read(m);
if(n > m) n ^= m ^= n ^= m;
write(solve(n, m)), pc('\n');
}
return 0;
}
// A.S.
Luogu P2257 YY的GCD
Description
\(T\) 组数据，每组给定 \(n,m\)，求
\[\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)\in prime]
\]
\(T=10^4,n,m\le 10^7\)
Solution
默认 \(n\le m\)，\(\mathbb{P}\) 表示质数
\[\begin{align*}
\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)\in \mathbb{P}]
&=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\in \mathbb{P}}[\gcd(i,j)=d] \\
&=\sum_{d\in \mathbb{P}}\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac m d}[\gcd(i,j)=1] \\
&=\sum_{d\in \mathbb{P}}\sum_{k=1}^{\frac n d}\mu(k)\left\lfloor \frac n{dk} \right\rfloor\left\lfloor \frac m{dk} \right\rfloor \\
\end{align*}
\]
令 \(T=dk\)
\[\sum_{T=1}^n\sum_{d|T,d\in \mathbb{P}}\mu(\frac T d)\left\lfloor \frac n{T} \right\rfloor\left\lfloor \frac m{T} \right\rfloor
=\sum_{T=1}^n\left\lfloor \frac n{T} \right\rfloor\left\lfloor \frac m{T} \right\rfloor\sum_{d|T,d\in \mathbb{P}}\mu(\frac T d)
\]
预处理一下后面的东西就可以整除分块了
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 1e7 + 5;
int p[N], tot, mu[N], sum[N];
bool flag[N];
void init(int n)
{
mu[1] = 1;
for(int i = 2; i <= n; i++)
{
if(!flag[i]) p[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * p[j] <= n; j++)
{
flag[i * p[j]] = 1;
if(i % p[j] == 0)
{
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for(int i = 1; i <= tot; i++)
for(int j = p[i]; j <= n; j += p[i])
sum[j] += mu[j / p[i]];
for(int i = 1; i <= n; i++) sum[i] += sum[i - 1];
return;
}
ll solve(int n, int m)
{
if(n > m) n ^= m ^= n ^= m;
ll res = 0;
for(int l = 1, r; l <= n; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
res += 1ll * (sum[r] - sum[l - 1]) * (n / l) * (m / l);
}
return res;
}
int main()
{
init(1e7);
int T; read(T);
while(T--)
{
int n, m;
read(n), read(m);
write(solve(n, m)), pc('\n');
}
return 0;
}
// A.S.
$$A\ drop\ of\ tear\ blurs\ memories\ of\ the\ past.$$
