「SOL」同余方程(LOJ)

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

---

# 题面

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

\[a^2+b^2\equiv x\pmod P \]

数据规模：不超过 \(10^5\) 次询问，\(P\le 10^7\)。

---

# 解析

你可以不知道为什么，但是起码你要会打表嘛
——Tiw_Air_OAO

首先「\(P\) 没有平方因子」这个条件非常棒棒，对 CRT 非常友好。我们可以把 \(P\) 分解质因子，对每个质因子 \(p_i\) 求解 \(a^2+b^2\equiv x\pmod{p_i}\)。

假如现在我们已经对每个质因子求出了答案，怎么合并得到全局的答案？其实非常简单，所有的答案都乘起来就好了。为什么呢？

不妨从下面这个角度理解。我们来举一个栗子 \(P=42=2\times3\times7\)，看看 \(P\) 的剩余系如何构成——

根据 CRT，如下的线性同余方程组有唯一解：

\[\begin{cases} x\equiv x_0\pmod 2\\ x\equiv x_1\pmod 3\\ x\equiv x_2\pmod 7 \end{cases} \]

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

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

\[\begin{cases} a\equiv a_0\pmod{p_0}\\ a\equiv a_1\pmod{p_1}\\ \cdots\\ a\equiv a_t\pmod{p_t} \end{cases} \begin{cases} b\equiv b_0\pmod{p_0}\\ b\equiv b_1\pmod{p_1}\\ \cdots\\ b\equiv b_t\pmod{p_t} \end{cases} \]

那么 \((a_0,a_1,\cdots,a_t)\) 可以与模 \(P\) 意义下的 \(a\) 一一对应，\(b\) 同理。于是可以直接乘起来。

于是现在就把问题转化成了「\(P\) 为一个奇素数」的问题。

为什么让你打表呢？你会惊奇地发现，对于同一个 \(P\)，\(x\neq 0\) 的答案都是一样的！哇好神奇！但是既然写题解，肯定不能告诉你直接打表就能发现……

不妨暴力地枚举 \(i\equiv a^2\)，则 \(b^2\equiv x-i\)，然后我们求出 \(i\) 和 \(x-i\) 的二次剩余数量，就可以求得 \(a^2\equiv i\) 时的解 \((a,b)\) 数量了。

怎么统计二次剩余的数量呢？让我们又来看看这个老朋友——勒让德符号，根据其实际意义不难得到 \(\left(\frac{a}{p}\right)+1\) 就是 \(x^2\equiv a\pmod p\) 的解的数量。

于是写出答案的式子：

\[\begin{aligned} &\sum_{i=0}^{p-1}\Big[\big(\tfrac{i}{p}\big)+1\Big]\Big[\big(\tfrac{x-i}{p}\big)+1\Big]\\ =&\sum_{i=0}^{p-1}\big(\tfrac{i}{p}\big)\big(\tfrac{x-i}{p}\big)+2\sum_{i=0}^{p-1}\big(\tfrac{i}{p}\big)+p \end{aligned} \]

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

\[\sum_{i=1}^{p-1}\big(\tfrac{i}{p}\big)\big(\tfrac{x-i}{p}\big)+2\sum_{i=1}^{p-1}\big(\tfrac{i}{p}\big)+p \]

由勒让德符号的计算式 \(\big(\frac{x}{p}\big)=x^{\frac{p-1}{2}}\)，不难得到 \(\big(\frac{x}{p}\big)\big(\frac{y}{p}\big)=\big(\frac{xy}{p}\big)\)。又有二次剩余的数量和非二次剩余的数量相等，所以有 \(\sum\limits_{i=1}^{p-1}\big(\frac{i}{p}\big)=0\)。上式可以继续简化：

\[\sum_{i=1}^{p-1}\big(\tfrac{i(x-i)}{p}\big)+p \]

根据下面的推导：

\[\big(\tfrac{xi^{-1}-1}{p}\big)\big(\tfrac{i^2}{p}\big)=\big(\tfrac{i(x-i)}{p}\big)=\big(\tfrac{xi^{-1}-1}{p}\big) \]

则：

\[\sum_{i=1}^{p-1}\big(\tfrac{xi^{-1}-1}{p}\big)+p \]

最后我们就可以知道为什么 \(x=0\) 的答案就和其他 \(x\neq0\) 的答案不一样——

• 当 \(x=0\) 时，答案为 \((p-1)\big(\tfrac{-1}{p}\big)+p=p+(p-1)(-1)^{\frac{p-1}{2}}\)；

• 当 \(x\neq0\) 时，由 \(i\) 从 \(1\) 枚举至 \(p-1\)，则 \(i^{-1}\) 取值也取尽 \(1\sim p-1\)，\(xi^{-1}\) 同理。

  那么 \(xi^{-1}-1\) 取尽 \(0\sim p-2\)。答案为

  \[\sum_{i=0}^{p-2}\big(\tfrac{i}{p}\big)+p=\sum_{i=0}^{p-1}\big(\tfrac{i}{p}\big)-\big(\tfrac{p-1}{p}\big)+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\n",ans);
	}
	return 0;
}

---

THE END

Thanks for reading!

你是春山 你是岁酒
你是自由 你是误谬
你是颠沛流离之后 我绮丽的愁
你是沙鸥 你是滴漏
你是白昼 你是不朽

——《你是我遥不可及的梦》By 苍穹/papaw泡泡

> Link 你是我遥不可及的梦-Bilibili