Du Jiao sieve (basic chapter)

Posted by fredouille on Fri, 28 Jan 2022 09:35:11 +0100

Pre knowledge: Special topic of Mobius inversion (basic part)

Du Jiao sieve

Let \ (f(n) \) be an integral function
Find \ (S(n)=\sum_{i=1}^nf(i)\),\(n\leq 10^{10}\)
We consider prefixing and adding another integral function \ (f(n) \) on the \ (f(n) \) volume

\[\sum_{i=1}^n\sum_{d|i}g(d)f(\frac{i}{d}) \]

\[\sum_{d=1}^ng(d)\sum_{d|i}f(\frac{i}{d}) \]

\[\sum_{d=1}^ng(d)\sum_{d|i}f(\frac{i}{d}) \]

\[\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}f(i) \]

\[\sum_{d=1}^ng(d)S(\lfloor \frac{n}{d}\rfloor) \]

This leads to a sub problem
Next, let's focus on the term \ (d=1 \), that is \ (g(1)\times S(n) \)

\[g(1)S(n)=\sum_{i=1}^n\sum_{d|i}g(d)f(\frac{i}{d})-\sum_{d=2}^ng(d)S(\lfloor \frac{n}{d}\rfloor) \]

\[S(n)=\frac{\sum_{i=1}^n\sum_{d|i}g(d)f(\frac{i}{d})-\sum_{d=2}^ng(d)S(\lfloor \frac{n}{d}\rfloor)}{g(1)} \]

