【学术篇】分析矿洞 杜教筛

数论什么的都去死吧！

看着题解我都能化式子用完4页草纸。。。
另外吐槽一句出题人的拼音学的是真好, 不知道是不是故意的.
其实题解已经写得挺详细的了.
我就是提一些出题人觉得太easy没必要提但是做题还是需要的一些东西....(因为这些东西我基本都是现学的)

然而之前刚刚学完mobius反演就暂时性脱坑的我啥也不会啊..
看到前排dp和曲神在水luogu的欢(bao/du)乐(ling/liu)赛, 就想去看看.
然后就点了报名但是发现自己什么都不会.

去看了看T1. 就是这道题.
然后成功的化出了第一步的式子.
这样就可以水30pts了.
一眼看出应该是反演类型的题目, 但是真的tmd不会啊,, 80pts都水不到.
(部分分给的也是有点迷, 80pts和100pts完全不是一样东西好么= =)

30pts:
通过一眼看出法可以得到激光扫到的第一个点的坐标是

[(frac x{gcd(x,y)},frac y{gcd(x,y)}) ]

所以

[(frac{x+y}{frac x{gcd(x,y)}+frac y{gcd(x,y)}})^2=gcd^2(x,y) ]

于是很显然地就是要求出

[sum_{i=1}^Nsum_{j=1}^Nvarphi(gcd^2(i,j))cdots① ]

这个东西, 于是就变成了一道纯数论题. (本来就是一道纯数论题不是?!)

然后就(O(n^2logn))乱搞就30pts到手了.

80pts:
继续化式子.
对于(varphi)我们有这么一种操作:

[varphi(n)=prod_ip_i^{k_i-1}(p_i-1) ]

所以可以得到

[varphi(n^m)=prod_ip_i^{mk_i-1}(p_i-1)=prod_ip_i^{k_i-1}(p_i-1)*p_i^{(m-1)k_i}=varphi(n)*n^{m-1} ]

我们就可以把①中的gcd拆出来, 把平方去掉, 变成

[①=sum_{i=1}^{N}sum_{j=1}^{N}varphi(gcd(i,j))*gcd(i,j) ]

然后反演的常见套路之枚举公因数

[=sum_{i=1}^{N}sum_{j=1}^{N}sum_{d|i,d|j}varphi(d)*d*b[gcd(i,j)=d] ]

其中(b[x])表示(x)的真假性, (x)真则(b[x]=1), 否则(b[x]=0)
也就等价于在一个边长为(left lfloorfrac nd ight floor)的方阵中找互质的((i,j)), 然后对(d*varphi(d))求和.

[=sum_{d=1}^Nvarphi(d)*dsum_{i=1}^{left lfloorfrac Nd ight floor}sum_{j=1}^{left lfloorfrac Nd ight floor} b[gcd(i,j)=1]cdots② ]

后半边式子有点眼熟?? 我们好像在哪里见过这种形式.

仪仗队!!!!!!

我们可以得出

[sum_{i=1}^Nsum_{j=1}^Nb[gcd(i,j)=1] ]

这个式子的结果是(2*sum_{i=1}^Nvarphi(N)-1).
所以

[②=sum_{d=1}^Nvarphi(d)*d*(2*sum_{i=1}^{left lfloorfrac Nd ight floor}varphi(i)-1)cdots③ ]

这样这个式子就化到头了.
而此时我们枚举(d)就可以做到(O(nsqrt n))求解了..
这样就能水到80pts. (个人感觉80pts部分分给得略高了)

100pts
满分做法就要用到一种高端科技了..

杜教筛!

顾名思义是一种筛法. 但是要比线筛快一些.
举个栗子, 我们来求一下

[sum_{i=1}^Nvarphi(i)$$ (其实跟上面的式子是有联系的OvO) 那么我们看数据范围想算法: - $N<=1000$? - 枚举, 每次从头扫求一遍欧拉函数都能过. - $N<=10000000$? - $varphi(x)$是个积性函数, 直接线筛一下就好咯. - $N<=10000000000$? - 这个...... $O(n)$过不了啊.. 这就是说我们必须要想别的方法了. 比如**杜教筛**.. 我们先来化一波式子, 尽量把N变小到能做的范围. 对于$varphi$函数有一条: $$sum_{d|n}varphi(d)=n]

