嘛......开个坑吧 咕不咕就不知道了qwq......
数论分块
这个东西之前水过一篇来着：数论分块
狄利克雷卷积
定义
对于两个数论函数(f(n),g(n))，定义他们的狄利克雷卷积（Dirichlet卷积）为
注意柿子中(d,frac{n}{d})是对称的，所以
也就是说Dirichlet卷积满足交换律，同时狄利克雷卷积也满足结合律，证明的话暴力把柿子写开就好了，这里就不写了qaq
积性函数
对于一个数论函数(f(n))，如果对于两个互质的整数(a,b)，有(f(ab)=f(a)f(b))，则称(f(n))为积性函数。
如果对于所有的(a,b)都有(f(ab)=f(a)f(b))，那么(f(n))就是一个完全积性函数。
特别的对于两个积性函数(f(n),g(n))，他们的Dirichlet卷积也是积性函数。
常见的积性函数
(varepsilon(n)=[n=1])，这个函数很重要，它也是Dirichlet卷积的单位元（所有函数卷它都是原来的函数，类比乘法中的(1)）。
(1(n)equiv 1)，常数函数，它永远等于(1)。
(id_k(n)=n^k)，恒等函数，一般直接将(id_1(n))记作(id(n))。
(sigma_k(n)=sumlimits_{dmid n}d^k)，除数函数，(sigma_0(n))就是除数个数，一般写作( au(n))或(d(n))，(sigma_1(n))即约数之和，记作(sigma(n))。
(varphi(n))，Euler函数，即(1cdots n-1)内与(n)互质的数的个数（同时也只个积性函数），同时也是个积性函数。
(mu(n))，即莫比乌斯函数，可以理解为Dirichlet卷积中(1)的逆元（类比乘法逆元），也就是说((mu*1)(n)=sumlimits_{dmid n} mu(d)=varepsilon(n)=[n=1])。（定义在下面会讲）
以及他们的狄利克雷卷积的形式
莫比乌斯函数
终于讲到点东西了...
上面我们知道(sumlimits_{dmid n} mu(d)=[n=1])，那么(mu)的定义到底是啥呢？
先把(n)质因数分解
那么
鬼知道它为什么长这样......
特别的(mu(1)=1)
如果(n)有一个平方因子的话(mu(n))就是(0)，不然就是(-1)的质因数个数次方。
首先对于(n=1)，有(sumlimits_{dmid n} mu(d)=1)
那么其他情况呢？
令(n'=prod_{i=1}^k p_i)，显然有
考虑(sumlimits_{d mid n'} mu(d))，枚举质因子的个数，有
显然是(0)。
故(sumlimits_{dmid n} mu(d) =[n=1])
我们也得到了反演结论：
线性筛(mu)
(mu)是一个积性函数，所以我们考虑在线性筛的过程中算出(mu)（不知道的话可以去学一学线性筛qwq）。有如下代码
void getmu() {
mu[1]=1;
rep(i,2,n+1) {
if(!vis[i]) {mu[i]=-1;p[pn++]=i;}
for(int j=0;j<pn&&i*p[j]<=n;j++) {
vis[i*p[j]]=1;
if (i%p[j]) mu[i*p[j]]=mu[i]*mu[p[j]]; // 这里也可以写成mu[i*p[j]]=-mu[i]
else {mu[i*p[j]]=0;break;}
}
}
}
莫比乌斯反演
对于两个数论函数(f,g)
有
证明的话可以考虑Dirichlet卷积
两边同时乘(mu)
所以
其实反演的柿子并不会怎么用到2333
例题
还是来看题吧
LuoguP3455 [POI2007]ZAP-Queries
给(n,m,k)，让你求
多组数据。
不难发现这个柿子等价于
把这些两两互质的同时乘上(k)就是(gcd=k)的数。
那么简化问题，令
答案就是(S(leftlfloorfrac{n}{k} ight floor,leftlfloorfrac{m}{k} ight floor))。
发现后面中括号就是一个单位函数，于是
考虑把(d)提到外面
后面两个(Sigma)就是(leftlfloorfrac{n}{d} ight floorleftlfloorfrac{m}{d} ight floor)。
所以(S(n,m)=sumlimits_{d} mu(d)leftlfloorfrac{n}{d} ight floorleftlfloorfrac{m}{d} ight floor)，(d)最大显然是到(min(n,m))
线性筛出莫比乌斯函数预处理出前缀和，数论分块就好了。
复杂度(O(Tsqrt{n}))（(T)为数据组数）
#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i=(a);i<(b);++i)
#define per(i,a,b) for (int i=(a)-1;i>=(b);--i)
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef double db;
typedef long long ll;
typedef pair<int,int> PII;
typedef vector<int> VI;
const int N=5e4+10;
int mu[N],sum[N],vis[N];
int p[N],pn;
void init(int n) {
mu[1]=1;
rep(i,2,n+1) {
if(!vis[i]) {
mu[i]=-1;
p[pn++]=i;
}
for(int j=0;j<pn&&i*p[j]<=n;j++) {
vis[i*p[j]]=1;
if (i%p[j]) mu[i*p[j]]=mu[i]*mu[p[j]];
else {mu[i*p[j]]=0;break;}
}
}
rep(i,1,n+1) sum[i]=sum[i-1]+mu[i];
}
#define S(l,r) (sum[r]-sum[l-1])
int solve(int n,int m,int d) {
n/=d,m/=d; int tn=min(n,m);
ll ans=0;
for(int l=1,r=0;l<=tn;l=r+1) {
r=min(n/(n/l),m/(m/l));
ans+=1ll*S(l,r)*(n/l)*(m/l);
}
return ans;
}
int main() {
init(50000);
int _,n,m,d; for(scanf("%d",&_);_;_--) {
scanf("%d%d%d",&n,&m,&d);
printf("%d
",solve(n,m,d));
}
return 0;
}
LuoguP2522 [HAOI2011]Problem b
这个就是上面那题吧......类似二维前缀和容斥一下就没了，代码就不放了qwq实际上是为了提示双倍经验
LuoguP1829 [国家集训队]Crash的数字表格 / JZPTAB
给(n,m)，让你求
首先(operatorname{lcm}(i,j)=frac{ij}{gcd(i,j)})
那么就是
枚举(gcd)
把(d)提到外面，然后套路的变一下
下面令(S(n,m)=sumlimits_{i=1}^nsumlimits_{j=1}^m [gcd(i,j)=1] ij)
则(Ans=sumlimits_{d} d imes S(leftlfloorfrac{n}{d} ight floor,leftlfloorfrac{m}{d} ight floor))
我们想快一点求出(S)，然后对(d)数论分块。
推柿子
后面的两个(Sigma)就直接等差数列求和吧，然后预处理出(mu(d)d^2)的前缀和，数论分块就好了。
我们能在根号的时间内求出(S)，求(Ans)的话也用数论分块就好了。
这样就是外层数论分块套内层数论分块。复杂度不太清楚(qwq)，好像上界是(O(n))的，总之跑的过就是了。
#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i=(a);i<(b);++i)
#define per(i,a,b) for (int i=(a)-1;i>=(b);--i)
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef double db;
typedef long long ll;
typedef pair<int,int> PII;
typedef vector<int> VI;
const int maxn=1e7,N=maxn+10,P=20101009;
inline int add(int x,int y) {return (x+=y)>=P?x-P:x;}
inline int sub(int x,int y) {return (x-=y)<0?x+P:x;}
inline int normal(int x) {return x<0?x+P:x;}
inline int sqr(int x) {return 1ll*x*x%P;}
inline int fpow(int x,int y) {
int ret=1; for(;y;y>>=1,x=1ll*x*x%P)
if(y&1) ret=1ll*ret*x%P;
return ret;
}
const int i2=(P+1)/2,i6=fpow(6,P-2);
int mu[N],sum[N],vis[N],p[N],pn;
void init(int n) {
mu[1]=1;
rep(i,2,n+1) {
if(!vis[i]) {p[pn++]=i;mu[i]=-1;}
for(int j=0;j<pn&&p[j]*i<=n;j++) {
vis[i*p[j]]=1;
if(i%p[j]==0) {mu[i*p[j]]=0;break;}
else mu[i*p[j]]=mu[i]*mu[p[j]];
}
}
rep(i,1,n+1) mu[i]=normal(mu[i]);
rep(i,1,n+1) sum[i]=add(sum[i-1],1ll*mu[i]*sqr(i)%P);
}
inline int s1(int l,int r) {return 1ll*(r-l+1)*(l+r)%P*i2%P;}
inline int ss(int l,int r) {return sub(sum[r],sum[l-1]);}
int S(int n,int m) {
int tn=min(n,m),ans=0;
for(int l=1,r=0;l<=tn;l=r+1) {
r=min(n/(n/l),m/(m/l));
ans=add(ans,1ll*ss(l,r)*s1(1,n/l)%P*s1(1,m/l)%P);
}
return ans;
}
int solve(int n,int m) {
int ans=0,tn=min(n,m);
for(int l=1,r=0;l<=tn;l=r+1) {
r=min(n/(n/l),m/(m/l));
ans=add(ans,1ll*s1(l,r)*S(n/l,m/l)%P);
}
return ans;
}
int main() {
init(maxn);
int n,m;scanf("%d%d",&n,&m);
printf("%d
",solve(n,m));
return 0;
}
LuoguP3704 [SDOI2017]数字表格
让你求
多组数据
其中
第一次以为是Fibonacci然后调了半天样例怎么都过不去qwq......
枚举柿子中的(gcd)
发现上面的柿子我们在第一题中就推过了2333
令
原来的柿子就是
这个柿子我们已经可以数论分块套数论分块做到(O(n))的复杂度了，但是它还多测qaq，所以我们再推一下柿子
注意到(leftlfloorfrac{leftlfloorfrac{n}{a} ight floor}{b} ight floor=leftlfloorfrac{n}{ab} ight floor)（感性理解一下先每(a)个数中选一个出来再在选出的数中每(b)个中选一个，等价于在原来(n)个数中每(ab)个里面选一个出来）
枚举(i=dk)，有
指数上面的该数论分块还是数论分块，然后发现括号内的可以在调和级数复杂度内预处理出前缀的乘积，要取一段区间的乘积的话乘上逆元就好了。复杂度(O(n ln n +Tsqrt{n}))。
#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i=(a);i<(b);++i)
#define per(i,a,b) for (int i=(a)-1;i>=(b);--i)
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef double db;
typedef long long ll;
typedef pair<int,int> PII;
typedef vector<int> VI;
const int N=1e6+10,maxn=1e6,P=1e9+7;
inline int add(int x,int y) {return (x+=y)>=P?x-P:x;}
inline int sub(int x,int y) {return (x-=y)<0?x+P:x;}
inline int fpow(int x,int y) {
int ret=1; for(;y;y>>=1,x=1ll*x*x%P)
if(y&1) ret=1ll*ret*x%P;
return ret;
}
inline int normal(int x) {return x<0?x+P:x;}
int vis[N],p[N],pn,mu[N];
void getmu(int n) {
mu[1]=1;
rep(i,2,n+1) {
if(!vis[i]) {mu[i]=-1;p[pn++]=i;}
for(int j=0;j<pn&&i*p[j]<=n;j++) {
vis[i*p[j]]=1;
if(i%p[j]==0) {mu[i*p[j]]=0;break;}
else mu[i*p[j]]=-mu[i];
}
}
}
int fib[N],f[N],ifib[N];
void init(int n) {
getmu(n);
fib[0]=0,fib[1]=1;
rep(i,2,n+1) fib[i]=add(fib[i-1],fib[i-2]);
rep(i,0,n+1) ifib[i]=fpow(fib[i],P-2),f[i]=1;
rep(i,1,n+1) if(mu[i]!=0) {
for(int j=i;j<=n;j+=i)
f[j]=1ll*f[j]*(mu[i]==-1?ifib[j/i]:fib[j/i])%P;
}
rep(i,1,n+1) f[i]=1ll*f[i]*f[i-1]%P;
}
int solve(int n,int m) {
int tn=min(n,m),ans=1;
for(int l=1,r=0;l<=tn;l=r+1) {
r=min(n/(n/l),m/(m/l));
ans=1ll*ans*fpow(1ll*f[r]*fpow(f[l-1],P-2)%P,1ll*(n/l)*(m/l)%(P-1))%P;
}
return ans;
}
int main() {
#ifdef LOCAL
freopen("a.in","r",stdin);
#endif
init(maxn);
int n,m,_; for(scanf("%d",&_);_;_--) {
scanf("%d%d",&n,&m);
printf("%d
",solve(n,m));
}
return 0;
}