莫比烏斯反演+杜教篩
今天本身是要學習莫比烏斯反演的，然後題單裡面出現了一道需要用杜教篩的題，我還推到了隻剩杜教篩的地方
于是迫不得已學習了這種神仙算法，真是迫不得已啊。。。。
杜教篩是一種線上性複雜度以下的求積性函數字首和的進階算法，大概在\(O(n^{\frac{2}{3}})\)左右
前置知識
積性函數
積性函數：對于任意互質的整數\(a,b\)有\(f(ab)=f(a)f(b)\)則稱\(f(x)\)的數論函數
完全積性函數：對于任意整數\(a,b\)有\(f(ab)=f(a)f(b)\)的數論函數
常見的積性函數：\(\varphi,\mu,\sigma,d\)
常見的完全積性函數：\(\epsilon,I,id\)
這裡特殊解釋一下\(\epsilon,I,id\)分别是什麼意思：\(\epsilon(n)=[n=1],I(n)=1,id(n)=n\)
狄利克雷卷積
設\(f,g\)是兩個數論函數，它們的狄利克雷卷積卷積是：\((f*g)(n)=\sum\limits_{d|n}f(d)g(\frac{n}{d})\)
性質：滿足交換律，結合律
機關元：\(\epsilon\)（即\(f*\epsilon=f\)）
結合狄利克雷卷積得到的幾個性質：
\(\mu * I = \epsilon\)
\(\varphi * I = id\)
\(\mu * id = \varphi\)
然後一定要記住這些常用的卷積函數，關鍵時候要反應過來
杜教篩
設\(f\)為積性函數，我們現在要求\(\sum\limits_{i=1}^{n}f(i)\)
線性範圍一下顯然可以直接篩，但是一般題目要求範圍為\(n\in [1,1e10]\)
這樣我們需要複雜度線上性以下的算法
設\(S(n)=\sum\limits_{i=1}^{n}f(i)\)
那麼再找一個積性函數\(g\)，并把它和\(f\)卷起來，考慮這個函數的字首和
\(\begin{aligned} &= \sum\limits_{i=1}^{n} \sum \limits _{d|i} f(d)g(\frac{i}{d}) \\ &= \sum \limits _{d=1}^{n} g(d)\sum\limits _{i=1}^{\lfloor \frac{n}{d}\rfloor } f(i) \\ &= \sum \limits _{d=1}^{n} g(d) S(\lfloor \frac{n}{d} \rfloor) \end{aligned}\)
考慮\(g(1)S(n)\)就等于
\(\sum\limits_{i=1}^{n}g(i)S(\frac{n}{i})-\sum\limits_{i=2}^{n}g(i)S(\frac{n}{i})\)
若不特殊說明，本文中所寫的分數均是向下取整
那麼前面的式子可以化簡為\(\sum\limits_{i=1}^{n}(f*g)(i)\)
是以得到了\(\huge{大套路式}\)
\(g(1)S(n)=\sum\limits_{i=1}^{n}(f*g)(i)-\sum\limits_{i=2}^{n}g(i)S(\frac{n}{i})\)
這樣的話我們剩下的工作就是找到一個合适的\(g\)函數
要求盡量是可以在複雜度\(O(1)orO(\sqrt{n})\)的情況下出\(\sum\limits_{i=1}^{n}(f*g)(i)\)，否則無法遞歸求出解
找到\(g\)後剩下的就是數論分塊加遞歸求解加記憶化了
下面放一個模闆題
主要就是使用上述的
以及
code
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<unordered_map>
#define int long long
using namespace std;
const int NN=5000000;
namespace AE86{
auto read=[](){
int 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<<1)+(x<<3)+(ch^48);ch=getchar();}
return x*f;
};
}using namespace AE86;
int T,n,mu[NN+5],len,prime[NN+5];
int phi[NN+5];
auto getprime=[](){
mu[1]=1; phi[1]=1;
for(int i=2;i<=NN;i++){
if(!phi[i]) prime[++len]=i,mu[i]=-1,phi[i]=i-1;
for(int j=1;j<=len&&prime[j]*i<=NN;j++){
if(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<=NN;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
};
unordered_map<int,int>smu;
unordered_map<int,int>sphi;
inline int getmu(int n){
if(n<=NN) return mu[n];
if(smu.find(n)!=smu.end()) return smu[n];
int res=1,l=2,r;
while(l<=n){
r=n/(n/l);
res-=(r-l+1)*getmu(n/l);
l=r+1;
} return smu[n]=res;
}
inline int getphi(int n){
if(n<=NN) return phi[n];
if(sphi.find(n)!=sphi.end()) return sphi[n];
int res=(1+n)*n/2;int l=2,r;
while(l<=n){
r=n/(n/l);
res-=(r-l+1)*getphi(n/l);
l=r+1;
} return sphi[n]=res;
}
namespace WSN{
inline short main(){
T=read(); getprime();
while(T--){
n=read();
printf("%lld %lld\n",getphi(n),getmu(n));
}
return 0;
}
}
signed main(){return WSN::main();}
接下來就是今天做莫比烏斯反演時碰見的杜教篩題
簡單的數學題
\(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}ij\times gcd(i,j),(\mod{P})\)
\(\begin{aligned} & \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}ij\times gcd(i,j)\\
&=\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}id\times jd\times d[gcd(i,j)=1]\\
&=\sum\limits_{d=1}^{n}d^3\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}i\times j[gcd(i,j)=1] \\
&=\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}\sum\limits_{k|gcd(i,j)}\mu(k)ij \\
&=\sum\limits_{d=1}^{n}d^3\sum\limits_{k=1}^{\frac{n}{d}}\mu(k)k^2(1+2+...+\frac{\frac{n}{d}}{k})^2 \\
\end{aligned}\)
接下來我們令\(sum(n)=\sum\limits_{i=1}^{n}i\)，\(T=kd\)
則有，
\(\begin{aligned}&=\sum\limits_{d=1}^{n}d^3\sum\limits_{k=1}^{\frac{n}{d}}\mu(k)k^2sum^2(\frac{n}{T}) \\
&=\sum\limits_{T=1}^{n}sum^2(\frac{n}{T})\sum\limits_{d|T}d^3\mu(\frac{T}{d})(\frac{T}{d})^2 \\
&=\sum\limits_{T=1}^{n}sum^2(\frac{n}{T})T^2\sum\limits_{d|T}d\mu(\frac{T}{d}) \\
發現後面的\(\sum\limits_{d|T}d\mu(\frac{T}{d})\)是\(id*\mu=\varphi\)
那麼原式可以化為
\(\begin{aligned}&=\sum\limits_{T=1}^{n}sum^2(\frac{n}{T})T^2\varphi(T)\\
&=\sum\limits_{i=1}^{n}sum^2(\frac{n}{i})i^2\varphi(i)\\
後面的不會了，然後看題解發現要用杜教篩。。
是以學會杜教篩之後這題就成了簡單的數學題
然後我們需要找到一個複雜度線上性以下的算法計算\(\sum\limits_{i=1}^{n}i^2\varphi(i)\)
不難考慮杜教篩，那麼需要選擇一個\(g\)，發現可以選擇\(g=id^2\)，即\(g(n)=n^2\)
那麼設\(f(n)=n^2\varphi(n)\)，則
\((f*g)(n)=\sum\limits_{d|n}d^2\varphi(d)\times (\frac{n}{d})^2=n^2\sum\limits_{d|n}\varphi(d)\)
發現後面的\(\sum\limits_{d|n}\varphi(d)=(I*\varphi)(n)=id(n)=n\)
那麼\((f*g)(n)=n^3\)，然後杜教篩出字首和再數論分塊就可以解決了
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<unordered_map>
#define int long long
using namespace std;
const int NN=5000005;
int n,mod,v6,v2,ans;
namespace AE86{
auto read=[](){
int 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<<1)+(x<<3)+(ch^48);ch=getchar();}
return x*f;
};
auto write=[](int x,char opt='\n'){
char ch[20];short len=0;if(x<0)x=~x+1,putchar('-');
do{ch[len++]=x%10+(1<<5)+(1<<4);x/=10;}while(x);
for(short i=len-1;i>=0;--i)putchar(ch[i]);putchar(opt);
};
auto qmo=[](int a,int b,int ans=1){
int c=mod;for(;b;b>>=1,a=a*a%c)if(b&1)ans=ans*a%c;
return ans;
};
inline int sig(int n){return (1+n)%mod*n%mod*v2%mod;}
inline int pws(int n){return (2*n%mod+1)%mod*n%mod*(n+1)%mod*v6%mod;}
inline int pps(int n){return sig(n)*sig(n)%mod;}
}using namespace AE86;
namespace Prime{
signed len,prime[NN];
int phi[NN],F[NN];
bool vis[NN];
auto getprime=[](){
phi[1]=1; F[1]=1;
for(signed i=2;i<NN;++i){
if(!vis[i]) prime[++len]=i,phi[i]=i-1;
for(signed j=1;j<=len&&prime[j]*i<NN;++j){
vis[i*prime[j]]=1;
if(!(i%prime[j])){phi[i*prime[j]]=phi[i]*prime[j];break;}
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
F[i]=(phi[i]*i%mod*i%mod+F[i-1])%mod;
}
};
}using namespace Prime;
namespace dujiao_shai{
unordered_map<int,int>sF;
inline int getF(int n){
if(n<NN) return F[n];
if(sF[n]) return sF[n];
int res=pps(n);
for(int l=2,r;l<=n;l=r+1) r=n/(n/l),
res=(res-(pws(r)-pws(l-1)+mod)%mod*getF(n/l)%mod+mod)%mod;
res=(res%mod+mod)%mod;
return sF[n]=res;
}
}using namespace dujiao_shai;
namespace WSN{
inline short main(){
mod=read();n=read();getprime();
v6=qmo(6,mod-2);v2=qmo(2,mod-2);
for(int l=1,r;l<=n;l=r+1) r=n/(n/l),
ans=(ans+(getF(r)-getF(l-1)+mod)%mod*pps(n/l)%mod)%mod;
write(ans);
return 0;
}
}
signed main(){return WSN::main();}