二次剩余 Cipoila 模板 P5491

---

前言

本文介绍求解二次剩余的Cipolla模板
求解二次剩余问题一般有两个方向

比较常规的方向是用 Cipolla 算法 或 Legendre 算法 或 Tonelli-Shanks 算法进行计算
这三个算法中Cipolla使用最为广泛, 严谨证明过程可参考oiwiki-Cipolla
Legendre算法记载不详、Tonelli算法用到了原根相关知识
而Cipolla则关键思想是用到了模p域的虚根

另一个方向是利用BSGS大小步算法配合原根解决，可参考oiwiki BSGS进阶篇

一、例题

以洛谷P5491为例 链接: P5491

二、思路及代码

1.思路

我们的问题是求$$x^2\equiv N(\mathrm{m}\mathrm{o}\mathrm{d}p)$$的解，Cipolla方法是基于虚根的构造。我们首先找一个 a 使得$${\omega }^2=a^2-N$$为非二次剩余 ，这样的 a 很好找，[1, p-1]之间有一半是满足的。那么则有：$${\omega }^{p-1}\equiv (a^2-N{)}^{\frac{p-1}{2}}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$此时的 w 便是我们所构造的虚根，那我们断言：$$x\equiv (a+\omega {)}^{\frac{p+1}{2}}(\mathrm{m}\mathrm{o}\mathrm{d}p)$$即为我们的解，原因如下：$$\begin{gathered} x^2\equiv (a+w{)}^{p+1} \\ \equiv (a+\omega )(a^p+{\omega }^p) \\ \equiv (a+\omega )(a-\omega ) \\ \equiv (a^2-{\omega }^2)\equiv N(\mathrm{m}\mathrm{o}\mathrm{d}p) \end{gathered}$$
现在的问题是我们如何求出这个具体的 x ，方法是利用虚根，即：$$(a+b\omega )\ast (c+d\omega )=ac+bd{\omega }^2+(ad+bc)\omega$$我们最后只需取其实部即可，再将$${\omega }^2=a^2-N$$代入复数运算即可。

2.代码

代码如下（示例）：

#include <iostream>
#define int long long
using namespace std;
int p, ww;
struct complex {
  int x, y;
  complex(int _x = 0, int _y = 0) : x(_x), y(_y) {}
};
complex operator+(complex &a, complex &b) {
  return complex((a.x + b.x) % p, (a.y + b.y) % p);
}
complex operator-(complex &a, complex &b) {
  return complex((a.x - b.x + p) % p, (a.y - b.y + p) % p);
}
complex operator*(complex &a, complex &b) {
  return complex((a.x * b.x % p + a.y * b.y % p * ww % p + p) % p,
                 (a.x * b.y + a.y * b.x) % p);
}  // (a + bw) * (c + dw) = ac + bdww + w(ad + bc)
complex quickPow(complex x, int n) {
  complex ans = complex(1, 0);
  while (n) {
    if (n & 1) ans = ans * x;
    x = x * x;
    n >>= 1;
  }
  return ans;
}
int quickPow(int x, int n) {  // 计算 Legendre
  int ans = 1;
  while (n) {
    if (n & 1) ans = ans * x % p;
    x = x * x % p;
    n >>= 1;
  }
  return ans;
}
int Legendre(int x) {  // 1: 二次剩余, p-1: 非二次剩余
  return quickPow((x % p + p) % p, (p - 1) / 2);
}
complex Cipolla(int n) {
  int a = rand() % p;  // w 为单位根 ww = w^2,  w ^ (p - 1) = -1 (mod p)
  while (Legendre(a * a - n) == 1) a = rand() % p;
  ww = ((a * a - n) % p + p) % p;  // ww = a^2 - n (mod p)
  complex ans = quickPow(complex(a, 1), (p + 1) / 2);
  return ans;  // ans = (a + w) ^ (p + 1) / 2
}
signed main() {
  //   freopen("in.txt", "r", stdin);
  //   freopen("out.txt", "w", stdout);
  int t, n;
  scanf("%lld", &t);
  while (t--) {
    scanf("%lld%lld", &n, &p);
    if (!n) {
      printf("0\n");
      continue;
    }
    if (p == 2) {
      printf("%lld\n", n);
      continue;
    }
    if (Legendre(n) == p - 1) {
      printf("Hola!\n");
      continue;
    }
    complex ans = Cipolla(n);
    int x1 = ans.x;
    int x2 = p - x1;
    if (x1 > x2) swap(x1, x2);
    if (x1 != x2)
      printf("%lld %lld\n", x1, x2);
    else
      printf("%lld\n", x1);
  }
}