那么

[sum_{d|n,d<n}varphi(d)+varphi(n)=n\ varphi(n)=n-sum_{d|n,d<n}varphi(d)]

我们$$phi(i)=
sum_{i=1}^{Nvarphi(i)=sum_{i=1}}N(i-sum_{d|i,d<i}varphi(d))=sum_{i=1}^Ni-sum_{i=2}Nsum_{d|i,d<i}varphi(d)
=sum_{i=1}^Ni-sum_{frac id=2}^Nsum_{d=1}{leftlfloorfrac n{frac id} ight floor}varphi(d)
令j=frac id,则phi(i)sum_{i=1}^{Nvarphi(i)=sum_{i=1}}Ni-sum_{j=2}^Nsum_{d=1}{leftlfloorfrac nj ight floor}varphi(d)=sum_{i=1}^N-sum_{j=2}Nphi(leftlfloorfrac nj ight floor)$$
其中减号前面的显然是可以(O(1))计算的(别说你不会), 后面的值是不会超过(sqrt n)个的, 我们枚举因数递归计算即可.
代码太长而且基本是这题代码的子集就不糊在这里了..留个传送门自己去看吧.

然后我们回到这道题. 我们化出了③式, 为了防止忘掉, 我们再贴一遍.

[sum_{d=1}^Nvarphi(d)*d*(2*sum_{i=1}^{left lfloorfrac Nd ight floor}varphi(i)-1) ]

这样括号里面的刚刚学习了怎么筛(所以说是子集嘛), 所以问题就集中在了前面的

[sum_{d=1}^Nvarphi(d)*d ]

怎么快速的筛出来. 而这个题解已经说的挺清楚了的(说你是不是懒得继续化了←_←)
我们令(f(i)=sum_{d=1}^Nvarphi(d)*d)
这个(varphi(d)*d=varphi(d^2)), 我们就猜测和(sum_{i=1}^Ni^2)(这个式子可以用平方和公式(O(1))求哟~)有什么联♂系.
(这个理由是蒙的, 比赛的时候怎么凑知道该怎么凑出来 不是很急但是会在线等= =)
那我们就化一下

[sum_{i=1}^Ni^2=sum_{i=1}^N(i*sum_{d|i}varphi(d))=sum_{i=1}^N(frac id*sum_{d|n}varphi(d)*d) ]

这样凑出了(varphi(d)*d)的形式. 但是还没有(sum_{i=1}^n)的形式, 考虑枚举.
我们令(t=frac id), 然后枚举(t). 因为要求1..N的和, 所以如果有(t*d>N)的(d)显然不能对答案做出贡献, 所以我们枚举(d)的时候枚举到(leftlfloorfrac Nt ight floor)即可. 也就是说

[sum_{i=1}^Ni^2=sum_{t=1}^Nt*sum_{d=1}^{leftlfloorfrac Nt ight floor}varphi(d)*d=sum_{d=1}^Nvarphi(d)*d+sum_{t=2}^Nt*sum_{d=1}^{leftlfloorfrac Nt ight floor}varphi(d)*d ]

那我们就可以得到

[f(i)=sum_{i=1}^Ni^2-sum_{t=2}^Nt*sum_{d=1}^{leftlfloorfrac Nt ight floor}varphi(d)*d=sum_{i=1}^Ni^2-sum_{t=2}^Nt*f(leftlfloorfrac Nt ight floor) ]

就这么得到了一个杜教筛的形式, 就可以仿照上面做咯~

md化式子化到吐系列……

