杜教筛(进阶篇)

一道更比一道毒瘤

[51 nod 1227] 平均最小公倍数

其实就是求
a n s = A n s ( b ) − A n s ( a − 1 ) ans=Ans(b)-Ans(a-1) ans=Ans(b)Ans(a1)
因此我们只需要求出函数 A n s ( n ) Ans(n) Ans(n)就行了
A n s ( n ) = ∑ i = 1 n 1 i ∑ j = 1 i l c m ( j , i ) Ans(n)=\sum_{i=1}^n\frac{1}{i}\sum_{j=1}^ilcm(j,i) Ans(n)=i=1ni1j=1ilcm(j,i)
= ∑ i = 1 n 1 i ∑ j = 1 i j i g c d ( j , i ) =\sum_{i=1}^n\frac{1}{i}\sum_{j=1}^i\frac{ji}{gcd(j,i)} =i=1ni1j=1igcd(j,i)ji
= ∑ i = 1 n ∑ j = 1 i j g c d ( j , i ) =\sum_{i=1}^n\sum_{j=1}^i\frac{j}{gcd(j,i)} =i=1nj=1igcd(j,i)j
= ∑ i = 1 n ∑ j = 1 i ∑ d = 1 [ g c d ( i , j ) = = d ] j d =\sum_{i=1}^n\sum_{j=1}^i\sum_{d=1}[gcd(i,j)==d]\frac{j}{d} =i=1nj=1id=1[gcd(i,j)==d]dj
= ∑ d = 1 n 1 d ∑ d ∣ i n ∑ d ∣ j i [ g c d ( i , j ) = = d ] j =\sum_{d=1}^n\frac{1}{d}\sum_{d|i}^n\sum_{d|j}^i[gcd(i,j)==d]j =d=1nd1dindji[gcd(i,j)==d]j
= ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 i [ g c d ( i , j ) = = 1 ] j =\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^i[gcd(i,j)==1]j =d=1ni=1dnj=1i[gcd(i,j)==1]j
= ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 i ∑ k ∣ i , k ∣ j μ ( k ) j =\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^i\sum_{k|i,k|j}\mu(k)j =d=1ni=1dnj=1iki,kjμ(k)j

F ( n ) = ∑ i = 1 n ∑ k ∣ i , k ∣ n μ ( k ) i F(n)=\sum_{i=1}^n\sum_{k|i,k|n}\mu(k)i F(n)=i=1nki,knμ(k)i
F ( n ) = ∑ k ∣ n μ ( k ) k ∑ i = 1 n k i F(n)=\sum_{k|n}\mu(k)k\sum_{i=1}^{\frac{n}{k}}i F(n)=knμ(k)ki=1kni
F ( n ) = ∑ k ∣ n μ ( k ) k ( 1 + n k ) n k 1 2 F(n)=\sum_{k|n}\mu(k)k(1+\frac{n}{k})\frac{n}{k}\frac{1}{2} F(n)=knμ(k)k(1+kn)kn21
F ( n ) = 1 2 ∑ k ∣ n μ ( k ) ( 1 + n k ) n F(n)=\frac{1}{2}\sum_{k|n}\mu(k)(1+\frac{n}{k})n F(n)=21knμ(k)(1+kn)n
F ( n ) = 1 2 ( ∑ k ∣ n μ ( k ) n + ∑ d ∣ n μ ( k ) n n k ) F(n)=\frac{1}{2}(\sum_{k|n}\mu(k)n+\sum_{d|n}\mu(k)n\frac{n}{k}) F(n)=21(knμ(k)n+dnμ(k)nkn)
而我们知道
∑ k ∣ n μ ( k ) n = [ n = = 1 ] \sum_{k|n}\mu(k)n=[n==1] knμ(k)n=[n==1]
我们也知道
∑ d ∣ n μ ( k ) n k = μ ∗ i d = φ \sum_{d|n}\mu(k)\frac{n}{k}=\mu*id=\varphi dnμ(k)kn=μid=φ
其中 ∗ * 表示狄利克雷卷积
因此变成
F ( n ) = 1 2 ( e + n φ ( n ) ) F(n)=\frac{1}{2}(e+n\varphi(n)) F(n)=21(e+nφ(n))
带进整个式子
= 1 2 ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ( e + i φ ( i ) ) =\frac{1}{2}\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}(e+i\varphi(i)) =21d=1ni=1dn(e+iφ(i))
考虑把 e e e拆出来,那么只有 i i i等于1是会有贡献,一共有n个1,所以总贡献是n
= 1 2 ( n + ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ i × φ ( i ) ) =\frac{1}{2}(n+\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}i\times \varphi(i)) =21(n+d=1ni=1dni×φ(i))
问题转化为求 i × φ ( i ) i\times \varphi(i) i×φ(i)的前缀和,我们考虑卷上一个 i d ( n ) = n id(n)=n id(n)=n,化简可以得到 n 2 n^2 n2,上杜教筛就行了
S ( n ) = ∑ i = 1 n n 2 − ∑ i = 2 n i S ( n / i ) S(n)=\sum_{i=1}^nn^2-\sum_{i=2}^niS(n/i) S(n)=i=1nn2i=2niS(n/i)
复杂度 O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32)