In this formula, \ (g(1) \) is easy to find, and the prefix sum of the convolution of \ (g \) and \ (f \) can also be easy to find by constructing an appropriate \ (g \), and the latter can be solved by integer division, blocking and recursion
If we preprocess the division of \ (S \) in advance, it can be proved that the complexity of calculating \ (S(n) \) is about \ (O(n^{\frac{2}{3}) \)

[BZOJ3944]--Sum

Most common

\[\sum_{i=1}^{n}\mu(i) \]

\[\sum_{i=1}^{n}\varphi(i) \]

First look at \ (mu(n) \), our task is to construct an appropriate \ (g \)
We know that \ (\ mu \times I = e \) (\ (\ times \) represents Dirichlet convolution)
So we construct \ (g(n)=I(n)=1 \)
Then the prefix sum of \ (g \) and \ (\ mu \) is still \ (e \), so the sum of this block is 1
in other words

\[S(n)=\frac{1-\sum_{d=2}^nS(\lfloor \frac{n}{d}\rfloor)}{1} \]

It can be solved by Du jiaosieve
Look at \ (\ varphi(n) \), we know \ (\ varphi(n)\times I=id \)
Then we also consider \ (g(n)=I(n)=1 \) on the volume
Then the value of the formula is

\[S(n)=\frac{\sum_{i=1}^ni-\sum_{d=2}^nS(\lfloor \frac{n}{d}\rfloor)}{1} \]

On the left is the sequence of equal differences, which is easy to sum

#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

The title gives a function \ (F(n)=n^2-3n+2 \)
And he guarantees that \ (F(n)=\sum_{d|n}f(d)\)
Find \ (\ sum {I = 1} ^ NF (I) \)
We first find the recurrence formula of \ (f(n) \)
Mobius can be obtained by inversion
\(f(n)=\sum_{d|n}F(d)\mu(\frac{n}{d})\)
This is obviously an integral function
However, the scope of \ (n \) is very large. Consider Du Jiao sieve
It is found that the \ (f (n) \) given by the topic can be regarded as the convolution of \ (f (n) \) and \ (I(n) \)
So the answer is

\[S(n)=\frac{\sum_{i=1}^ni^2-3i+2-\sum_{d=2}^nS(\lfloor \frac{n}{d}\rfloor)}{1} \]

The three items in the left half can be calculated separately, of which
\(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

Title Requirements

\[\sum_{i=1}^n\sum_{j=1}^nf(gcd(i,j))f(lcm(i,j)) \]

Where \ (f(n)=\sum_{d|n}\mu(d)d\)
At first, I wanted to reverse this formula, but I found that I couldn't do it anymore
Consider analyzing the properties of a seat \ (f(n) \)
This function is obviously an integral function, so we consider each prime factor separately
\(f(p^k)=1-p \), that is to say, for each prime factor, its index has nothing to do with it, so the \ (f \) of a number is only related to the number of types of other prime factors
Then we analyze \ (gcd(i,j) \) and \ (lcm(i,j) \)
1: For the common factors of \ (i,j \), they are multiplied twice, let them be A
2: For the respective factors of \ (i,j \), they are multiplied only once, assuming that they are B and C respectively
That is, A(BCA)
We can get (AB)(AC) by exchanging the order
That is \ (f(gcd(i,j))\times f(lcm(i,j))=f(i)\times f(j) \)
The formula becomes

\[\sum_{i=1}^n\sum_{j=1}^nf(i)f(j) \]

\[\sum_{i=1}^nf(i)\sum_{j=1}^nf(j) \]

\[(\sum_{i=1}^nf(i))^2 \]

The problem becomes to find the prefix and sum of \ (f \)
that is

\[\sum_{i=1}^n\sum_{d|i}\mu(d)d \]

\[\sum_{d=1}^n\mu(d)d\sum_{d|i}^n1 \]

\[\sum_{d=1}^n\mu(d)d\lfloor \frac{n}{d}\rfloor \]

The second half can be divided into blocks
The problem is further transformed into finding the prefix and sum of \ (\ mu(d)d \)
Consider Du Jiao sieve
We consider constructing a \ (g(n)=Id(n)=n \)
That is, $$g\times (\mu(d)d)$$

\[=\sum_{d|n}\;\mu(d)d\;g(\frac{n}{d}) \]

\[=\sum_{d|n}\;\mu(d)d\;\frac{n}{d} \]

\[=n\sum_{d|n}\;\mu(d) \]

\[=e \]

So the prefix and are

\[S(n)=\frac{1-\sum_{d=2}^niS(\lfloor \frac{n}{d}\rfloor)}{1} \]

Just take it in

#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 sum of least common multiples V3

cancer
seek

\[\sum_{i=1}^n\sum_{j=1}^nlcm(i,j) \]

For more details, see Special topic of Mobius inversion (basic part) [BZOJ2693] jzptab in
Let's use the last formula directly here

\[= \sum_{T=1}^nS(\lfloor \frac{n}{T} \rfloor)S(\lfloor \frac{m}{T} \rfloor)T\sum_{k|T}^T\mu(k)k \]

Where \ (S(n)=\sum_{i=1}^ni\)
So the key problem is to screen out \ (f (n) = n \ sum_ Prefix and of {d|n} \ mu (d) d \)

\[\sum_{i=1}^n\sum_{d|i}i\mu(d)d \]

\[\sum_{d=1}^n\mu(d)d\sum_{d|i}i \]

\[\sum_{d=1}^n\mu(d)d\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}id \]

\[\sum_{d=1}^n\mu(d)d^2\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}i \]

Obviously, in the front part, we can roll up a \ (g(n)=n^2 \) to teach the sieve. The whole formula can be divided into blocks
However, don't forget that there is an integer division block outside, so it becomes \ (O(n) \)
So we can only consider other methods
By observing and typing tables, we find that \ (f\times (n^2) = (n) \)
As for the specific proof, I'm sorry I don't understand, so we put a \ (g(n)=n^2 \) on the volume
Prefix and change to

\[S(n)=\frac{\sum_{i=1}^ni-\sum_{d=2}^ni^2S(\lfloor \frac{n}{d}\rfloor)}{1} \]

It can be solved

#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] sum of maximum common divisor V3

See the above for details Special topic of Mobius inversion (basic part) [NOI2010day1]T1 energy acquisition in
The final formula is

\[=\sum_{d=1}^n\varphi(d)\lfloor \frac{n}{d} \rfloor\lfloor \frac{m}{d} \rfloor \]

Just process the prefix and of \ (\ varphi(n) \) directly

#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's number theory

See the above for details Special topic of Mobius inversion (basic part) The sum of approximate numbers of [SDOI2015] in

\[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 \) can be processed by Du Jiao sieve, and the rear part can be divided into blocks
In other words, didn't it become \ (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;
}