SPOJ DIVCNT2 - Counting Divisors (square)

版权声明:本文为博主原创文章,未经博主允许必须转载。 https://blog.csdn.net/C20181220_xiang_m_y/article/details/85036860

题目描述:

i = 1 n σ 0 ( i 2 ) \sum_{i=1}^n\sigma_0(i^2)
其中 σ 0 \sigma_0 表示约数个数, n 1 0 12 n\le10^{12} ,时间限制20s

题目分析:

σ 0 ( i ) = d i σ 0 ( i 2 ) = d i 2 ω ( d ) \large\sigma_0(i)=\sum_{d|i}\\\sigma_0(i^2)=\sum_{d|i}2^{\omega(d)}
其中, ω ( d ) \omega(d) 表示 d d 的质因子个数,上式可以理解为从i中选出一些质因子,每个质因子可以取当前次幂或者加上最高次幂两种选择。

可看出 2 ω ( d ) 2^{\omega(d)} 为d的无平方因子个数,所以
2 ω ( d ) = k d μ ( k ) \large2^{\omega(d)}=\sum_{k|d}|\mu(k)|
整理一下,记 g ( i ) = μ ( k ) g(i)=|\mu(k)| ,那么 σ 0 ( i 2 ) = ( g I ) I = g ( I I ) = g σ 0 \sigma_0(i^2)=(g*I)*I=g*(I*I)=g*\sigma_0
所以 i = 1 n σ 0 ( i 2 ) = i = 1 n d i μ ( i d ) σ 0 ( d ) = i = 1 n μ ( i ) d = 1 n i σ 0 ( d ) \sum_{i=1}^n\sigma_0(i^2)=\sum_{i=1}^n\sum_{d|i}|\mu(\frac id)|*\sigma_0(d) \\=\sum_{i=1}^n|\mu(i)|\sum_{d=1}^{\frac ni}\sigma_0(d)
显然需要处理 μ ( i ) |\mu(i)| σ 0 ( i ) \sigma_0(i) 的前缀和
线性筛五千万到一亿(少了会TLE),剩下的可以这么算:

  • 从" n n 以内无平方因子数的个数"的意义出发,考虑容斥,可以得出:
    i = 1 n μ ( i ) = i = 1 n μ ( i ) n i 2 \sum_{i=1}^n|\mu(i)|=\sum_{i=1}^{\sqrt n}\mu(i)*\left\lfloor\frac n{i^2}\right\rfloor
    O ( n ) O(\sqrt n) 计算
  • 约数个数的前缀和很好算:
    i = 1 n σ 0 ( i ) = i = 1 n n i \sum_{i=1}^n\sigma_0(i)=\sum_{i=1}^n\left\lfloor\frac ni\right\rfloor
    O ( n ) O(\sqrt n) 分块计算

与杜教筛的时间复杂度分析类似,总时间 O ( n 2 3 ) \large O(n^{\frac 23})

#include<cstdio>
#include<algorithm>
#define LL long long
using namespace std;
const int N = 50000000;
int p[N/10],mu[N+5],sd[N+5],sm[N+5];
bool v[N+5];
void Prime(int N)
{
	mu[1]=sd[1]=sm[1]=1;int cnt=0;
	for(int i=2;i<=N;i++)
	{
		if(!v[i]) p[++cnt]=i,mu[i]=-1,sm[i]=sd[i]=2;//sm now replace ak to use
		for(int j=1,k;j<=cnt&&p[j]*i<=N;j++)
		{
			v[k=p[j]*i]=1;
			if(i%p[j]==0) {mu[k]=0,sd[k]=sd[i]/sm[i]*(sm[i]+1),sm[k]=sm[i]+1;break;}
			mu[k]=-mu[i],sd[k]=sd[i]*sd[p[j]],sm[k]=2;
		}
	}
	for(int i=1;i<=N;i++) sd[i]+=sd[i-1],sm[i]=sm[i-1]+abs(mu[i]);
}
LL Summu(LL n)
{
	if(n<=N) return sm[n];
	LL ret=0;
	for(LL i=1;i*i<=n;i++) ret+=mu[i]*n/(i*i);
	return ret;
}
LL Sumd(LL n)
{
	if(n<=N) return sd[n];
	LL ret=0;
	for(LL i=1,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		ret+=(j-i+1)*(n/i);
	}
	return ret;
}
LL solve(LL n)
{
	LL ans=0;
	for(LL i=1,j,pre=0,tmp;i<=n;i=j+1)
	{
		j=n/(n/i);
		tmp=Summu(j);
		ans=ans+Sumd(n/i)*(tmp-pre);
		pre=tmp;
	}
	return ans;
}
int main()
{
	LL n[10005];int T;
	scanf("%d",&T);
	for(int i=1;i<=T;i++) scanf("%lld",&n[i]);
	if(*max_element(n+1,n+1+T)<=10000) Prime(10000);
	else Prime(N);
	for(int i=1;i<=T;i++) printf("%lld\n",solve(n[i]));
}

猜你喜欢

转载自blog.csdn.net/C20181220_xiang_m_y/article/details/85036860