#include<bits/stdc++.h>
using namespace std;
const int N = 1e6+7;
typedef long long LL;
const LL mod = 1e9+7;
LL phi[N];
LL f[N];
int v[N],prime[N],tot=0;
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 inv6=Pow(6,mod-2);
LL inv2=Pow(2,mod-2);
LL sqr(LL n)
{
	n%=mod;
	return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
LL sum(LL n)
{
	return 1ll*(1+n)*n%mod*inv2%mod;
}
void init(int n)
{
	phi[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=i;
			prime[++tot]=i;
			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)
			{
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			else phi[i*prime[j]]=phi[i]*(prime[j]-1);
		}
	}
	for(int i=1;i<=n;i++)
	phi[i]=(phi[i-1]+phi[i]*i%mod)%mod;
}
unordered_map<LL,LL> s,vis;
LL Sum(LL n)
{
	if(n<=1e6) return phi[n];
	if(vis[n]) return s[n];
	LL res=sqr(n);
	LL l=2,r;
	for(;l<=n;l=r+1)
	{
		r=(n/(n/l));
		res=(res-(sum(r)-sum(l-1)+mod)%mod*Sum(n/l)%mod+mod)%mod;
	}
	vis[n]=1;
	s[n]=res;
	return res;
}
LL calc(LL n)
{
	LL res=n;
	LL l=1,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=(res+1ll*(r-l+1)%mod*Sum(n/l)%mod)%mod;
	}
	return res*inv2%mod;	
}
int main()
{
	freopen("minave.in","r",stdin);
	freopen("minave.out","w",stdout);
	init(1e6);
	LL a,b;
	cin>>a>>b;
	cout<<(calc(b)-calc(a-1)+mod)%mod;
	return 0;
} 

SP20173 DIVCNT2 - Counting Divisors (square)

考虑先分析 σ ( n 2 ) \sigma(n^2) σ(n2)的性质
显然这是一个积性函数,我们考虑从每一个质因子分别考虑
σ ( ( p k ) 2 ) = 2 ∗ k + 1 = k + ( k + 1 ) = σ ( p k ) + σ ( p k − 1 ) \sigma((p^k)^2)=2*k+1=k+(k+1)=\sigma(p^k)+\sigma(p^{k-1}) σ((pk)2)=2k+1=k+(k+1)=σ(pk)+σ(pk1)
我们继续观察, p k p^k pk? p k − 1 p^{k-1} pk1?
也就是 p k p_k pk除以他们的质因子的次幂是0或1
那么有没有什么是和指数相关的的呢?
没错,就是 μ \mu μ,但是 μ \mu μ有负数,怎么办呢
没错,加个平方就行了
事实上和我们分析的一样, σ ( n 2 ) = ∑ d ∣ n σ ( d ) μ 2 ( n d ) \sigma(n^2)=\sum_{d|n}\sigma(d)\mu^2(\frac{n}{d}) σ(n2)=dnσ(d)μ2(dn)

