Updated on 2020.8.6
巨幅更新,对积性函数和狄利克雷卷积部分进行重构。
新增对一类特殊求和式的讲解。
Updated on 2020.9.9
添加了几个例题。
前置知识
小碎骨
艾佛森括号 ([P] = egin{cases} 1 & ext{If P is true}\ 0 & ext{Otherwise} end{cases})
此处 (P) 是一个可真可假的命题。
引理1
证明
数论分块
内容独立了出来,详细内容见 数论分块 - Luckyblock
对于一类含有(leftlfloorfrac{n}{i} ight floor)的求和式 ((n) 为常数),由于(leftlfloorfrac{n}{i} ight floor)单调不增,故存在多个区间([l,r]), 使得(leftlfloorfrac{n}{i} ight floor = leftlfloorfrac{n}{j} ight floor(i,jin [l,r])) 成立。
对于任意一个(i),最大的满足上式的 (j=leftlfloor{dfrac{n}{leftlfloor{dfrac{n}{i}} ight floor}} ight floor)
积性函数
定义
若 (gcd(x,y) = 1) 且 (f(xy)=f(x)f(y)),则 (f(n)) 为积性函数。
性质
若 (f(x)),(g(x))均为积性函数,则以下函数也为积性函数:
常见积性函数
- 单位函数 (e(n) = [n = 1])。
- 幂函数 (operatorname{Id}_{k}(n) = n^k), (operatorname{id}_1(n)) 通常简记为(operatorname{id}(n))。
- 常数函数 (1(n) = 1)。
- 因数个数 (operatorname{d}(n) = sumlimits_{dmid n} 1)。
- 除数函数 (sigma_{k}(n) = sumlimits_{dmid n} d^k)。
(k=1) 时为因数和函数,通常简记为 (sigma(n))。
(k=0) 时为因数个数函数 (sigma_{0}(n))。 - 欧拉函数 (varphi(n) = sumlimits_{i=1}^{n} [gcd(i,n) = 1])。
- 莫比乌斯函数 (mu(n) = egin{cases}1 &n=1\0 &n ext{含有平方因子}\(-1)^k &k ext{为} n ext{的本质不同质因子个数} end{cases})
不是很懂上面写的什么玩意?
不用深究,有个印象继续往下看就好。
莫比乌斯函数
定义
(mu) 为莫比乌斯函数,定义为
解释
令 (n = prodlimits_{i=1}^{k} p_{i}^{c_i}),其中(p_i)为质因子,(c_ige 1)。
- (n=1)时,(mu (n) = 1)。
- (n
ot ={1})时 ,
- (exist iin [1,k], c_i > 1) 时,(mu (n) = 0)。
当某质因子出现次数大于(1)时,(mu (n) = 0)。 - (forall iin [1,k], c_i = 1) 时,(mu (n) = (-1)^k)。
当每个质因子只出现一次时,即(n = prodlimits_{i=1}^{k}p_i),({p_i})中元素唯一。
(mu (n) = (-1)^k),此处(k)为质因子的种类数。
- (exist iin [1,k], c_i > 1) 时,(mu (n) = 0)。
性质
莫比乌斯函数是积性函数,且具有以下性质
证明,设 (n = prodlimits_{i=1}^{k}{p_i^{c_i}}, n' = prodlimits_{i=1}^{k}{p_i})。
- 根据莫比乌斯函数定义,则有:(sumlimits_{dmid n}{mu(d)} = sumlimits_{dmid n'}{mu(d)})。
- 若 (n') 的某因子 (d),有 (mu (d) = (-1)^i),则它由 (i) 个 本质不同的质因子组成。
由于质因子总数为 (k),满足上式的因子数为 (C_{k}^{i})。 - 对于原求和式,转为枚举 (mu(d)) 的值。
则 (sumlimits_{dmid n'}{mu(d)} = sumlimits_{i=0}^{k}{C_{k}^{i} imes (-1)^i} = sumlimits_{i=0}^{k}{C_{k}^{i} imes (-1)^i imes 1^{k-i}})
根据二项式定理,上式 (= (1+(-1))^k),
易知该式在 (k=0),即 (n=0) 时为 (1),否则为 (0)。
反演常用结论
一个反演常用结论:
证明 1:
设 (n = gcd(i,j)),则右(= sumlimits_{dmid n} {mu (d)} = [n = 1] = [gcd(i,j) = 1]=) 左。
证明 2:
暴力展开:([gcd(i,j) = 1] = e(gcd(i,j)) = sumlimits_{dmid gcd(i,j)}mu(d))
线性筛求莫比乌斯函数
(mu) 为积性函数,因此可以线性筛莫比乌斯函数。
int pnum, mu[kMaxn], p[kMaxn];
bool vis[kMaxn];
void Euler(int lim_) {
vis[1] = true, mu[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = - mu[i];
}
}
}
狄利克雷(Dirichlet)卷积
建议阅读 算法学习笔记(35): 狄利克雷卷积 By: Pecco。
定义两个数论函数 (f,g) 的狄利克雷卷积为
性质
- 显然满足 交换律,结合律,分配律。
- (f ast g = g ast f)
- ((f ast g) ast h = fast (gast h))
- (fast (g+h) = fast g + fast h)
- (e) 为狄利克雷卷积的单位元,有((fast e)(n) = f(n))。
- 若 (f,g) 为积性函数,则 (fast g) 为积性函数。
关于单位元 (e)
有:
证明
- 对于([dfrac{n}{d} = 1]),当且仅当 (dfrac{n}{d}=1),即 (d=n) 时为 (1),否则为(0)。
- 则当 (d=n) 时,(f(d)left[dfrac{n}{d} = 1
ight] = f(n))。
当 (d ot ={n}) 时,(f(d)left[dfrac{n}{d} = 1 ight] = 0)。
综上,((fast e)(n) = sumlimits_{dmid n} f(d)left[dfrac{n}{d} = 1
ight] = f(n)),满足单位元的性质。
则 (e = mu ast 1) 成立。
除数函数与幂函数
幂函数 (operatorname{Id}_{k}(n) = n^k)。
除数函数 (sigma_{k}(n) = sumlimits_{dmid n} d^k)。
显然有:
当 (k=0) 时,(operatorname{Id_0}=1),(sigma_0) 为因数个数函数,有:
欧拉函数与恒等函数
对于一式,当 (n=p^m) 时((p) 为质数),有:
(p^i) 的因子有 (p^{i-1}) 个,为 (1sim p^{i-1}),故 (varphi(p^i) = p^i-p^{i-1})。
又 ((varphi ast 1)(n)) 为积性函数,则对于任意正整数 (n),有:
得证。
对于 2 式,在 1 式基础上两侧同时 (ast mu) 即得。
左侧变为 (varphi ast 1ast mu = varphi ast e = varphi)。
计算
考虑枚举 (n) 的因子,将贡献累加即可。
显然可以使用埃氏筛筛出所有前缀狄利克雷卷积,复杂度 (O(nklog n)),其中 (k) 是计算函数值的复杂度。
莫比乌斯反演
反演是个啥?反演。
公式
设(f(n),g(n))为两个数论函数。
如果有
那么有
证明
方法一:运用卷积。
原问题为:已知 (f = gast 1),证明 (g = fast mu)。
易知如下转化:(fast mu = gast 1 ast mu Longrightarrow fast mu = gast e = g)。
方法二:对原式进行数论变换。
-
用 (sumlimits_{dmid n}g(d)) 替换(fleft(dfrac{n}{d} ight))。
[sum_{dmid n}{mu(d)sum_{kmid frac{n}{d}}g(k)} ] -
变换求和顺序。
[sum_{kmid n}g(k)sum_{dmid frac{n}{k}}{mu(d)} ] -
依据 (sumlimits_{dmid n}{mu(d)} = [n=1]),仅当 (dfrac{n}{k} = 1) 时,(sumlimits_{dmid frac{n}{k}}{mu(d)} = 1),否则为 (0)。
此时(k=n),故上式等价于 (sumlimits_{kmid n} {[n=k]cdot g(k)} = g(n))。
举例
从 狄利克雷(Dirichlet)卷积 部分可以知道一些积性函数的关系。
尝试对它们进行反演:
关于一类求和式
一般套路
考虑构造出函数 (g),满足下列形式:
则求和式变为:
考虑算术基本定理,发现若满足 (dmid gcd (i,j)),则 (dmid i) 且 (dmid j) 成立。
考虑 (g(d)) 在何时做出贡献,调整枚举顺序:
(sumlimits_{i=1}^{n}[dmid i]) 等价于 (1sim n) 中 (d) 的倍数的个数,则上式等价于:
数论分块求解即可。
例 1
发现此题的 (f) 等价于 (operatorname{Id}),则上式等价于:
例 2
感悟
卷点什么东西,把 (g) 卷出来。
(g) 不一定是特殊意义的函数。
例题
[HAOI2011] Problem b
(n) 组询问,每次给定参数 (a,b,c,d,k),求:
[sumlimits_{x=a}^{b}sumlimits_{y=c}^{d}[gcd(x,y) = k] ](1le n,k,a,b,c,dle 5 imes 10^4),(ale b,cle d)。
令 (f(n,m) = sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j) = k])。
根据容斥原理,则原式等价于 (f(b,d) - f(a-1,d) - f(b,d-1) + f(a-1,d-1))。
(f) 变成了上述一类求和式的形式,考虑化简 (f)。
易知原式等价于
代入反演常用结论 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}),上式化为
变换求和顺序,先枚举(dmid gcd(i,j)),可得
对于上式后半的解释:当(dmid i) 且 (dmid j) 时,(dmid gcd(i,j))。
易知(1sim leftlfloordfrac{n}{k} ight floor) 中 (d) 的倍数有 (leftlfloordfrac{leftlfloordfrac{n}{k} ight floor}{d} ight floor = leftlfloordfrac{n}{kd} ight floor) 个(由引理 1 可知),原式变为
预处理 (mu) 后,显然可以数论分块求解,复杂度(O(n + Tsqrt{n}))。
代码
//知识点:莫比乌斯反演
/*
//By:Luckyblock
*/
#include <cstdio>
#include <cctype>
#include <algorithm>
#define ll long long
const int MARX = 6e4 + 10;
//=============================================================
int N, a, b, c, d, k;
int cnt, Mobius[MARX], Prime[MARX], Sum[MARX];
bool vis[MARX];
//=============================================================
inline int read()
{
int f = 1, w = 0; char ch = getchar();
for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = -1;
for(; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void GetMobius()
{
Mobius[1] = 1;
int MAX = MARX - 10;
for(int i = 2; i <= MAX; i ++)
{
if(! vis[i]) Mobius[i] = - 1, Prime[++ cnt] = i;
for(int j = 1; j <= cnt && i * Prime[j] <= MAX; j ++)
{
vis[i * Prime[j]] = true;
if(i % Prime[j] == 0) break;
Mobius[i * Prime[j]] = - Mobius[i];
}
}
for(int i = 1; i <= MAX; i ++) Sum[i] = Sum[i - 1] + Mobius[i];
}
ll Calc(int x, int y)
{
ll ans = 0ll; int max = std :: min(x, y);
for(int l = 1, r; l <= max; l = r + 1)
r = std :: min(x / (x / l), y / (y / l)),
ans += (1ll * x / (1ll * l * k)) * (1ll * y / (1ll * l * k)) * (Sum[r] - Sum[l - 1]);
return ans;
}
ll Solve()
{
a = read(), b = read(), c = read(), d = read(), k = read();
return Calc(b, d) - Calc(b, c - 1) - Calc(a - 1, d) + Calc(a - 1, c - 1);
}
//=============================================================
int main()
{
N = read(); GetMobius();
while(N --) printf("%lld
", Solve());
return 0;
}
[国家集训队]Crash的数字表格
给定 (n,m) , 求:
[sum_{i=1}^nsum_{j=1}^{m} operatorname{lcm}(i,j)mod 20101009 ](1le n,mle 10^7)。
易知原式等价于:
考虑枚举 (gcd(i,j)),设枚举量为 (d)。
则 (d=gcd(i,j)) 的充要条件是满足 (d|i, d|j) 且 (gcd(dfrac{i}{d},dfrac{j}{d}) = 1),则原式等价于:
先枚举 (d),则原式等价于:
这个 (d) 很烦人,把 (i,j) 中的 (d) 提出来,变为枚举 (frac{i}{d}),(frac{j}{d})。
消去 (dmid i),(dmid j) 的限定条件,则原式等价于:
单独考虑后半部分,设 (f(x,y) = sumlimits_{i=1}^{x} sumlimits_{j=1}^{y}[gcd(i,j)=1]ij)。
发现 (f(x,y)) 的形式与 [HAOI2011] Problem b 中的式子类似,代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}) 套路一波:
前半部分 (sumlimits_{d=1}mu(d) d^2),可以考虑筛出 (mu(d)) 后求前缀和。
后半部分是等差数列乘等差数列的形式,设 (g(p,q) = sumlimits_{i=1}^{p} sumlimits_{j=1}^{q}ij),(g_{p,q}) 可以通过下式 (O(1)) 计算:
则对于 (f(x,y)),有:
数论分块求解即可。
再看回原式,原式等价于:
又是一个可以数论分块求解的形式。
线性筛预处理后 数论分块套数论分块,复杂度 (O(n + m)),瓶颈是线性筛。
一些注意的点
处理时会出现求平方的运算,需要特别注意取模问题,ll 都会爆,被坑惨了。
在预处理前缀和的这个地方:
sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
注意先对平方取模,在乘上 (mu),否则就会爆掉。
以及可以仅令 (mu + mod)。
以及这个地方:
int g(int n_, int m_) {
return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
平方计算,注意随时取模。
代码
//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const ll kMod = 20101009;
const int kMaxn = 1e7 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], sum[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void Getmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
vis[1] = true, mu[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = - mu[i];
}
}
sum[1] = 1ll;
for (int i = 1; i <= lim_; ++ i) {
sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
}
}
int g(int n_, int m_) {
return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
int f(int n_, int m_) {
int lim = std :: min(n_, m_), ret = 0;
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n_ / (n_ / l), m_ / (m_ / l));
ret = (ret + 1ll * (sum[r] - sum[l - 1] + kMod) * g(n_ / l, m_ / l) % kMod) % kMod;
}
return ret;
}
int Sum(ll l, ll r) {
return (1ll * (r - l + 1ll) * (l + r) / 2ll) % kMod;
}
//=============================================================
int main() {
int n = read(), m = read();
int lim = std :: min(n, m), ans = 0;
Euler(lim);
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n / (n / l), m / (m / l));
ans = (ans + 1ll * Sum(l, r) * f(n / l, m / l) % kMod) % kMod;
}
printf("%d", ans);
return 0;
}
/*
7718820 8445880
*/
SP5971 LCMSUM - LCM Sum
(T) 次询问,每次询问给定 (n),求:
[sum_{i=1}^{n}operatorname{lcm}(i,n) ](1<Tle 3 imes 10^5),(1le nle 10^6)。
法一:无脑暴力
先拆 (operatorname{lcm}),原式等价于:
套路的枚举 (gcd(i,n)),调换枚举顺序,原式等价于:
把 (i,n) 中的 (d) 提出来,变为枚举 (frac{i}{d}),消去整除的条件,原式等价于:
代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}),原式等价于:
值得注意的是 (t) 的上界为 (frac{n}{d}),(dtle n)。
调换枚举顺序,先枚举 (t),原式等价于:
套路地消去整除的条件,把 (i) 中的 (t) 提出来,原式等价于:
对于最后的一个求和项,设 (g(x) = sumlimits_{i=1}^{x}i = frac{x(x+1)}{2}),显然可以 (O(1)) 求解,原式等价于:
考虑枚举 (T = dt),显然 (Tle n)。
(mu(t)t) 与 (d) 无关,可以直接考虑枚举 (t|T),原式等价于:
前半块是一个数论分块的形式,可以 (O(sqrt{n})) 求解。
考虑后半块,设 (f(n)=sumlimits_{d|n}mu(d)d),发现它是一个积性函数,可以线性筛筛出,具体地:
其中 (p) 为 (n) 的最小质因子。
此时已经可以线性筛 + 数论分块求解,复杂度 (O(n+Tsqrt{n})),比较菜鸡,时限 500ms 过不了。
考虑筛出 (f) 后再用埃氏筛预处理 (sumlimits_{T=1}^{n}g(leftlfloordfrac{n}{T}
ight
floor)f(T)),输出时乘上 (n),复杂度变为 (O(nlog^2 n + n))。
法二:
同样先拆 (operatorname{lcm}),枚举 (gcd(i,n)),调换枚举顺序,原式等价于:
把 (i,n) 中的 (d) 提出来,变为枚举 (frac{i}{d}),消去整除的条件,原式等价于:
调整枚举对象,上式等价于:
考虑 (sumlimits_{i=1}^{d}[gcd(i,d) = 1]i) 的实际意义,表示 ([1,d]) 中与 (d) 互质的数的和。
当 (d>1) 时,与 (d) 互质的数总是成对存在,即若 (gcd(i,d)=1) 成立,则 (gcd(d-i,d)=1) 成立。
每对这样的数的和为 (d),共有 (frac{varphi(d)}{2}) 对这样的数。
则原式等价于:
可以直接预处理答案。
预处理时先线性筛出 (varphi),再埃氏筛枚举 (i) 的倍数,令它们的答案加上 (frac{varphi(i)i}{2}),最后输出时乘上 (n)。
复杂度 (O(nlog^2 n + T))。
法二代码
//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 1e6;
//=============================================================
ll phi[kMaxn + 10], ans[kMaxn + 10];
int pnum, p[kMaxn + 10];
bool flag[kMaxn + 10];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void GetMax(int &fir, int sec) {
if (sec > fir) fir = sec;
}
void GetMin(int &fir, int sec) {
if (sec < fir) fir = sec;
}
void GetPrime() {
phi[1] = 1, flag[1] = true; //注意初始化
for (int i = 2; i <= kMaxn; ++ i) {
if (! flag[i]) {
p[++ pnum] = i;
phi[i] = i - 1ll;
}
for (int j = 1; j <= pnum && i * p[j] <= kMaxn; ++ j) {
flag[i * p[j]] = true;
if (i % p[j]) {
phi[i * p[j]] = phi[i] * phi[p[j]];
} else {
phi[i * p[j]] = phi[i] * p[j];
break;
}
}
}
for (int i = 1; i <= kMaxn; ++ i) {
for (int j = 1; i * j <= kMaxn; ++ j) {
ans[i * j] += (i == 1 ? 1 : 1ll * phi[i] * i / 2);
}
}
}
//=============================================================
int main() {
GetPrime();
int T = read();
while (T --) {
int n = read();
printf("%lld
", 1ll * ans[n] * n);
}
return 0;
}
[SDOI2015]约数个数和
(T) 次询问,每次询问给定 (n,m)。
定义 >(operatorname{d}(i)) 为 (i) 的约数个数,求:[sum_{i=1}^{n}sum_{j=1}^moperatorname{d}(ij) ](1<T,nle 5 imes 10^4)。
一个结论:
证明
先考虑 (i = p^a),(j=p^b(pin mathrm{primes})) 的情况,有:
对于等式左侧,(p^{a+b}) 的约数个数为 (a+b+1)。
对于等式右侧,在保证 (gcd(x,y)=1) 成立的情况下,有贡献的数对 ((x,y)) 只能是下列三种形式:
- (x>0,y-0),(x) 有 (a) 种取值方法。
- (x=0,y>0),(y) 有 (b) 种取值方法。
- (x=0,y=0)。
则等式右侧贡献次数为 (a+b+1) 次,等于 (p^{a+b}) 的约数个数。
则当 (i = p^a),(j=p^b(pin mathrm{primes})) 时等式成立。
又不同质因子间相互独立,上述情况可拓展到一般的情况。
对 (operatorname{d}(i,j)) 进行化简,代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}),原式等价于:
调换枚举顺序,先枚举 (d),原式等价于:
把各项中的 (d) 提出来,消去整除的条件,原式等价于:
将 (operatorname{d}(ij)) 代回原式,原式等价于:
调换枚举顺序,先枚举 (d),原式等价于:
把 (i,j) 中的 (d) 提出来,变为枚举 (frac{i}{d}, frac{j}{d}),消去整除的条件,原式等价于:
考虑预处理 (S(x) = sumlimits_{i=1}^{x}operatorname{d}(i)),则原式等价于:
线性筛预处理 (mu,operatorname{d}),数论分块求解即可,复杂度 (O(n+Tsqrt{n}))。
代码
//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 5e4 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], num[kMaxn], d[kMaxn]; //num 为最小质因子的次数
ll summu[kMaxn], sumd[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void Getmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
vis[1] = true;
mu[1] = d[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
num[i] = 1;
d[i] = 2;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
num[i * p[j]] = num[i] + 1;
d[i * p[j]] = 1ll * d[i] / num[i * p[j]] * (num[i * p[j]] + 1ll);
break;
}
mu[i * p[j]] = - mu[i];
num[i * p[j]] = 1;
d[i * p[j]] = 2ll * d[i]; //
}
}
for (int i = 1; i <= lim_; ++ i) {
summu[i] = summu[i - 1] + mu[i];
sumd[i] = sumd[i - 1] + d[i];
}
}
//=============================================================
int main() {
Euler(kMaxn - 10);
int T = read();
while (T --) {
int n = read(), m = read(), lim = std :: min(n, m);
ll ans = 0ll;
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n / (n / l), m / (m / l));
ans += 1ll * (summu[r] - summu[l - 1]) * sumd[n / l] * sumd[m / l]; //
}
printf("%lld
", ans);
}
return 0;
}
/*
in
1
32 43
*/
/*
out
15420
*/
P3768 简单的数学题
给定参数 (n)、(p),求:
[left(sum_{i=1}^nsum_{j=1}^nicdot jcdot gcd(i,j) ight)mod p ](nleq10^{10}),(5 imes10^8leq pleq1.1 imes10^9),(pin mathrm{primes})。
时限 4s。
无脑套路暴力。
考虑先枚举 (gcd(i,j)),原式等价于:
提出 (d),变为枚举 (frac{i}{d}) 与 (frac{j}{d}),消去整除的条件,原式等价于:
代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}),原式等价于:
值得注意的是 (t) 的上界为 (frac{n}{d}),(dtle n)。
调换枚举顺序,先枚举 (t),原式等价于:
和上面一样,提出 (t),套路地消去整除的条件,原式等价于:
发现后面两个求和是等差数列乘等差数列的形式。
设 (g(p,q) = sumlimits_{i=1}^{p} sumlimits_{j=1}^{q}ij),(g_{p,q}) 可以通过下式 (O(1)) 计算:
代入原式,原式等价于:
考虑枚举 (T = dt),显然 (Tle n)。
再考虑枚举 (d|T),即可得到 (t = frac{T}{d}),原式等价于:
对于后面这一坨,用 (sumlimits_{d|T}dcdot mu{left(frac{T}{d} ight)} = operatorname{Id} ast mu(T)= varphi(T)) 反演,则原式等价于:
后半块可以数论分块,考虑前半块。
发现前半段即为 (operatorname{Id}^2(T)varphi(T)),又是前缀和形式,考虑杜教筛。
有:
考虑找到一个函数 (g),构造函数 (h = fast g) 使其便于求值,有:
看到同时存在 (d) 和 (frac{n}{d}),考虑把 (d^2) 消去。
令 (g = operatorname{Id}^2),有:
又 (varphi ast 1 = operatorname{Id}),则有:
找到了合适的 (g),套杜教筛的公式。
前一项是自然数的立方和,有 (sumlimits_{i=1}^{n} i^3 = (frac{n(n+1)}{2})^2)。证明详见:自然数前n项平方和、立方和公式及证明 - 百度文库。
后一项直接等差数列求和 + 数论分块求解即可。
「SDOI2017」数字表格
设 (f_{i}) 表示斐波那契数列的第 (i) 项。
(T) 组数据,每次给定参数 (n,m),求:[prod_{i=1}^{n}prod_{j=1}^{m}f_{gcd(i,j)} pmod {10^9 + 7} ](1le Tle 10^3),(1le n,mle 10^6)。
5S,256MB。
以下钦定 (nge m)。
大力化式子,先套路地枚举 (gcd(i,j)),用初中知识把两个 (prod) 化到指数位置,原式等于:
分母套路一波,有:
代回原式,原式等于:
考虑再暴力拆一波,原式等于:
做不动了,但发现变量仅有 (k,d,kd),考虑更换枚举对象改为枚举 (t = kd) 与 (d),则原式等于:
枚举对象变成了约数形式。从后面的式子推前面的式子是比较显然的,可以认为这种枚举 (t=kd) 的形式是一种逆向操作。
设:
(g(t)) 可以用类似埃氏筛的方法 (O(nlog ^2 n)) 地预处理出来。再把 (g) 代回原式,原式等于:
可以考虑预处理 (g(t)) 的前缀积,数论分块枚举指数求解即可。
总时间复杂度 (O(nlog ^2 n + Tsqrt n)),轻微卡常可以过。
//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define LL long long
const LL mod = 1e9 + 7;
const int kN = 1e6;
//=============================================================
LL n, m, ans;
int p_num, p[kN + 10];
bool vis[kN + 10];
LL mu[kN + 10], f[kN + 10], g[kN + 10], prod[kN + 10];
LL invf[kN + 10], invp[kN];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) {
w = (w << 3) + (w << 1) + (ch ^ '0');
}
return f * w;
}
void Chkmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
x_ %= mod;
LL ret = 1;
for (; y_; y_ >>= 1ll) {
if (y_ & 1) ret = ret * x_ % mod;
x_ = x_ * x_ % mod;
}
return ret;
}
void Euler() {
vis[1] = true, mu[1] = 1;
for (int i = 2; i <= kN; ++ i) {
if (! vis[i]) {
p[++ p_num] = i;
mu[i] = -1;
}
for (int j = 1; j <= p_num && i * p[j] <= kN; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
}
void Prepare() {
g[1] = g[2] = 1;
f[1] = f[2] = 1;
invf[1] = invf[2] = 1;
for (int i = 3; i <= kN; ++ i) {
g[i] = 1;
f[i] = (f[i - 1] + f[i - 2]) % mod;
invf[i] = QPow(f[i], mod - 2);
}
Euler();
for (int d = 1; d <= kN; ++ d) {
for (int j = 1; d * j <= kN; ++ j) {
if (mu[j] == 1) {
g[d * j] = g[d * j] * f[d] % mod;
} else if (mu[j] == -1) {
g[d * j] = g[d * j] * invf[d] % mod;
}
}
}
invp[0] = prod[0] = 1;
for (int i = 1; i <= kN; ++ i) {
prod[i] = prod[i - 1] * g[i] % mod;
invp[i] = QPow(prod[i], mod - 2);
}
}
//=============================================================
int main() {
Prepare();
int T = read();
while (T -- ){
n = read(), m = read(), ans = 1;
if (n < m) std::swap(n, m);
for (LL l = 1, r = 1; l <= m; l = r + 1) {
r = std::min(n / (n / l), m / (m / l));
ans = (ans * QPow(prod[r] * invp[l - 1] % mod, 1ll * (n / l) * (m / l))) % mod;
}
printf("%lld
", ans);
}
return 0;
}
P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题
给定参数 (p),有 (T) 组数据,每次给定参数 (A,B,C),求:
[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}left(dfrac{operatorname{lcm}(i,j)}{gcd(i,k)} ight)^{f(type)} ]其中 (f(type)) 的取值如下:
[f(type) = egin{cases} 1 &type = 0\ i imes j imes k &type = 1\ gcd(i,j,k) &type = 2 end{cases}](1le A,B,Cle 10^5),(10^7le ple 1.05 imes 10^9),(pin mathbb{P}),(T=70)。
2.5S,128MB。
先化下原式,原式等于:
发现每一项仅与两个变量有关,设:
发现 (prod) 可以随意交换,则原式等价于:
考虑在 (type) 取值不同时,如何快速求得 (f_1) 与 (f_2)。
一共有 (6) 个需要推导的式子,不妨就叫它们 (1sim 6) 面了(
type = 0
对于 1 面,显然有:
预处理阶乘 + 快速幂即可,单次计算时间复杂度 (O(log n))。
再考虑 2 面,套路地枚举 (gcd),显然有:
指数是个套路,可以看这里:P3455 [POI2007]ZAP-Queries。于是有:
代回原式,略做处理,则原式等于:
像「SDOI2017」数字表格 一样,考虑枚举 (t=kd) 和 (d),则原式等于:
设:
线性筛预处理 (mu) 后,(g_0(t)) 可以用埃氏筛预处理,时间复杂度 (O(nlog n))。再代回原式,原式等于:
预处理 (g_0(t)) 的前缀积和前缀积的逆元,复杂度 (O(nlog n))。
数论分块 + 快速幂计算即可,单次时间复杂度 (O(sqrt nlog n))。
type = 1
考虑 3 面,把 (prod k) 扔到指数位置,有:
再把 (prod j) 也扔到指数位置,引入 (operatorname{sum}(n) = sum_{i=1}^{n} i = frac{n(n+1)}{2}),原式等于:
预处理 (i^i) 的前缀积,复杂度 (O(nlog n))。
指数可以 (O(1)) 算出,然后快速幂,单次时间复杂度 (O(log n))。
根据费马小定理,指数需要对 (p - 1) 取模。注意 (p-1) 不是质数,计算 (operatorname{sum}) 时不能用逆元,但乘不爆 LL
,直接算就行。
再考虑 4 面,发现 (k) 与 (gcd) 无关,则同样把 (prod k) 扔到指数位置,则有:
套路地枚举 (gcd),原式等于:
大力化简指数,有:
指数化不动了,代回原式,原式等于:
同 2 面的情况,先展开一下,再枚举 (t=kd) 和 (d),原式等于:
与二面相同,设:
(g_1(t)) 可以用埃氏筛套快速幂筛出,时间复杂度 (O(nlog^2 n))。再代回原式,原式等于:
同样预处理 (g_1(t)) 的前缀积及其逆元,时间复杂度 (O(nlog n))。
整除分块 + 快速幂即可,单次时间复杂度 (O(sqrt nlog n))。
注意指数的取模。
type = 2
考虑 5 面,手段同上,大力反演化简一波,再调换枚举对象,则有:
引入 (operatorname{fac}(n) = prod_{i=1}^{n} i),再根据枚举对象调整一下指数,原式等于:
指数中出现了一个经典的狄利克雷卷积的形式,对其进行反演。
将 ((operatorname{Id}ast mu) (n)= varphi (n)) 代入原式,原式等于:
预处理 (t^{varphi(t)}) 的前缀积及逆元,阶乘的前缀积及阶乘逆元,(pmod {p-1}) 下的 (varphi(t)) 的前缀和(指数
),时间复杂度 (O(nlog n))。
同样整除分块 + 快速幂即可,单次时间复杂度 (O(sqrt nlog n))。
然后是最掉 sans 的 6 面。有 (gcd(i,j,k) = gcd(gcd(i,j), k)),考虑先枚举 (gcd(i,j)),然后套路化式子,则有:
先考虑最外面的指数,这也是个套路,可以参考 一个例子。用 (operatorname{Id} = varphi ast 1) 反演,显然有:
再考虑里面的指数,发现这式子在 2 面已经推了一遍了,于是直接拿过来用,有:
将化简后的两个指数代入原式,原式等于:
与 2、4 面同样套路地,考虑枚举 (t=yd) 和 (d),再略作调整,原式等于:
发现要同时枚举 (d) 和 (x),化不动了。
从题解里学到一个比较神的技巧,考虑把 (d) 拆成 (x) 和 (frac{d}{x}) 分别计算贡献再相乘,即分别计算下两式:
先考虑 (x) 的情况,首先把枚举 (x) 调整到最外层。设 (operatorname{lim}=max(a,b,c)),则原式等于:
把 (prod {t}) 挪到指数位置,原式等于:
指数中又出现了一个经典的狄利克雷卷积的形式,对其进行反演。
将 ((mu ast 1) (n)= epsilon (n)=[n=1]) 代入原式,原式等于:
得到了一个非常优美的式子,而且发现这个式子是 5 面最终答案的一部分。同 5 面的做法,直接整除分块即可。
再考虑 (frac{d}{x}) 的情况,同上先把枚举 (x) 放到最外层,并调整一下指数,则原式等于:
考虑枚举 (dx),替换原来的 (d),注意一下这里的倍数关系。原式等于:
发现最内层的式子 (prod_{d|t}d^{muleft(frac{t}{d} ight)}),即为二面处理过的 (g_0(t))。直接代入,原式等于:
一个小结论,证明可以看 这里:
则原式等于:
于是可以先对外层整除分块,再对内层整除分块。
然后就做完了,哈哈哈。
一些实现上的小技巧:
- 逆元能预处理就预处理。
- 注意对指数取模,模数为 (p-1)。
//知识点:莫比乌斯反演
/*
By:Luckyblock
用了比较清晰易懂的变量名,阅读体验应该会比较好。
vsc 的自动补全真是太好啦!
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
using std::min;
using std::max;
#define LL long long
const int Lim = 1e5;
const int kN = 1e5 + 10;
//=============================================================
LL A, B, C, mod, ans;
int T, p_num, p[kN];
bool vis[kN];
LL mu[kN], phi[kN], fac[kN], g[2][kN];
LL sumphi[kN], prodt_phi[kN], prodi_i[kN], prodg[2][kN];
LL inv[kN], inv_fac[kN], inv_prodt_phi[kN], inv_prodg[2][kN];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) {
w = (w << 3) + (w << 1) + (ch ^ '0');
}
return f * w;
}
void Chkmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
x_ %= mod;
y_ %= mod - 1;
LL ret = 1;
for (; y_; y_ >>= 1ll) {
if (y_ & 1) ret = ret * x_ % mod;
x_ = x_ * x_ % mod;
}
return ret;
}
LL Inv(LL x_) {
return QPow(x_, mod - 2);
}
LL Sum(LL n_) {
return ((n_ * (n_ + 1ll)) / 2ll) % (mod - 1);
}
void Euler() {
vis[1] = true, mu[1] = phi[1] = 1; //初值
for (int i = 2; i <= Lim; ++ i) {
if (! vis[i]) {
p[++ p_num] = i;
mu[i] = -1;
phi[i] = i - 1;
}
for (int j = 1; j <= p_num && i * p[j] <= Lim; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
phi[i * p[j]] = phi[i] * p[j];
break;
}
mu[i * p[j]] = -mu[i];
phi[i * p[j]] = phi[i] * (p[j] - 1);
}
}
}
void Prepare() {
Euler();
inv[1] = fac[0] = prodt_phi[0] = prodi_i[0] = 1;
for (int i = 1; i <= Lim; ++ i) {
g[0][i] = g[1][i] = 1;
fac[i] = 1ll * fac[i - 1] * i % mod;
sumphi[i] = (sumphi[i - 1] + phi[i]) % (mod - 1);
prodi_i[i] = prodi_i[i - 1] * QPow(i, i) % mod;
if (i > 1) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
prodt_phi[i] = prodt_phi[i - 1] * QPow(i, phi[i]) % mod;
inv_prodt_phi[i] = Inv(prodt_phi[i]);
}
for (int d = 1; d <= Lim; ++ d) {
for (int j = 1; d * j <= Lim; ++ j) {
int t = d * j;
if (mu[j] == 1) {
g[0][t] = g[0][t] * d % mod;
g[1][t] = g[1][t] * QPow(1ll * d, 1ll * t * t) % mod;
} else if (mu[j] == -1) {
g[0][t] = g[0][t] * inv[d] % mod;
g[1][t] = g[1][t] * Inv(QPow(1ll * d, 1ll * t * t)) % mod;
}
}
}
inv_prodg[0][0] = prodg[0][0] = 1;
inv_prodg[1][0] = prodg[1][0] = 1;
inv_prodt_phi[0] = 1;
for (int i = 1; i <= Lim; ++ i) {
for (int j = 0; j <= 1; ++ j) {
prodg[j][i] = prodg[j][i - 1] * g[j][i] % mod;
inv_prodg[j][i] = Inv(prodg[j][i]);
}
}
}
LL f1(LL a_, LL b_, LL c_, int type) {
if (! type) return QPow(fac[a_], b_ * c_);
if (type == 1) return QPow(prodi_i[a_], Sum(b_) * Sum(c_));
LL ret = 1, lim = min(min(a_, b_), c_);
for (LL l = 1, r = 1; l <= lim; l = r + 1) {
r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
ret = ret * QPow(fac[a_ / l], (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
}
return ret;
}
LL f2_2(LL a_, LL b_) {
LL ret = 1;
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
ret = ret * QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l)) % mod;
}
return ret;
}
LL f2(LL a_, LL b_, LL c_, int type) {
LL ret = 1;
if (! type) {
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
LL val = QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l));
ret = (ret * QPow(val, c_)) % mod;
}
} else if (type == 1) {
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
LL val = QPow(prodg[1][r] * inv_prodg[1][l - 1], Sum(a_ / l) * Sum(b_ / l));
ret = (ret * QPow(val, Sum(c_))) % mod;
}
} else {
LL lim = min(min(a_, b_), c_);
for (LL l = 1, r = 1; l <= lim; l = r + 1) {
r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
ret = ret * QPow(f2_2(a_ / l, b_ / l), (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (c_ / l)) % mod;
ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
}
}
return ret;
}
//=============================================================
int main() {
T = read(), mod = read();
Prepare();
while (T -- ) {
A = read(), B = read(), C = read();
for (int i = 0; i <= 2; ++ i) {
ans = f1(A, B, C, i) * f1(B, A, C, i) % mod;
ans = ans * Inv(f2(A, B, C, i)) % mod * Inv(f2(A, C, B, i)) % mod;
printf("%lld ", ans);
}
printf("
");
}
return 0;
}
写在最后
参考资料:
Oi-Wiki-莫比乌斯反演
算法学习笔记(35): 狄利克雷卷积 By: Pecco
题解 SP5971 【LCMSUM - LCM Sum】 - BJpers2 的博客
题解 SP5971 【LCMSUM - LCM Sum】 - Venus 的博客
题解 P3327 【[SDOI2015]约数个数和】 - suncongbo 的博客
把 Oi-Wiki 上的内容进行了复制 整理扩充。
我是个没有感情的复读机(大雾)