題目鏈接

\(Description\)
令\(f(d)\)表示空間中到原點距離為\(d\)的整點個數，給定\(L,R,k,p\)，求

\(L,R\leq 10^{13},R-L+1\leq 10^6\)。

\(Solution\)
必然是找規律。OEIS可以直接找到規律，或者打表。

\(f\)都是\(6\)的倍數，令\(g(x)=\frac{f(x)}{6}\)，可以發現\(g(x)\)為積性函數，即\(g(\prod p_i^{k^i})=\prod g(p_i^{k_i})\)。
列出\(g(p)\ （p為質數）\)可以知道，若\(p\%4=1\)，\(g(p)=p\)；若\(p\%4=3\)，\(g(p)=p+2\)。
再看\(g(p^2)\ （p為質數）\)可以知道，若\(p\%4=1\)，\(g(p^2)=p^2\)；若\(p\%4=3\)，\(g(p^2)=p^2+h(p)\)，其中\(h(p)=8k\)，\(k\)是\(p\)為第幾個\(\%4=3\)的質數。
考慮如何更好地表示\(h(p)\)。
由OEIS結論可發現\((p-1)h(p)=2(p^2-1)\)，即\(g(p^2)=p^2+\frac{2(p^2-1)}{p-1}\)，且該結論對\(p^k\)同樣成立。
所以

求出\(\sqrt n\)以內的質數，對每個數就可以\(O(\log n)\)分解質因數求了。
復雜度\(O(Tn\log n)\)。

PS：注意&的運算優先級比==低。。&運算要加括號！
\(ans\)求和過程中最多是\(10^6\)個\(2\times 10^{13}\)的數，所以還是要取模！

//1261ms	20MB
#include <bits/stdc++.h>
#define pc putchar
#define gc() getchar()
#define pb emplace_back
typedef long long LL;
const int N=3162300,M=1e6+5;
const LL LIM=2e18;

int cnt,P[230000];
bool not_P[N];

inline LL read()
{
	LL now=0,f=1; char c=gc();
	for(;!isdigit(c);c=='-'&&(f=-1),c=gc());
	for(;isdigit(c);now=now*10+c-48,c=gc());
	return now*f;
}
void Init(const int n)
{
	for(int i=2; i<=n; ++i)
	{
		if(!not_P[i]) P[++cnt]=i;
		for(int j=1; j<=cnt && i*P[j]<=n; ++j)
		{
			not_P[i*P[j]]=1;
			if(!(i%P[j])) break;
		}
	}
}
void Solve()
{
	static LL A[M],f[M];

	LL L=read(),R=read(),K=read(),mod=read(),ans=0;
	if(!L) ans+=K, ++L;
	int n=(int)(R-L+1);
	for(int i=1; i<=n; ++i) A[i]=i+L-1,f[i]=1;

	for(int i=1,lim=::cnt; i<=lim&&P[i]<=R; ++i)
	{
		int p=P[i];
		LL now=L/p*p; now<L&&(now+=p);
		for(now-=L-1; now<=n; now+=p)
		{
			LL pk=1;
			while(!(A[now]%p)) A[now]/=p, pk*=p;
			if(p==2) pk=1;
			else if((p&3)==3) pk+=2*(pk-1)/(p-1);
			f[now]*=pk;
		}
	}
	for(int i=1; i<=n; ++i)
	{
		if((A[i]&3)==3) f[i]*=(A[i]+2);//本身/剩下為大質數 
		else f[i]*=A[i];
		ans+=((f[i]*6)^K)%mod;
		ans>=LIM&&(ans%=mod);//取模！ 
	}
	printf("%lld\n",ans%mod);
}

int main()
{
	Init(N-2);
	for(int T=read(); T--; Solve());

	return 0;
}