当然和 σ ( n 2 ) = ∑ d ∣ n σ ( n d ) μ 2 ( d ) \sigma(n^2)=\sum_{d|n}\sigma(\frac{n}{d})\mu^2(d) σ(n2)=dnσ(dn)μ2(d)是一样的
a n s = ∑ i = 1 n ∑ d ∣ n σ ( i d ) μ 2 ( i ) ans=\sum_{i=1}^n\sum_{d|n}\sigma(\frac{i}{d})\mu^2(i) ans=i=1ndnσ(di)μ2(i)
莫反一下
a n s = ∑ d = 1 n μ 2 ( d ) ∑ d ∣ i σ ( i d ) ans=\sum_{d=1}^n\mu^2(d)\sum_{d|i}\sigma(\frac{i}{d}) ans=d=1nμ2(d)diσ(di)
a n s = ∑ d = 1 n μ 2 ( d ) ∑ i = 1 ⌊ n d ⌋ σ ( i ) ans=\sum_{d=1}^n\mu^2(d)\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sigma(i) ans=d=1nμ2(d)i=1dnσ(i)
这个东西可以用整除分块做一下,现在的问题变为得到
μ 2 ( n ) \mu^2(n) μ2(n) σ ( n ) \sigma(n) σ(n)的前缀和
根据一些容斥的知识,我们可以得到
∑ i = 1 n μ 2 ( i ) = ∑ i = 1 n μ ( i ) ⌊ n d 2 ⌋ \sum_{i=1}^n\mu^2(i)=\sum_{i=1}^{\sqrt n}\mu(i)\lfloor \frac{n}{d^2}\rfloor i=1nμ2(i)=i=1n μ(i)d2n
这个可以用整除分块做到 O ( n ) O(\sqrt n) O(n )
而有一个众所众知的式子
∑ i = 1 n σ ( i ) = ∑ i = 1 n ⌊ n i ⌋ \sum_{i=1}^n\sigma(i)=\sum_{i=1}^n\lfloor \frac{n}{i}\rfloor i=1nσ(i)=i=1nin
证明的话就是考虑计算每个数在n个数中有多少倍数
这个也可以整除分块做到 O ( n ) O(\sqrt n) O(n )
不过两个 n n n套起来就是 O ( n ) O(n) O(n)
我们考虑先预处理 μ 2 ( n ) \mu^2(n) μ2(n) σ ( n ) \sigma(n) σ(n)的前 n 2 3 n^{\frac{2}{3}} n32项的前缀和,询问的时候直接调用
复杂度和杜教筛类似是 O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32)
当然,这个题有毒
它的 n n n 1 0 12 10^{12} 1012,预处理要处理到 1 0 8 10^8 108
因此不仅空间大,常数也大
这份代码在SPOJ上并不能跑过

#include<bits/stdc++.h>
using namespace std;
const int N = 5e7+5;
const int M = 1e7+7;
typedef long long LL;
LL d[N];
int mu[N];
int t[N];
bool v[N];
int prime[M],tot=0;
void init(int n)
{
	mu[1]=1;
	d[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=1;
			prime[++tot]=i;
			d[i]=2;
			mu[i]=-1;
			t[i]=2;
		}
		for(int j=1;j<=tot;j++)
		{
			if(i*prime[j]>n) break;
			v[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				t[i*prime[j]]=t[i]+1;
				d[i*prime[j]]=d[i]/t[i]*(t[i]+1);
				mu[i*prime[j]]=0;
				break;
			}
				t[i*prime[j]]=2;
				d[i*prime[j]]=d[i]*d[prime[j]];
				mu[i*prime[j]]=mu[i]*mu[prime[j]];
		}
	}
	for(int i=1;i<=n;i++)
	{
		d[i]+=d[i-1];
		t[i]=t[i-1]+abs(mu[i]);
	}
}
LL R;
LL SumU(LL n)
{
	if(n<=R) return t[n];
	LL res=0;
	for(LL i=1;i*i<=n;i++)
	res=(res+mu[i]*(n/i/i));
	return res;
}
LL SumD(LL n)
{
	if(n<=R) return d[n];
	LL res=0;
	LL l=1,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=res+(r-l+1)*(n/l);		
	}
	return res;
}
void solve(LL n)
{
	LL res=0;
	LL l=1,r;
	LL last=0;
	LL m=sqrt(n);
	for(LL i=1;i<=m;i++)
	if(mu[i]!=0) res=res+SumD(n/i);
	l=m+1;last=SumU(m);
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		LL now=SumU(r);
		res=(res+(now-last)*SumD(n/l));
		last=now;
	}
	printf("%lld\n",res);
}
LL a[20000];
LL INF = 1e12;

