题目大意
求
(
∑
d
=
L
R
(
f
d
xor
K
)
)
(
m
o
d
P
)
\Big(\sum\limits_{d=L}^{R} (f_d\text{ xor } K)\Big)\pmod{P}
(d=L∑R(fd xor K))(modP)
其中
f
d
f_d
fd表示在空间直角坐标系下,圆心
(
0
,
0
,
0
)
(0,0,0)
(0,0,0)且半径为
d
d
d的球上的整点个数
多组数据(
T
≤
10
T\le 10
T≤10),其中
0
≤
L
≤
R
≤
1
0
13
,
0
≤
K
≤
1
0
18
,
R
−
L
+
1
≤
1
0
6
0\le L\le R\le 10^{13}, 0\le K\le 10^{18},R-L+1\le 10^6
0≤L≤R≤1013,0≤K≤1018,R−L+1≤106,质数P满足
P
≤
3
×
1
0
13
P\le 3\times 10^{13}
P≤3×1013
分析
通过打表/OEIS发现
f
d
6
\frac {f_d} 6
6fd是积性函数
且当p为质数且
p
=
4
k
+
1
p=4k+1
p=4k+1时,
f
(
p
e
)
=
p
e
f(p^e)=p^e
f(pe)=pe;
当p为质数且
p
=
4
k
+
3
p=4k+3
p=4k+3时,
f
(
p
e
)
=
p
e
+
2
(
p
e
−
1
)
p
−
1
f(p^e)=p^e+\frac{2(p^e-1)}{p-1}
f(pe)=pe+p−12(pe−1)
注意到区间长度 n = R − L + 1 ≤ 1 0 6 n=R-L+1\le 10^6 n=R−L+1≤106,因此通过预处理 R \sqrt R R以内的素数,之后通过区间筛的方式求出 f x ( L ≤ x ≤ R ) f_x(L\le x \le R) fx(L≤x≤R),然后就可以 O ( n ) O(n) O(n)求上面那个式子了
总复杂度 O ( T n log n ) O(Tn\log n) O(Tnlogn)
代码
#include<bits/stdc++.h>
#define LL long long
using namespace std;
LL L,R,k,p;
int T;
const int lim=4000000;
int prime[lim/10+5];
LL f[lim+5],id[lim+5];
bool vis[lim+5];
void init()
{
for (int i=2;i<=lim;++i)
{
if (!vis[i]) prime[++prime[0]]=i;
for (int j=1;j<=prime[0]&&prime[j]*i<=lim;++j)
{
vis[prime[j]*i]=1;
if (i%prime[j]==0) break;
}
}
}
void get(LL L,LL R)//区间筛[L,R]以内的f
{
for (int i=1;i<=R-L+1;++i) id[i]=i+L-1,f[i]=1;
for (int i=1;i<=prime[0]&&prime[i]<=R;++i)
{
LL last=(L/prime[i])*prime[i],mk;
if (last<L) last+=prime[i];
for (;last<=R;last+=prime[i])
{
mk=1;
while (id[last-L+1]%prime[i]==0)
id[last-L+1]/=prime[i],
mk*=prime[i];
if (prime[i]==2) mk=1;
else if (prime[i]%4==3) mk+=(2*(mk-1)/(prime[i]-1));
f[last-L+1]*=mk;
}
}
for (int i=1;i<=R-L+1;++i)
{
if (id[i]>1)
if (id[i]%4==3)
f[i]*=(id[i]+2);
else
f[i]*=id[i];
f[i]*=6;
f[i]^=k;
f[i]%=p;
}
}
main()
{
init();
for(scanf("%d",&T);T;--T)
{
scanf("%lld%lld%lld%lld",&L,&R,&k,&p);
LL ans=0;
if (L==0) ++L,(ans+=k)%=p;
get(L,R);
for (int i=1;i<=R-L+1;++i) ans=(ans+f[i])%p;
printf("%lld\n",ans);
}
}