【BZOJ4916】神犇和蒟蒻(杜教筛水题)

点此看题面

大致题意:(sum_{i=1}^nmu(i^2))(sum_{i=1}^nphi(i^2))

前言

今天中午回班稍微早了一点,还没到吃饭时间,反正没事干就开始推式子。

结果突然就推出来了。。。?

然后吃午饭的时候又无聊地把杜教筛的式子凭借残余的印象在脑海里推了一遍。

于是发现这道题好像还真挺水的。

(sum_{i=1}^nmu(i^2))

。。。。。。

你确定这个式子放在这里不是用来搞笑的吗。。。

众所周知,对于一个数(n),只要它含平方因子,(mu(n))就等于(0),而(mu(i^2))的话。。。

显然,这个式子值始终为(1)(因为只有(mu(1)=1)啊!),甚至"良心"的样例中已经给出了答案。

(sum_{i=1}^nphi(i^2))

根据(phi)的定义,我们有:

[phi(i^2)=phi(i)i ]

考虑(n)这么大,肯定要用杜教筛。因此,我们需要找到一个函数(g(x))和一个函数(h(x)),使得:

[sum_{d|n}phi(d)dcdot g(frac nd)=h(n) ]

然后我们注意到式子中的(d)(frac nd),不难发现如果(g(x))(id(x))函数,(d)就刚好可以约掉!

[sum_{d|n}phi(d)dcdot id(frac nd)=nsum_{d|n}phi(d)=n^2 ]

因此,(g(x))(id(x)),而(h(x))(id^2(x))

于是我们就能愉快地套上杜教筛的式子了:(式子的推导可以看一看这篇博客:杜教筛入门

[S(n)=frac{sum_{i=1}^nh(i)-sum_{d=2}^ng(d)·S(lfloorfrac nd floor)}{g(1)} ]

[S(n)=frac{n(n+1)(2n+1)}6-sum_{d=2}^ndcdot S(lfloorfrac nd floor) ]

然后只要按照杜教筛的套路,除法分块+记忆化搞一下就行了。

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define SN 32000
#define X 1000000007
using namespace std;
int n,sn;
I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
class DuSieve
{
	private:
		#define SZ 1000000
		int Pt,P[SZ+5],phi[SZ+5],s[SZ+5],f1[SN+5],f2[SN+5];
	public:
		I DuSieve()//线性筛预处理
		{
			phi[1]=1;for(RI i=2,j,t;i<=SZ;++i)
				for(!P[i]&&(phi[P[++Pt]=i]=i-1),j=1;j<=Pt&&(t=i*P[j])<=SZ;++j)
					if(P[t]=1,i%P[j]) phi[t]=phi[i]*(P[j]-1);else {phi[t]=phi[i]*P[j];break;}
			for(RI i=1;i<=SZ;++i) s[i]=(1LL*phi[i]*i+s[i-1])%X;//初始化小数据
		}
		I int SPhi(CI x)//杜教筛
		{
			if(x<=SZ) return s[x];//如果较小,直接返回预处理的答案
			if(x<=sn) {if(f1[x]) return f1[x];}else {if(f2[n/x]) return f2[n/x];}//用两个数表记忆化
			RI l,r,t=0;for(l=2;l<=x;l=r+1)//除法分块
				r=x/(x/l),t=((1LL*(l+r)*(r-l+1)/2)%X*SPhi(x/l)+t)%X;
			return (x<=sn?f1[x]:f2[n/x])=(1LL*x*(x+1)%X*(2*x+1)%X*166666668-t+X)%X;//计算最终结果,并记忆化存储
		}
}D;
int main()
{
	return scanf("%d",&n),sn=sqrt(n),printf("1
%d
",D.SPhi(n)),0;
}
原文地址:https://www.cnblogs.com/chenxiaoran666/p/BZOJ4916.html