Today, I want to learn Mobius inversion, and then there is a problem that needs Du Jiao sieve in the problem list, and I pushed it to the place where there is only Du Jiao sieve left
So I had to learn this fairy algorithm. I really had to....
Du Jiao sieve is an advanced algorithm for prefix sum of quadrature functions below linear complexity, which is about \ (O(n^{\frac{2}{3}) \)
Pre knowledge
Integral function
Integral function: for any coprime integer \ (a,b \) with \ (f(ab)=f(a)f(b) \), it is called the number theoretic function of \ (f(x) \)
Complete integrable function: for any integer \ (a,b \), a number theoretic function with \ (f(ab)=f(a)f(b) \)
Common integrable functions: \ (\ varphi,\mu,\sigma,d \)
Common fully integrated functions: \ (\ epsilon,I,id \)
Here, we will explain what \ (\ epsilon,I,id \) means respectively: \ (\ epsilon(n)=[n=1],I(n)=1,id(n)=n \)
Dirichlet convolution
Let \ (f,g \) be two number theoretic functions whose Dirichlet convolution is: \ ((f*g)(n)=\sum\limits_{d|n}f(d)g(\frac{n}{d}) \)
Properties: satisfy commutative law, associative law
Unit element: \ (\ epsilon \) (i.e. \ (f*\epsilon=f \))
Several properties obtained by combining Dirichlet convolution:
\(\mu * I = \epsilon\)
\(\varphi * I = id\)
\(\mu * id = \varphi\)
Then be sure to remember these commonly used convolution functions and react at the key time
Du Jiao sieve
Let \ (f \) be an integral function. We now require \ (\ sum \ limits {I = 1} ^ {n} f (I) \)
It is obvious that the linear range can be directly screened, but the required range of general topics is \ (n\in [1,1e10] \)
In this way, we need algorithms with complexity below linearity
Set \ (s (n) = \ sum \ limits {I = 1} ^ {n}f (I) \)
Then find another integrable function \ (g \), roll it up with \ (f \), and consider the prefix and prefix of this function
\(\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}\)
Consider that \ (g(1)S(n) \) is equal to
\(\sum\limits_{i=1}^{n}g(i)S(\frac{n}{i})-\sum\limits_{i=2}^{n}g(i)S(\frac{n}{i})\)
Unless otherwise specified, the scores written in this article are rounded down
Then the previous formula can be reduced to \ (\ sum \ limits {I = 1} ^ {n} (f * g) (I) \)
So we get \ (\ huge {big routine} \)
\(g(1)S(n)=\sum\limits_{i=1}^{n}(f*g)(i)-\sum\limits_{i=2}^{n}g(i)S(\frac{n}{i})\)
In this case, the rest of our work is to find a suitable \ (g \) function
It is required that \ (\ sum \ limits {I = 1} ^ {n} (f * g) (I) \) can be given under the complexity of \ (O(1)orO(\sqrt{n}) \), otherwise the solution cannot be solved recursively
After finding \ (g \), the rest is number theory blocking, recursive solution and memorization
Put one below Template question
Mainly using the above
\(\varphi * I = id\)
as well as
\(\mu * I = \epsilon\)
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();}
The following is the Du Jiao screening problem we met when doing Mobius inversion today
\(\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}\)
Next, let's make \ (sum (n) = \ sum \ limits {I = 1} ^ {n} I \), \ (T=kd \)
There are,
\(\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}) \\ \end{aligned}\)
It is found that the following \ (\ sum\limits_{d|T}d\mu(\frac{T}{d}) \) is \ (id*\mu=\varphi \)
Then the original formula can be transformed into
\(\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)\\ \end{aligned}\)
The latter won't, and then look at the solution and find that we need to use Du Jiao sieve..
Therefore, after learning Du Jiao sieve, this problem has become a simple mathematical problem
Then we need to find an algorithm whose complexity is less than linear \ (\ sum \ limits {I = 1} ^ {n}i ^ 2 \ varphi (I) \)
It is not difficult to consider Du Jiao sieve, so you need to select a \ (g \), and it is found that you can select \ (g=id^2 \), that is \ (g(n)=n^2 \)
Then let \ (f(n)=n^2\varphi(n) \), then
\((f*g)(n)=\sum\limits_{d|n}d^2\varphi(d)\times (\frac{n}{d})^2=n^2\sum\limits_{d|n}\varphi(d)\)
Found the following \ (\ sum\limits_{d|n}\varphi(d)=(I*\varphi)(n)=id(n)=n \)
Then \ ((f*g)(n)=n^3 \), and then Du Jiao sifts out the prefix and Re number theory blocks
code#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();}