前置知识：莫比乌斯反演专题(基础篇)
杜教筛
设\(f(n)\)为一积性函数
求\(S(n)=\sum_{i=1}^nf(i)\)，\(n\leq 10^{10}\)
我们考虑给\(f(n)\)卷上另一个积性函数\(f(n)\)再做前缀和
这样就得到了一个子问题
接着我们特别看\(d=1\)这一项，也就是\(g(1)\times S(n)\)
这个式子中，\(g(1)\)很容易求，\(g\)和\(f\)卷积的前缀和可以通过构造合适的\(g\)让他也很容易求，后边的可以通过整除分块+递归解决
如果我们事先预处理除一部分的\(S\)那么可以证明计算\(S(n)\)的复杂度约为\(O(n^{\frac{2}{3}})\)
[BZOJ3944]--Sum
最常见的
\[\sum_{i=1}^{n}\mu(i) \] \[\sum_{i=1}^{n}\varphi(i) \]先看\(mu(n)\)，我们的任务是构造一个合适的\(g\)
我们知道\(\mu \times I = e\)（\(\times\) 表示狄利克雷卷积）
因此我们构造\(g(n)=I(n)=1\)
那么\(g\)与\(\mu\)的前缀和还是\(e\)，因此这一块的和为1
也就是说
可以杜教筛解决
再看\(\varphi(n)\)，我们知道\(\varphi(n)\times I=id\)
那么我们还考虑卷上\(g(n)=I(n)=1\)
那么式子的值为
左边是等差数列求和很容易算
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 4641600;
LL s[N];
LL mu[N];
int v[N],prime[N],tot=0;
LL phi[N];
LL p[N];
void init(int n)
{
phi[1]=1;
mu[1]=1;
for(int i=2;i<=n;i++)
{
if(!v[i])
{
prime[++tot]=i;
v[i]=i;
mu[i]=-1;
phi[i]=i-1;
}
for(int j=1;j<=tot;j++)
{
if(prime[j]>v[i]||i*prime[j]>n) break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
else
{
mu[i*prime[j]]=mu[i]*mu[prime[j]];
phi[i*prime[j]]=1ll*phi[i]*(prime[j]-1);
}
}
}
for(int i=1;i<=n;i++)
mu[i]=mu[i-1]+mu[i];
for(int i=1;i<=n;i++)
phi[i]=phi[i-1]+phi[i];
}
bool vis[N];
LL S(LL n,int k)
{
if(n<=N-10) return mu[n];
if(vis[k]) return s[k];
LL res=1;
LL l=2,r;
for(;l<=n;l=r+1)
{
r=n/(n/l);
res=res-1ll*(r-l+1)*S(n/l,1ll*k*l);
}
vis[k]=1;
s[k]=res;
return res;
}
LL P(LL n,int k)
{
if(n<=N-10) return phi[n];
if(vis[k]) return p[k];
LL res=1ll*(1ll+n)*n/2;
LL l=2,r;
for(;l<=n;l=r+1)
{
r=n/(n/l);
res=res-1ll*(r-l+1)*P(n/l,k*l);
}
vis[k]=1;
p[k]=res;
return res;
}
int main()
{
freopen("sum.in","r",stdin);
freopen("sum.out","w",stdout);
init(N-10);
int T;
cin>>T;
while(T--)
{
LL n;
cin>>n;
cout<<P(n,1)<<' ';
memset(vis,0,sizeof(vis));
cout<<S(n,1)<<endl;
memset(vis,0,sizeof(vis));
}
return 0;
}
[HDOJ5608]function
题目给出了一个函数\(F(n)=n^2-3n+2\)
并且他保证\(F(n)=\sum_{d|n}f(d)\)
求\(\sum_{i=1}^nf(i)\)
我们先求出\(f(n)\)的递推式
莫比乌斯反演一下可得
\(f(n)=\sum_{d|n}F(d)\mu(\frac{n}{d})\)
这显然是一个积性函数
但是\(n\)的范围很大，考虑杜教筛
发现题目给的这个\(F(n)\)正好可以看成\(f(n)\)和\(I(n)\)的卷积
因此答案为
其中左半部分的三项可以分开算，其中
\(1^2+2^2+3^2+……n^2=n(n+1)(2n+1)/6\)
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e6+7;
const LL mod = 1e9+7;
LL Pow(LL a,LL b)
{
LL res=1;
while(b)
{
if(b&1) res=1ll*res*a%mod;
a=1ll*a*a%mod;
b>>=1;
}
return res;
}
LL inv(LL n)
{
return Pow(n,mod-2);
}
LL inv6=Pow(3,mod-2),inv2=Pow(2,mod-2);
LL F(LL n){return (1ll*n*(n+1)/2%mod*(2ll*n+1)%mod*inv6%mod-3*(1+n)*n/2%mod+2*n+mod)%mod;}
int mu[N],prime[N],tot=0;
int v[N];
LL f[N];
void init(int n)
{
mu[1]=1;
f[1]=0;
for(int i=2;i<=n;i++)
{
if(!v[i])
{
v[i]=i;
mu[i]=-1;
prime[++tot]=i;
}
for(int j=1;j<=tot;j++)
{
if(prime[j]>v[i]||i*prime[j]>n) break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
else mu[i*prime[j]]=mu[i]*mu[prime[j]];
}
}
for(LL d=1;d<=n;d++)
for(LL T=d;T<=n;T+=d)
f[T]=(f[T]+(d*d-3*d+2)%mod*mu[T/d]+mod)%mod;
for(int i=1;i<=n;i++)
f[i]=(f[i-1]+f[i])%mod;
}
LL bin[N],top=0;
int tag[N];
int cnt=0;
unordered_map<LL,LL>s,vis;
LL S(LL n)
{
if(n<=1e6) return f[n];
if(vis[n]) return s[n];
LL res=F(n);
LL l=2,r;
for(;l<=n;l=r+1)
{
r=n/(n/l);
res=res-(1ll*(r-l+1)*S(n/l)%mod);
res=(res%mod+mod)%mod;
}
vis[n]=1;
s[n]=res;
return res;
}
inline int read()
{
int X=0; bool flag=1; char ch=getchar();
while(ch<'0'||ch>'9') {if(ch=='-') flag=0; ch=getchar();}
while(ch>='0'&&ch<='9') {X=(X<<1)+(X<<3)+ch-'0'; ch=getchar();}
if(flag) return X;
return ~(X-1);
}
void write(LL x)
{
if(x < 0) {
putchar('-');
x = -x;
}
if(x > 9)
write(x/10);
putchar(x % 10 + '0');
return;
}
int main()
{
freopen("a.in","r",stdin);
freopen("a.out","w",stdout);
init(1e6);
LL T;
cin>>T;
LL n;
while(T--)
{
n=read();
if(n<=1e6)
{
write(f[n]);
printf("\n");
continue;
}
write(S(n));
printf("\n");
}
return 0;
}
【51nod 2026】Gcd and Lcm
题目要求
\[\sum_{i=1}^n\sum_{j=1}^nf(gcd(i,j))f(lcm(i,j)) \]其中\(f(n)=\sum_{d|n}\mu(d)d\)
一开始想把这个式子莫反一下，但是发现做不下去了……
考虑分析一席\(f(n)\)的性质
这个函数显然是积性函数，因此我们对每个质因子分别考虑
\(f(p^k)=1-p\)，也就是说对于每个质因子而言，他的指数是多少根本没有关系，因此一个数的\(f\)只跟他质因子的种类数有关
接着我们分析\(gcd(i,j)\)和\(lcm(i,j)\)
1：对于\(i,j\)共有的因子，他们被乘了两次,设他们是A
2：对于\(i,j\)各自有的因子，他们都只乘了一次，假设分别是B,C
也就是 A(BCA)
我们交换一下顺序可以得到(AB)(AC)
也就是\(f(gcd(i,j))\times f(lcm(i,j))=f(i)\times f(j)\)
式子变为
问题变成了求\(f\)的前缀和
也就是
后半部分可以整除分块
问题进一步转化为求\(\mu(d)d\)的前缀和
考虑杜教筛
我们考虑构造一个\(g(n)=Id(n)=n\)
也就是那么$$g\times (\mu(d)d)$$
因此前缀和就是
\[S(n)=\frac{1-\sum_{d=2}^niS(\lfloor \frac{n}{d}\rfloor)}{1} \]带进去算就行了
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int mod = 1e9+7;
const int N =1e6+7;
int mu[N];
LL F[N];
int v[N],prime[N],tot=0;
void init(int n)
{
mu[1]=1;
for(int i=2;i<=n;i++)
{
if(!v[i])
{
v[i]=i;
prime[++tot]=i;
mu[i]=-1;
}
for(int j=1;j<=tot;j++)
{
if(v[i]<prime[j]||i*prime[j]>n) break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
else mu[i*prime[j]]=mu[i]*mu[prime[j]];
}
}
for(LL i=1;i<=n;i++)
F[i]=(F[i-1]+mu[i]*i%mod)%mod;
}
unordered_map<LL,LL> s,vis;
LL S(LL n)
{
if(n<=1e6) return F[n];
if(vis[n]) return s[n];
LL res=1;
LL l=2,r;
for(;l<=n;l=r+1)
{
r=n/(n/l);
res=(res-(l+r)*(r-l+1)/2%mod*S(n/l)%mod+mod)%mod;
}
vis[n]=1;
s[n]=res;
return res;
}
LL f(LL n)
{
LL l=1,r;
LL ans=0;
for(;l<=n;l=r+1)
{
r=(n/(n/l));
ans=(ans+(S(r)-S(l-1)+mod)%mod*(n/l)%mod)%mod;
}
return ans;
}
int main()
{
freopen("c.in","r",stdin);
freopen("c.out","w",stdout);
init(1e6);
LL n;
cin>>n;
cout<<f(n)*f(n)%mod;
return 0;
}
【51Nod 1238最小公倍数之和 V3
毒瘤
求
莫反的部分详见莫比乌斯反演专题(基础篇)中的【BZOJ2693】jzptab
我么这里直接用最后的的式子
其中\(S(n)=\sum_{i=1}^ni\)
所以关键的问题还是筛出来\(f(n)=n\sum_{d|n}\mu(d)d\)的前缀和
前边的部分显然我们卷上一个\(g(n)=n^2\)就可以杜教筛了，整个式子整除分块一下就行了
但是，别忘了外边还有一个整除分块，这样就成\(O(n)\)的了
因此我们只能考虑别的方法
通过观察打表 等方式我们发现\(f\times (n^2) = (n)\)
至于具体证明抱歉我不懂，因此我们卷上一个\(g(n)=n^2\)
前缀和变为
就可以解决了
#include<bits/stdc++.h>
using namespace std;
const int N = 4641600;
typedef unsigned long long LL;
int mu[N],prime[N],tot=0;
int v[N];
LL F[N];
const LL mod = 1e9+7;
void init(LL n)
{
F[1]=1;
for(int i=2;i<=n;i++)
{
if(!v[i])
{
v[i]=i;
F[i]=1LL*i*((1-i+mod)%mod)%mod%mod;
prime[++tot]=i;
}
for(int j=1;j<=tot;j++)
{
if(prime[j]>v[i]||i*prime[j]>n) break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
{
F[i*prime[j]]=1LL*F[i]*prime[j]%mod;
break;
}
else F[i*prime[j]]=1LL*F[i]*F[prime[j]]%mod;
}
}
for(int i=1;i<=n;i++)
F[i]=(F[i-1]+F[i])%mod;
}
unordered_map<LL,LL>s,vis;
LL Pow(LL a,LL b)
{
LL res=1;
while(b)
{
if(b&1) res=1LL*res*a%mod;
a=1LL*a*a%mod;
b>>=1LL;
}
return res;
}
LL inv6=Pow(6,mod-2),inv2=Pow(2,mod-2);
LL sqr(LL n){n%=mod;return 1LL*n*(n+1)%mod*(2LL*n+1)%mod*inv6%mod;}
LL sum(LL n){n%=mod;return 1LL*(1+n)%mod*n%mod*inv2%mod;}
bool flag=1;
LL Sum(LL n)
{
if(n<=N-10) return F[n];
if(vis[n]) return s[n];
LL res=sum(n);
LL l=2,r;
for(;l<=n;l=r+1)
{
r=n/(n/l);
res=(res-1LL*((sqr(r)-(sqr(l-1)))+mod)%mod*Sum(n/l)%mod+mod)%mod;
}
vis[n]=1;
s[n]=res;
return res;
}
int main()
{
freopen("d.in","r",stdin);
freopen("d.out","w",stdout);
init(N-10);
LL n;
cin>>n;
LL ans=0;
LL l=1,r;
for(;l<=n;l=r+1)
{
r=(n/(n/l));
ans=(ans+1LL*sum(n/l)*sum(n/l)%mod*((Sum(r)-Sum(l-1)+mod)%mod)%mod)%mod;
}
cout<<ans;
return 0;
}
51【Nod 1237】最大公约数之和 V3
同上莫反部分详见莫比乌斯反演专题(基础篇)中的 [NOI2010day1]T1-能量采集
最后的式子为
直接杜教筛处理\(\varphi(n)\)的前缀和就行了
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod = 1e9+7;
const LL N = 4641600;
LL phi[N];
LL v[N],prime[N],tot=0;
void init(LL n)
{
phi[1]=1;
for(LL i=2;i<=n;i++)
{
if(!v[i])
{
v[i]=i;
phi[i]=i-1;
prime[++tot]=i;
}
for(LL j=1;j<=tot;j++)
{
if(prime[j]>v[i]||i*prime[j]>n) break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
{
phi[i*prime[j]]=1LL*phi[i]*prime[j]%mod;
break;
}
else phi[i*prime[j]]=1LL*phi[i]*(prime[j]-1)%mod;
}
}
for(LL i=1;i<=n;i++)
phi[i]=(phi[i-1]+phi[i])%mod;
}
LL Pow(LL a,LL b)
{
LL res=1;
while(b)
{
if(b&1) res=1LL*res*a%mod;
a=1LL*a*a%mod;
b>>=1;
}
return res;
}
LL inv2=Pow(2ll,mod-2);
unordered_map<LL,LL>vis,s;
LL mul(LL a,LL b)
{
LL res=0;
while(b)
{
if(b&1) res=(res+a)%mod;
a=(a+a)%mod;
b>>=1;
}
return res;
}
LL Sum(LL n)
{
if(n<=N-10) return phi[n];
if(vis[n]) return s[n];
LL res=mul(mul(n,n+1),inv2);
LL l=2,r;
for(;l<=n;l=r+1)
{
r=(n/(n/l));
res=(res-1LL*(r-l+1)%mod*Sum(n/l)%mod)%mod;
}
res=(res+mod)%mod;
vis[n]=1;
s[n]=res;
return res;
}
int main()
{
freopen("d.in","r",stdin);
freopen("d.out","w",stdout);
init(N-10);
LL n;
cin>>n;
LL res=0;
LL l=1,r;
for(;l<=n;l=r+1)
{
r=n/(n/l);
res=(res+mul(Sum(r)-Sum(l-1),mul(n/l,n/l))+mod)%mod;
}
cout<<res;
return 0;
}
[bzoj 4176] Lucas的数论
同上莫反部分详见莫比乌斯反演专题(基础篇)中的 [SDOI2015] 约数个数和
\[ans=\sum_{d=1}^n\mu(d)\sum_{x=1}^{\lfloor \frac{n}{d} \rfloor}\lfloor \frac{n}{dx} \rfloor \sum_{y=1}^{\lfloor \frac{m}{d} \rfloor}\lfloor \frac{m}{dy} \rfloor \]\(\mu\)的部分可以杜教筛处理，后边的部分可以整除分块话说这样不是变成\(O(n)\)了吗
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 2500005;
const int mod = 1e9+7;
int mu[N],prime[N],tot=0;
LL v[N];
void init(int n)
{
mu[1]=1;
for(int i=2;i<=n;i++)
{
if(!v[i])
{
v[i]=i;
mu[i]=-1;
prime[++tot]=i;
}
for(int j=1;j<=tot;j++)
{
if(prime[j]>v[i]||i*prime[j]>n) break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
else mu[i*prime[j]]=mu[i]*mu[prime[j]];
}
}
for(int i=1;i<=n;i++)
mu[i]=(mu[i]+mu[i-1]);
}
unordered_map<LL,LL> s;
LL Sum(LL n)
{
if(n<=N-10) return mu[n];
if(s[n]) return s[n];
LL res=1;
LL l=2,r;
for(;l<=n;l=r+1)
{
r=n/(n/l);
res=(res-1ll*(r-l+1)*Sum(n/l)%mod)%mod;
}
res=(res+mod)%mod;
s[n]=res;
return res;
}
LL S(LL n)
{
int l=1,r;
LL res=0;
for(;l<=n;l=r+1)
{
r=n/(n/l);
res=(res+(LL)(n/l)*(r-l+1)%mod)%mod;
}
return res;
}
int main()
{
freopen("math.in","r",stdin);
freopen("math.out","w",stdout);
init(N-10);
LL n;
cin>>n;
LL res=0;
int l=1,r;
for(;l<=n;l=r+1)
{
r=n/(n/l);
LL val=S(n/l);
res=(res+1ll*val*val%mod*(Sum(r)-Sum(l-1)+mod)%mod)%mod;
}
cout<<res;
return 0;
}