「SOL」同余方程(LOJ)

真的就是普普通通的同余方程 [doge]

---

# 题面

多组询问，每组询问给定正整数 (P,x)，其中 (P) 为偶数且不包含平方因子。对每组询问，求下列方程在 (ain[0,P),bin[0,P)) 有多少个解 ((a,b))：

[a^2+b^2equiv xpmod P ]

数据规模：不超过 (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，如下的线性同余方程组有唯一解：

[egin{cases} xequiv x_0pmod 2\ xequiv x_1pmod 3\ xequiv x_2pmod 7 end{cases} ]

这意味着我们可以用三元组 ((x_0,x_1,x_2)) 与 (P=42) 的剩余系内的任意元素一一对应。

再回到本问题，假如我们对 (p_i) 求得了解 ((a_i,b_i))，根据 CRT，我们可以列出同余方程组：

[egin{cases} aequiv a_0pmod{p_0}\ aequiv a_1pmod{p_1}\ cdots\ aequiv a_tpmod{p_t} end{cases} egin{cases} bequiv b_0pmod{p_0}\ bequiv b_1pmod{p_1}\ cdots\ bequiv b_tpmod{p_t} end{cases} ]

那么 ((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) 的解的数量。

于是写出答案的式子：

[egin{aligned} &sum_{i=0}^{p-1}Big[ig( frac{i}{p}ig)+1Big]Big[ig( frac{x-i}{p}ig)+1Big]\ =&sum_{i=0}^{p-1}ig( frac{i}{p}ig)ig( frac{x-i}{p}ig)+2sum_{i=0}^{p-1}ig( frac{i}{p}ig)+p end{aligned} ]

由 (ig(frac{0}{p}ig)=0)，可以改变 (i) 的枚举起点到 (1)：

[sum_{i=1}^{p-1}ig( frac{i}{p}ig)ig( frac{x-i}{p}ig)+2sum_{i=1}^{p-1}ig( frac{i}{p}ig)+p ]

由勒让德符号的计算式 (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)。上式可以继续简化：

[sum_{i=1}^{p-1}ig( frac{i(x-i)}{p}ig)+p ]

根据下面的推导：

[ig( frac{xi^{-1}-1}{p}ig)ig( frac{i^2}{p}ig)=ig( frac{i(x-i)}{p}ig)=ig( frac{xi^{-1}-1}{p}ig) ]

则：

[sum_{i=1}^{p-1}ig( frac{xi^{-1}-1}{p}ig)+p ]

最后我们就可以知道为什么 (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