听说这题卡常数? 但也没怎么卡嘛 感觉随便一跑就轻松第一页了?
(好像出题人改了题, 所以只算改题之后的话应该就rank1了= =之前q1~q4数组开大了memset废掉好多时间= =
用了一些小trick 比如全程能开int绝对不开long long... (所以因为少%了一次第一次交80pts嘤嘤嘤...
比如把结果记忆化一下. 看到std这一步是用map做的.
但是因为都是(N)的因数, 我们就可以把这些因数分比一个数大的和比这个数小的两类分
这个数大约取(sqrt N)即可. (大点也能存下, 但是小点好像可能会出现冲突)

const SQ=1e6/7; 	//142857
int p1[SQ],p2[SQ];
int& getaddr(LL x){
	if(x<SQ) return p1[x];
	return p2[n/x];
} //返回储存地址

这样就能少个log, 就会非常快...
开int的话一定要记得:步步取模, 强转long long!
Emmmm还有一种非常无良的针对数据的压常trick就是n<=1e5的时候线筛可以少筛一点...(这样就稳稳rank1了)

代码(为什么不去看std呢)

#include <cstdio>
#include <cstring>
const int N=5e6+5;
const int P=1e9+7;
const int SQ=1e6/7;
typedef long long LL;
LL n;
int p1[SQ],p2[SQ],p3[SQ],p4[SQ];
int prime[N],euler[N],mu[N],eusum[N],eumul[N],tot;
bool notp[N];
void shai(){
	notp[1]=euler[1]=mu[1]=eusum[1]=eumul[1]=1;
	for(int i=2;i<=5e6;++i){
		if(!notp[i]){
			prime[++tot]=i;
			euler[i]=i-1;
			mu[i]=-1;
		}
		for(int j=1;j<=tot&&i*prime[j]<5e6;++j){
			int k=i*prime[j];
			notp[k]=1;
			if(i%prime[j]==0){
				mu[k]=0;
				euler[k]=euler[i]*prime[j];
				break;
			}
			else{
				mu[k]=-mu[i];
				euler[k]=euler[i]*(prime[j]-1);
			}
		}
		eusum[i]=(eusum[i-1]+euler[i])%P;
		eumul[i]=(eumul[i-1]+1LL*euler[i]*i%P)%P;
	}
}
inline int qpow(int a,int b,int s=1){
	for(;b;b>>=1,a=1LL*a*a%P)
		if(b&1) s=1LL*s*a%P;
	return s;
}
int inv2=qpow(2,P-2),inv6=qpow(6,P-2);
inline int& getaddr(LL x,bool flag){
	if(flag){
		if(x<SQ) return p1[x];
		return p2[n/x];
	}
	if(x<SQ) return p3[x];
	return p4[n/x];	
}
int eulersum(LL x){
	if(x<=5e6) return eusum[x];
	int& addr=getaddr(x,1);
	if(addr!=-1) return addr;
	int ans; LL last;
	ans=1LL*(x%P)*(x%P+1)%P*inv2%P;
	for(LL i=2;i<=x;i=last+1){
		last=x/(x/i);
		ans=(ans-1LL*(last-i+1)*eulersum(x/i)%P)%P;
	}
	return addr=(ans%P+P)%P;
}
int eumulsum(LL x){
	if(x<=5e6) return eumul[x];
	int& addr=getaddr(x,0);
	if(addr!=-1) return addr;
	int ans; LL last;
	ans=1LL*(x%P)*(x%P+1)%P*((x+x+1)%P)%P*inv6%P;
	for(LL i=2;i<=x;i=last+1){
		last=x/(x/i);
		ans=(ans-1LL*(i+last)%P*(last-i+1)%P*inv2%P*eumulsum(x/i)%P)%P;
	}
	return addr=(ans%P+P)%P;
}
int main(){
	memset(p1,-1,sizeof(p1));
	memset(p2,-1,sizeof(p2));
	memset(p3,-1,sizeof(p3));
	memset(p4,-1,sizeof(p4));
	shai(); scanf("%lld",&n);
	int s=0; LL last;
	for(LL i=1;i<=n;i=last+1){
		last=n/(n/i);
		s=(s+1LL*(eumulsum(last)-eumulsum(i-1)+P)*((2*eulersum(n/i)%P-1+P)%P)%P)%P;
	}
	printf("%d",s);
}

真是要吐了 完结撒花~