int main()
{
	freopen("a.in","r",stdin);
	freopen("a.out","w",stdout);
	int T;
	cin>>T;
	LL r;
	for(int i=1;i<=T;i++)
	{
		scanf("%lld",&a[i]);
		r=max(r,a[i]);	
	}
	if(r<=10000) R=10000;
	else R=N-10;
	init(R);
	for(int i=1;i<=T;i++)
	solve(a[i]);
	return 0;
}

[51nod1222]最小公倍数计数

还是先做成前缀相减的形式
问题转化为求
F ( n ) = ∑ i = 1 n ∑ j = 1 n [ l c m ( i , j ) ≤ n ] F(n)=\sum_{i=1}^n\sum_{j=1}^n[lcm(i,j)\leq n] F(n)=i=1nj=1n[lcm(i,j)n]
F ( n ) = ∑ i = 1 n ∑ j = 1 n ∑ d ∣ i , d ∣ j [ g c d ( i , j ) = = d ] [ i j d ≤ n ] F(n)=\sum_{i=1}^n\sum_{j=1}^n\sum_{d|i,d|j}[gcd(i,j)==d][\frac{ij}{d}\leq n] F(n)=i=1nj=1ndi,dj[gcd(i,j)==d][dijn]
= ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ n d ⌋ [ g c d ( i , j ) = = 1 ] [ i j ≤ n d ] =\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor }\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor}[gcd(i,j)==1][ij\leq \frac{n}{d}] =d=1ni=1dnj=1dn[gcd(i,j)==1][ijdn]
= ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ n d ⌋ ∑ k ∣ i , k ∣ j μ ( k ) [ i j ≤ n d ] =\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor }\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{k|i,k|j}\mu(k)[ij\leq \frac{n}{d}] =d=1ni=1dnj=1dnki,kjμ(k)[ijdn]
= ∑ d = 1 n ∑ k = 1 ⌊ n d ⌋ μ ( k ) ∑ i = 1 ⌊ n d k ⌋ ∑ j = 1 ⌊ n d k ⌋ [ i j k 2 ≤ n d ] =\sum_{d=1}^n\sum_{k=1}^{\lfloor \frac{n}{d}\rfloor }\mu(k)\sum_{i=1}^{\lfloor \frac{n}{dk} \rfloor }\sum_{j=1}^{\lfloor \frac{n}{dk}\rfloor}[ijk^2\leq \frac{n}{d}] =d=1nk=1dnμ(k)i=1dknj=1dkn[ijk2dn]
= ∑ k = 1 n μ ( k ) ∑ d = 1 ⌊ n k ⌋ ∑ i = 1 ⌊ n d k ⌋ ∑ j = 1 ⌊ n d k ⌋ [ i j d ≤ n k 2 ] =\sum_{k=1}^n\mu(k)\sum_{d=1}^{\lfloor \frac{n}{k}\rfloor }\sum_{i=1}^{\lfloor \frac{n}{dk} \rfloor }\sum_{j=1}^{\lfloor \frac{n}{dk}\rfloor}[ijd\leq \frac{n}{k^2}] =k=1nμ(k)d=1kni=1dknj=1dkn[ijdk2n]
发现 k k k最大不会超过 ( n ) \sqrt(n) ( n)
因此
= ∑ k = 1 ( n ) μ ( k ) ∑ d = 1 ⌊ n k ⌋ ∑ i = 1 ⌊ n d k ⌋ ∑ j = 1 ⌊ n d k ⌋ [ i j d ≤ n k 2 ] =\sum_{k=1}^{\sqrt(n)}\mu(k)\sum_{d=1}^{\lfloor \frac{n}{k}\rfloor }\sum_{i=1}^{\lfloor \frac{n}{dk} \rfloor }\sum_{j=1}^{\lfloor \frac{n}{dk}\rfloor}[ijd\leq \frac{n}{k^2}] =k=1( n)μ(k)d=1kni=1dknj=1dkn[ijdk2n]
可以考虑枚举 k k k,计算后面的值
观察后三项的上指标
会发现,如果 a ≥ ⌊ n k ⌋ a\geq \lfloor \frac{n}{k} \rfloor akn,后面的式子就为0
类似的, i , j i,j i,j的上指标都可以被后面的约束掉
因此,后边部分等于
∑ d ∑ i ∑ j [ i j k ≤ n k ] \sum_{d}\sum_{i}\sum_{j}[ijk\leq\frac{n}{k}] dij[ijkkn]
我们考虑强制让 d ≤ i ≤ j d\leq i \leq j dij 然后乘一个系数就行了
考虑不同的 d , i , j d,i,j d,i,j的状态
1: d , d , d d,d,d d,d,d,即三个相等,枚举到的合法的d就加一就行了
2: d , i , i d,i,i d,i,i,即后两个相等,算出合法的 d , i d,i d,i的数量乘以3
3: d , d , i d,d,i d,d,i 与2类似
4: d , i , j d,i,j d,i,j
发现 d d d只需要枚举到n的三次根就行了,i只需要枚举到n的平方根就行了
总复杂度与杜教筛类似为 O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32),它的系数是6
可以再计算后面的时候一块计算

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N =1e6+7;
LL mu[N],v[N];
LL prime[N],tot=0;
void init(LL n)
{
	mu[1]=1;
	for(LL i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=i;
			mu[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)
			{
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=mu[i]*mu[prime[j]];
		}
	}
}
LL S(LL n)
{
	LL res=0;
	for(LL k=1;k*k<=n;k++)
	{
		if(!mu[k]) continue;
		LL cnt=0;
		LL limit=n/(k*k);
		for(LL i=1;i*i*i<=limit;i++)
		{
			for(LL j=i+1;i*j*j<=limit;j++)
			cnt+=(LL)(limit/(i*j)-j)*6+3;//i,j,k+i,j,j 
			cnt+=(LL)(limit/(i*i)-i)*3;//i,i,j
			cnt++;//i,i,i
		}
		res=res+mu[k]*cnt;
	}
	return (res+n)/2;
}
int main()
{
	freopen("a.in","r",stdin);
	freopen("a.out","w",stdout);
	init(1e6);
	LL l,r;
	cin>>l>>r;
	cout<<S(r)-S(l-1);
	return 0;
}

