真的就是普普通通的同余方程 [doge]
# 题面
多组询问，每组询问给定正整数 (P,x)，其中 (P) 为偶数且不包含平方因子。对每组询问，求下列方程在 (ain[0,P),bin[0,P)) 有多少个解 ((a,b))：
数据规模：不超过 (10^5) 次询问，(Ple 10^7)。
# 解析
你可以不知道为什么，但是起码你要会打表嘛
——Tiw_Air_OAO
首先「(P) 没有平方因子」这个条件非常棒棒，对 CRT 非常友好。我们可以把 (P) 分解质因子，对每个质因子 (p_i) 求解 (a^2+b^2equiv xpmod{p_i})。
假如现在我们已经对每个质因子求出了答案，怎么合并得到全局的答案？其实非常简单，所有的答案都乘起来就好了。为什么呢？
不妨从下面这个角度理解。我们来举一个栗子 (P=42=2 imes3 imes7)，看看 (P) 的剩余系如何构成——
根据 CRT，如下的线性同余方程组有唯一解：
这意味着我们可以用三元组 ((x_0,x_1,x_2)) 与 (P=42) 的剩余系内的任意元素一一对应。
再回到本问题，假如我们对 (p_i) 求得了解 ((a_i,b_i))，根据 CRT，我们可以列出同余方程组：
那么 ((a_0,a_1,cdots,a_t)) 可以与模 (P) 意义下的 (a) 一一对应，(b) 同理。于是可以直接乘起来。
于是现在就把问题转化成了「(P) 为一个奇素数」的问题。
为什么让你打表呢？你会惊奇地发现，对于同一个 (P)，(x
eq 0) 的答案都是一样的！哇好神奇！但是既然写题解，肯定不能告诉你直接打表就能发现……
不妨暴力地枚举 (iequiv a^2)，则 (b^2equiv x-i)，然后我们求出 (i) 和 (x-i) 的二次剩余数量，就可以求得 (a^2equiv i) 时的解 ((a,b)) 数量了。
怎么统计二次剩余的数量呢？让我们又来看看这个老朋友——勒让德符号，根据其实际意义不难得到 (left(frac{a}{p} ight)+1) 就是 (x^2equiv apmod p) 的解的数量。
于是写出答案的式子：
由 (ig(frac{0}{p}ig)=0)，可以改变 (i) 的枚举起点到 (1)：
由勒让德符号的计算式 (ig(frac{x}{p}ig)=x^{frac{p-1}{2}})，不难得到 (ig(frac{x}{p}ig)ig(frac{y}{p}ig)=ig(frac{xy}{p}ig))。又有二次剩余的数量和非二次剩余的数量相等，所以有 (sumlimits_{i=1}^{p-1}ig(frac{i}{p}ig)=0)。上式可以继续简化：
根据下面的推导：
则：
最后我们就可以知道为什么 (x=0) 的答案就和其他 (x eq0) 的答案不一样——
-
当 (x=0) 时，答案为 ((p-1)ig( frac{-1}{p}ig)+p=p+(p-1)(-1)^{frac{p-1}{2}})；
-
当 (x eq0) 时，由 (i) 从 (1) 枚举至 (p-1)，则 (i^{-1}) 取值也取尽 (1sim p-1)，(xi^{-1}) 同理。
那么 (xi^{-1}-1) 取尽 (0sim p-2)。答案为
[sum_{i=0}^{p-2}ig( frac{i}{p}ig)+p=sum_{i=0}^{p-1}ig( frac{i}{p}ig)-ig( frac{p-1}{p}ig)+p=p-(-1)^{frac{p-1}{2}} ]
# 源代码
/*Lucky_Glass*/
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
inline int rin(int &r){
int b=1,c=getchar();r=0;
while(c<'0' || '9'<c) b=c=='-'?-1:b,c=getchar();
while('0'<=c && c<='9') r=(r<<1)+(r<<3)+(c^'0'),c=getchar();
return r*=b;
}
#define con(type) const type &
const int N=1e7+10;
int prm[N/5],mndv[N],nprm;
int solve(con(int)varx,con(int)varp){
if(varx) return varp-(((varp-1)>>1&1)?-1:1);
else return varp+(((varp-1)>>1&1)?-1:1)*(varp-1);
}
void init(){
for(int i=2;i<N;i++){
if(!mndv[i]) prm[++nprm]=mndv[i]=i;
for(int j=1;j<=nprm && prm[j]*i<N;j++){
mndv[prm[j]*i]=prm[j];
if(i%prm[j]==0) break;
}
}
}
int main(){
init();
int varx,varp,ncas;
rin(ncas);
while(ncas--){
rin(varp),rin(varx);
long long ans=1;
while(varp>1){
ans*=solve(varx%mndv[varp],mndv[varp]);
varp/=mndv[varp];
}
printf("%lld
",ans);
}
return 0;
}
THE END
Thanks for reading!
你是春山 你是岁酒
你是自由 你是误谬
你是颠沛流离之后 我绮丽的愁
你是沙鸥 你是滴漏
你是白昼 你是不朽
——《你是我遥不可及的梦》By 苍穹/papaw泡泡
> Link 你是我遥不可及的梦-Bilibili