二次剩余Cipolla法

翻着翻着居然翻到了两个退役选手zy和mw的博客，想到自己可能也会像他们一样省选失利……

---

由于我比较水，目前只会模数是质数的情况，之后会了可能会补。

很多具体的证明我也不会，限于感性认识，所以请看zy的博客：
https://blog.youkuaiyun.com/a_crazy_czy/article/details/51959546

---

二次剩余要做什么：

给出$n,p$，p为一质数。

要求$x$,满足$x^2=n(modp)$

判断无解，或者输出所有的解。

---

以下的口胡均不考虑n=0或p=2的情况。

首先明确一个点：

如果有$x^2=n(modp)$，那么也有$(p-x{)}^2=n(modp)$，并且不会有其它的解。

这个可以证明，假设有其它的解$y$

那么$x^2-y^2=0(modp)$
$(x+y)(x-y)=0(modp)$，所以有$y=\pm x$

---

勒让德符号：
$\left(\frac{a}{p}\right)=\{\begin{array}{ll}1,& a为p的二次剩余\\ -1,& a为p的二次非剩余\\ 0,& a\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}p)\end{array}$

且$\left(\frac{a}{p}\right)=a^{(p-1)/2}(modp)$
1.
费马小定理有任意$n^{p-1}=1(modp)$
所以若有$x^2=n(modp)$，那么$x^{p-1}=n^{(p-1)/2}=1$
这明显是个必要条件，但是充分我不会证。
2.
反过来，如果$n^{(p-1)/2}=-1$,那一定没有解

---

做法：
随机一个$a$，设$w=a^2-n(modp)$，若w没有二次剩余。

那么$x1=(a+\sqrt{w}{)}^{(p+1)/2}$。

结论1：
$(a+\sqrt{w}{)}^p$
$={\sum }_{i=0}^p{C}_p^i\ast a^i\ast {\sqrt{w}}^{p-i}$
$=a^p+{\sqrt{w}}^p$
$=a-\sqrt{w}$

$(a+\sqrt{w}{)}^{p+1}$
$=(a+\sqrt{w})\ast (a-\sqrt{w})$
$=a^2-w=n$

那么$n^{1/2}=(a+\sqrt{w}{)}^{(p+1)/2}$

如何计算：
$\sqrt{w}$显然是不存在的，那么可以当作$i$来做复数乘法。

据说做完之后$\sqrt{w}$的系数一定是0，我又不会证明了，hhh…

随机a的次数是O(2)的。

例题：
https://www.codechef.com/problems/FN

这个还需要bsgs。

#include<map>
#include<ctime>
#include<cmath>
#include<cstdio>
#include<algorithm>
#define pp printf
#define ll long long
#define fo(i, x, y) for(int i = x; i <= y; i ++)
using namespace std;

int rand(int x, int y) {
	return ((ll) RAND_MAX * rand() + rand()) % (y - x + 1) + x;
}

ll mo;
int Q;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

int pd(ll n) {
	if(!n) return 0;
	return ksm(n, (mo - 1) / 2) == 1;
}

ll w;

struct P {
	ll x, y;
	P() {}
	P(ll _x, ll _y) { x = _x, y = _y;}
};

P operator +(P a, P b) {
	return P((a.x + b.x) % mo, (a.y + b.y) % mo);
}
P operator *(P a, P b) {
	return P((a.x * b.x + a.y * b.y % mo * w) % mo, (a.x * b.y + a.y * b.x) % mo);
}

P ksm(P x, ll y) {
	P s = P(1, 0);
	for(; y; y /= 2, x = x * x)
		if(y & 1) s = s * x;
	return s;
}

ll Sqrt(ll n) {
	if(!n) return 0;
	if(mo == 2) return 1;
	if(!pd(n)) return -1;
	while(1) {
		ll a = rand(0, mo - 1);
		w = (a * a % mo - n + mo) % mo;
		if(pd(w)) continue;
		return ksm(P(a, 1), (mo + 1) / 2).x;
	}
}

map<int, int> q[2];

ll bsgs(ll a, ll b, int k) {
	q[0].clear(); q[1].clear();
	int m = ceil(sqrt(mo));
	ll x = b, am = ksm(a, m);
	fo(i, 1, m) {
		x = x * a % mo;
		q[i & 1][x] = i;
	}
	x = 1;
	fo(i, 1, m) {
		x = x * am % mo;
		if(q[k ^ ((i * m) & 1)][x]) return i * m - q[k ^ ((i * m) & 1)][x];
	}
	return -1;
}

ll a, t, n, ni2, ans;

int main() {
	srand(time (0));
	freopen("a.in", "r", stdin);
	for(scanf("%d", &Q); Q; Q --) {
		scanf("%lld %lld", &n, &mo);
		ni2 = ksm(2, mo - 2);
		a = Sqrt(5); t = (1 + a) * ni2 % mo;
		ans = mo;
		n = a * n % mo;
		ll dat = n * n - 4; dat = (dat % mo + mo) % mo;
		ll q = Sqrt(dat);
		if(q != -1) {
			ll T = (n + q) * ni2 % mo;
			ll x = bsgs(t, T, 1);
			if(x != -1) ans = x;
			T = (n - q + mo) * ni2 % mo;
			x = bsgs(t, T, 1);
			if(x != -1) ans = min(ans, x);
		}
		dat = n * n + 4; dat = (dat % mo + mo % mo);
		q = Sqrt(dat);
		if(q != -1) {
			ll T = (n + q) * ni2 % mo;
			ll x = bsgs(t, T, 0);
			if(x != -1) ans = min(ans, x);
			T = (n - q + mo) * ni2 % mo;
			x = bsgs(t, T, 0);
			if(x != -1) ans = min(ans, x);
		}
		pp("%lld\n", ans == mo ? -1 : ans);
	}
}