[51nod1220] 约数之和


∑ i = 1 n ∑ j = 1 n σ ( i j ) \sum_{i=1}^n\sum_{j=1}^n\sigma(ij) i=1nj=1nσ(ij)
其中 σ ( n ) \sigma(n) σ(n)表示 n n n的因子之和
首先有引理
σ ( i j ) = ∑ x ∣ i ∑ y ∣ j [ g c d ( x , y ) = = 1 ] x j y \sigma(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]\frac{xj}{y} σ(ij)=xiyj[gcd(x,y)==1]yxj
证明与SDOI约数个数和类似,不赘述了
直接莫反
∑ i = 1 n ∑ j = 1 n ∑ x ∣ i ∑ y ∣ j [ g c d ( x , y ) = = 1 ] x j y \sum_{i=1}^n\sum_{j=1}^n\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]\frac{xj}{y} i=1nj=1nxiyj[gcd(x,y)==1]yxj
∑ x = 1 n x ∑ y = 1 n ∑ i = 1 ⌊ n x ⌋ ∑ j = 1 ⌊ n y ⌋ [ g c d ( x , y ) = = 1 ] j \sum_{x=1}^nx\sum_{y=1}^n\sum_{i=1}^{\lfloor \frac{n}{x}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{y}\rfloor}[gcd(x,y)==1]j x=1nxy=1ni=1xnj=1yn[gcd(x,y)==1]j
∑ x = 1 n x ∑ y = 1 n [ g c d ( x , y ) = = 1 ] ⌊ n x ⌋ ∑ j = 1 ⌊ n y ⌋ j \sum_{x=1}^nx\sum_{y=1}^n[gcd(x,y)==1]\lfloor \frac{n}{x}\rfloor\sum_{j=1}^{\lfloor \frac{n}{y}\rfloor}j x=1nxy=1n[gcd(x,y)==1]xnj=1ynj
∑ x = 1 n x ∑ y = 1 n ∑ d ∣ x , d ∣ y μ ( d ) ⌊ n x ⌋ ∑ j = 1 ⌊ n y ⌋ j \sum_{x=1}^nx\sum_{y=1}^n\sum_{d|x,d|y}\mu(d)\lfloor \frac{n}{x}\rfloor\sum_{j=1}^{\lfloor \frac{n}{y}\rfloor}j x=1nxy=1ndx,dyμ(d)xnj=1ynj
∑ d = 1 n μ ( d ) d ∑ x = 1 ⌊ n d ⌋ x ∑ y = 1 ⌊ n d ⌋ ⌊ n d x ⌋ ∑ j = 1 ⌊ n d y ⌋ j \sum_{d=1}^n\mu(d)d\sum_{x=1}^{\lfloor \frac{n}{d}\rfloor}x\sum_{y=1}^{\lfloor \frac{n}{d}\rfloor}\lfloor \frac{n}{dx}\rfloor\sum_{j=1}^{\lfloor \frac{n}{dy}\rfloor}j d=1nμ(d)dx=1dnxy=1dndxnj=1dynj
∑ d = 1 n μ ( d ) d ∑ x = 1 ⌊ n d ⌋ x ⌊ n d x ⌋ ∑ y = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ n d y ⌋ j \sum_{d=1}^n\mu(d)d\sum_{x=1}^{\lfloor \frac{n}{d}\rfloor}x\lfloor \frac{n}{dx}\rfloor\sum_{y=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{dy}\rfloor}j d=1nμ(d)dx=1dnxdxny=1dnj=1dynj
f ( n ) = ∑ i = 1 n ∑ j = 1 ⌊ n i ⌋ i f(n)=\sum_{i=1}^n\sum_{j=1}^{\lfloor \frac{n}{i}\rfloor}i f(n)=i=1nj=1ini, g ( n ) = ∑ i = 1 n i ⌊ n i ⌋ g(n)=\sum_{i=1}^ni\lfloor \frac{n}{i} \rfloor g(n)=i=1niin
a n s = ∑ d = 1 n μ ( d ) d × g ( ⌊ n d ⌋ ) f ( ⌊ n d ⌋ ) ans=\sum_{d=1}^n\mu(d)d\times g(\lfloor \frac{n}{d}\rfloor)f(\lfloor \frac{n}{d}\rfloor) ans=d=1nμ(d)d×g(dn)f(dn)
考虑 f ( n ) f(n) f(n) g ( n ) g(n) g(n)的关系
会发现, g ( n ) g(n) g(n)其实就是把 f ( n ) f(n) f(n)中每个数 i i i的权值乘上他出现的次数
因此两个式子是等价的,问题进一步转化为
a n s = ∑ d = 1 n μ ( d ) d × g ( ⌊ n d ⌋ ) 2 ans=\sum_{d=1}^n\mu(d)d\times g(\lfloor \frac{n}{d}\rfloor)^2 ans=d=1nμ(d)d×g(dn)2
我们只需要求出 g ( n ) g(n) g(n)的前缀和就可以整除分块了
接着会发现 ∑ i = 1 n σ ( i ) = ∑ i = 1 n i ⌊ n i ⌋ = g ( n ) \sum_{i=1}^n\sigma(i)=\sum_{i=1}^ni\lfloor \frac{n}{i} \rfloor=g(n) i=1nσ(i)=i=1niin=g(n)
因此 g ( n ) = ∑ i = 1 n σ ( i ) g(n)=\sum_{i=1}^n\sigma(i) g(n)=i=1nσ(i)
众所众知, σ ( n ) \sigma(n) σ(n)是积性函数,可以线性筛求出来
类似于杜教筛,预处理除 n 2 3 n^{\frac{2}{3}} n32的前缀和,剩下的整除分块暴力做,可以做到 O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32)
别忘了最后还剩一个 μ ( d ) d \mu(d)d μ(d)d,卷上一个 i d ( n ) = n id(n)=n id(n)=n,就可以杜教筛了

