二次剩余 Cipoila 模板 P5491

前言

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

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

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

一、例题

以洛谷P5491为例 链接: P5491

二、思路及代码 1.思路

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

2.代码

代码如下（示例）：

#include 
#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);
  }
}