LOJ.6053.简单的函数(Min_25筛)

题目链接
Min_25筛见这里：
https://www.cnblogs.com/cjyyb/p/9185093.html
https://www.cnblogs.com/zhoushuyu/p/9187319.html
https://www.cnblogs.com/SovietPower/p/10101811.html

---

(Description)

给定(n)，求积性函数(f(p^c)=poplus c)的前缀和。(oplus)表示异或运算。
(nleq 10^{10})。

(Solution)

所求积性函数为(f(p^c)=poplus c,quad pin Prime)。
先考虑质数的贡献。因为除(2)以外的质数(p)都是奇数，所以(f(p)=p-1)，(f(2))就是(2+1=3)了。

不妨先把(f(2))也看做(f(2)=p-1)。
这样(f(p)=p-1)，但还不是积性函数，但是我们可以拆成两个积性函数的和：(f(p)=g(p)-h(p))，其中(g(p)=p, h(p)=[p是质数])。分别计算这两个函数的前缀和。

然后，首先计算初值(g(n,0))，把所有合数看做质数，那么(g(n),h(n))的前缀和也就是初值分别是(frac{n(n+1)}{2}-1)，(n-1)（不考虑(f(1))）。

然后计算(g(x,|P|))。在外层枚举(j)把第二维滚动掉。
另外每次从(frac{n}{P_j})转移，都是整除，结果只有(O(sqrt(n)))种（好像是(2sqrt(n))？不知道为什么...），所以可以先离散化。
然后套式子得到(g(x,|P|))。

然后套式子计算(S(n,1))。当(j=1)时给结果加个(2)，因为(f(P_1)=f(2))是按(2-1=1)计算的，应该是(3)。
这里直接递归算就行了，而且不需要记忆化。

//1214ms	4.12M
#include <cmath>
#include <cstdio>
#include <algorithm>
#define mod 1000000007
#define Mod(x) x>=mod&&(x-=mod)
typedef long long LL;
const int N=2e5+5;

int Sqr,cnt,P[N>>2],sp[N],id1[N],id2[N],g[N],h[N];
LL n,w[N];
bool notP[N];

void Init(int n)
{
	notP[1]=1;
	for(int i=2; i<=n; ++i)
	{
		if(!notP[i]) P[++cnt]=i, sp[cnt]=sp[cnt-1]+i, Mod(sp[cnt]);
		for(int j=1; j<=cnt && i*P[j]<=n; ++j)
		{
			notP[i*P[j]]=1;
			if(!(i%P[j])) break;
		}
	}
}
int S(LL x,int y)
{
	if(x<=1||(y!=cnt+1&&P[y]>x)) return 0;//注意y=cnt+1时也需要计算！
	int k=x<=Sqr?id1[x]:id2[n/x];
	LL res=g[k]-sp[y-1]-h[k]+y-1;//g-h
	if(y==1) res+=2;//f(2)还是就放到里面算吧 否则要判下n<2。
	for(int i=y; /*i<=cnt&&*/1ll*P[i]*P[i]<=x; ++i)
	{
		LL p=P[i],p1=p,p2=p*p;
		for(int e=1; p2<=x; ++e,p1=p2,p2*=p)
			res+=1ll*(p^e)*S(x/p1,i+1)%mod+(p^(e+1));
		res%=mod;
	}
	return res%mod;
}

 main()
{
	scanf("%lld",&n); Init(Sqr=sqrt(n+0.5));
	int m=0;
	for(LL i=1,j; i<=n; i=j+1)
	{
		w[++m]=n/i, j=n/w[m];
		if(w[m]<=Sqr) id1[w[m]]=m;
		else id2[j]=m;
		g[m]=w[m]&1?w[m]%mod*((w[m]+1>>1)%mod)%mod-1:(w[m]>>1)%mod*(w[m]%mod+1)%mod-1;//w乘之前要取模！
		h[m]=(w[m]-1)%mod;
	}
	P[cnt+1]=1e9, w[m+1]=-1;
	for(int j=1; j<=cnt; ++j)
	{
		int pj=P[j]; LL lim=1ll*pj*pj;
		for(int i=1; lim<=w[i]; ++i)
		{
			int k=w[i]/pj<=Sqr?id1[w[i]/pj]:id2[n/(w[i]/pj)];//n/(w/pj)! id1[x]=x，但id2[x]的编号是对[n/x]的。
			(g[i]-=1ll*pj*(g[k]-sp[j-1])%mod)%=mod;//g[k]-sp[j-1]有可能是负的，如果+mod会爆int！
			h[i]+=mod-h[k]+j-1, Mod(h[i]);
		}
	}
	printf("%d
",((S(n,1)+1)%mod+mod)%mod);

	return 0;
}