#include<bits/stdc++.h>
using namespace std;
const int N = 1e6+7;
typedef long long LL;
int mu[N];
LL f[N];
int v[N],prime[N],tot=0;
const int mod = 1e9+7;
LL d[N],t[N],g[N];
void init(LL n)
{
	mu[1]=1;
	d[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=i;
			prime[++tot]=i;
			mu[i]=-1;
			d[i]=i+1;
			g[i]=i;
			t[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)
			{
				g[i*prime[j]]=g[i]*prime[j];
				t[i*prime[j]]=t[i]+g[i*prime[j]];
				d[i*prime[j]]=d[i]/t[i]*t[i*prime[j]];
				mu[i*prime[j]]=0;
				break;
			}
			else 
			{
				mu[i*prime[j]]=mu[i]*mu[prime[j]];
				d[i*prime[j]]=d[i]*d[prime[j]];
				g[i*prime[j]]=prime[j];
				t[i*prime[j]]=1+prime[j];
			}
		} 
	}
	for(int i=1;i<=n;i++)
	{
		f[i]=(f[i-1]+mu[i]*i%mod+mod)%mod;
		d[i]=(d[i]+d[i-1])%mod;
	}
}
unordered_map<LL,LL> vis,s;
LL Sum(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-1ll*(l+r)*(r-l+1)/2%mod*Sum(n/l)%mod+mod)%mod;
	}
	vis[n]=1;
	s[n]=res;
	return res;
}
LL D(LL n)
{
	if(n<=1e6) return d[n];
	LL l=1,r;
	LL res=0;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=(res+(l+r)*(r-l+1)/2%mod*(n/l)%mod)%mod;
	}
	return res;
}
int main()
{
	freopen("a.in","r",stdin);
	freopen("a.out","w",stdout);
	init(1e6);
	LL n;
	cin>>n;
	LL res=0;
	LL l=1,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		LL val=D(n/l);
		res=(res+(Sum(r)-Sum(l-1)+mod)%mod*val%mod*val%mod)%mod;
	}
	cout<<res;
	return 0;
}

[51nod1584]加权约数和


∑ i = 1 n ∑ j = 1 n m a x ( i , j ) σ ( i j ) \sum_{i=1}^n\sum_{j=1}^nmax(i,j)\sigma(ij) i=1nj=1nmax(i,j)σ(ij)
我们知道这个 n ∗ n n*n nn的答案矩阵是对称的,因此答案即为
2 ∑ i = 1 n ∑ j = 1 i i σ ( i j ) − ∑ i = 1 n σ ( i 2 ) 2\sum_{i=1}^n\sum_{j=1}^ii\sigma(ij)-\sum_{i=1}^n\sigma(i^2) 2i=1nj=1iiσ(ij)i=1nσ(i2)
先看前半部分
和上一题类似
σ ( i j ) = ∑ x ∣ i ∑ y ∣ j [ g c d ( x , y ) = = 1 ] x j y \sigma(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]\frac{xj}{y} σ(ij)=xiyj[gcd(x,y)==1]yxj
带进式子
∑ i = 1 n i ∑ j = 1 i ∑ x ∣ i ∑ y ∣ j [ g c d ( x , y ) = = 1 ] x j y \sum_{i=1}^ni\sum_{j=1}^i\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]\frac{xj}{y} i=1nij=1ixiyj[gcd(x,y)==1]yxj
∑ i = 1 n i ∑ j = 1 i ∑ x ∣ i x ∑ y ∣ j ∑ d ∣ x , d ∣ y μ ( d ) j y \sum_{i=1}^ni\sum_{j=1}^i\sum_{x|i}x\sum_{y|j}\sum_{d|x,d|y}\mu(d)\frac{j}{y} i=1nij=1ixixyjdx,dyμ(d)yj
∑ d = 1 n μ ( d ) d ∑ i = 1 ⌊ n d ⌋ i ∑ x ∣ i x d ∑ j = 1 i ∑ y ∣ j j y \sum_{d=1}^n\mu(d)d\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}i\sum_{x|i}xd\sum_{j=1}^i\sum_{y|j}\frac{j}{y} d=1nμ(d)di=1dnixixdj=1iyjyj
∑ d = 1 n μ ( d ) d ∑ i = 1 ⌊ n d ⌋ i ∑ x ∣ i x d ∑ j = 1 i ∑ y ∣ j y \sum_{d=1}^n\mu(d)d\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}i\sum_{x|i}xd\sum_{j=1}^i\sum_{y|j}y d=1nμ(d)di=1dnixixdj=1iyjy
∑ d = 1 n μ ( d ) d 2 ∑ i = 1 ⌊ n d ⌋ i σ ( i ) ∑ j = 1 i σ ( j ) \sum_{d=1}^n\mu(d)d^2\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}i\sigma(i)\sum_{j=1}^i\sigma(j) d=1nμ(d)d2i=1dniσ(i)j=1iσ(j)
f ( n ) = n σ ( n ) ∑ i = 1 n σ ( i ) f(n)=n\sigma(n)\sum_{i=1}^n\sigma(i) f(n)=nσ(n)i=1nσ(i)
a n s = ∑ d = 1 n μ ( d ) d 2 ∑ i = 1 ⌊ n d ⌋ f ( i ) ans=\sum_{d=1}^n\mu(d)d^2\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}f(i) ans=d=1nμ(d)d2i=1dnf(i)
我们可以预处理除 σ ( n ) \sigma(n) σ(n)的前缀和,这样就能快速计算 f ( i ) f(i) f(i)
但是询问的次数非常多,需要我们 O ( 1 ) O(1) O(1)回答
因此继续化简
a n s = ∑ T = 1 n ∑ k ∣ T μ ( k ) k 2 f ( T k ) ans=\sum_{T=1}^n\sum_{k|T}\mu(k)k^2f(\frac{T}{k}) ans=T=1nkTμ(k)k2f(kT)
这个式子可以通过调和级数 O ( n log ⁡ n ) O(n\log n) O(nlogn)的算出来,再计算前缀和就可以 O ( 1 ) O(1) O(1)查询了
接着看后半部分
问题是处理 σ ( n 2 ) \sigma(n^2) σ(n2)的前缀和
这个也是可以用线性筛筛出来的,毕竟这是一个积性函数

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N = 1e6+7;
const LL mod =1e9+7; 
LL d[N],t[N],p[N];
LL sd[N],st[N],sp[N];
LL s[N],g[N];
LL mu[N];
LL v[N],prime[N],tot=0;
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 f[N];
void init(LL n)
{
	mu[1]=1;
	sd[1]=1;
	d[1]=1;
	for(LL i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=i;
			prime[++tot]=i;
			mu[i]=-1;
			d[i]=i+1;
			t[i]=i+1;
			p[i]=i;
			sd[i]=(1+i+(LL)i*i%mod)%mod;
			st[i]=(1+i+(LL)i*i%mod)%mod;
			sp[i]=(LL)i*i%mod;
		}
		for(LL j=1;j<=tot;j++)
		{
			if(v[i]<prime[j]||i*prime[j]>n) break;
			LL k=i*prime[j];
			v[i*prime[j]]=prime[j];
			if(i%prime[j]==0)
			{
				mu[k]=0;
				p[k]=p[i]*prime[j]%mod;
				t[k]=(t[i]+p[k])%mod;
				d[k]=d[i]*inv(t[i])%mod*t[k]%mod;
				sp[k]=sp[i]*prime[j]%mod*prime[j]%mod;
				st[k]=(st[i]+sp[i]*prime[j]%mod+sp[k])%mod;
				sd[k]=(sd[i]*inv(st[i])%mod*st[k])%mod;
				break;
			}
			else
			{
				mu[k]=mu[i]*mu[prime[j]];
				p[i*prime[j]]=prime[j];
				t[i*prime[j]]=1+prime[j];
				d[i*prime[j]]=d[i]*d[prime[j]]%mod;
				sp[i*prime[j]]=prime[j]*prime[j]%mod;
				st[i*prime[j]]=(1+prime[j]+prime[j]*prime[j]%mod)%mod;
				sd[i*prime[j]]=sd[i]*sd[prime[j]]%mod;
			}
		}
	}
	for(LL i=1;i<=n;i++)
	s[i]=(d[i]+s[i-1])%mod;
	for(LL i=1;i<=n;i++)
	sd[i]=(sd[i-1]+i*sd[i]%mod)%mod;
	for(LL i=1;i<=n;i++)
	g[i]=1ll*i*d[i]%mod*s[i]%mod;
	for(LL d=1;d<=n;d++)
	for(LL T=d;T<=n;T+=d)
	f[T]=(f[T]+mu[d]*d%mod*d%mod*g[T/d]%mod+mod)%mod;
	for(LL i=1;i<=n;i++)
	f[i]=(f[i-1]+f[i])%mod;
}
LL calc(LL n)
{
	return (2ll*f[n]%mod-sd[n]+mod)%mod;
}
int main()
{
	freopen("a.in","r",stdin);
	freopen("a.out","w",stdout);
	init(1e6);
	LL T;
	cin>>T;
	int ca=0;
	while(T--)
	{
		LL n;
		ca++;
		scanf("%lld",&n);
		printf("Case #%d: %lld\n",ca,calc(n));
	}
	return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值