【模板】二次剩余：Cipolla 算法

给定 \(c, p\)，\(p\) 为奇质数，求解关于 \(x\) 的同余方程 \(x^2\equiv c\pmod p\)。

欧拉判别：对于任意 \(c\)，\(c^{(p-1)/2}\equiv \pm 1\pmod p\)。当且仅当 \(c\) 有二次剩余时，\(c^{(p-1)/2}\equiv 1\pmod p\)。

Cipolla 过程：随机一个 \(a\)，使得 \(a^2-c\) 没有二次剩余（期望 \(O(1)\) 次找到）。定义 \(\mathbf i\)，满足 \(\mathbf i^2 = a^2 - c\)（这个 \(\mathbf i\) 是不存在的，使用扩域，将数字表示为 \(u+v\mathbf i\)）。

则 \(x_1=(a+\mathbf i)^{(p+1)/2}\)，另一个解是它的相反数。使用快速幂解决，时间复杂度 \(O(\log p)\)。

证明：

引理：\(\mathbf i^p=-\mathbf i\)（注意 \(\mathbf i\) 不在剩余系中，自然不遵守费马小定理）。证明为 \(\mathbf i^p=\mathbf i(\mathbf i^2)^{(p-1)/2}=-\mathbf i\)（因为 \(\mathbf i^2=a^2-c\) 没有二次剩余，判别式为 \(-1\)）。

所以答案的平方为

\[(a+\mathbf i)^{p+1}=(a+\mathbf i)^p(a+\mathbf i)=(a^p+\mathbf i^p)(a+\mathbf i)=(a-\mathbf i)(a+\mathbf i)=a^2-\mathbf i^2=c \]

点击查看代码

#include <bits/stdc++.h>

#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) void(0)
#endif
typedef long long LL;
template <int id>
struct dynamic_modint {
  typedef dynamic_modint modint;
  static unsigned umod;
  static int mod;
  static void set_mod(int _mod) { umod = mod = _mod; }
  unsigned v;
  dynamic_modint() : v(0) {}
  template <class T>
  dynamic_modint(T x) {
    x %= mod;
    v = x < 0 ? x + umod : x;
  }
  modint operator+() const { return *this; }
  modint operator-() const { return modint(0) - *this; }
  friend int raw(const modint &self) { return self.v; }
  template <class T>
  friend modint qpow(modint a, T b) {
    modint r = 1;
    for (; b; b >>= 1, a *= a)
      if (b & 1) r *= a;
    return r;
  }
  modint &operator+=(const modint &rhs) {
    if (v += rhs.v, v >= umod) v -= umod;
    return *this;
  }
  modint &operator-=(const modint &rhs) {
    if (v -= rhs.v, v >= umod) v += umod;
    return *this;
  }
  modint &operator*=(const modint &rhs) {
    v = 1ull * v * rhs.v % umod;
    return *this;
  }
  modint &operator/=(const modint &rhs) {
    assert(rhs);
    return *this *= qpow(rhs, mod - 2);
  }
  friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
  friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
  friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
  friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
  bool operator==(const modint &rhs) const { return v == rhs.v; }
  bool operator!=(const modint &rhs) const { return v != rhs.v; }
};
template <int id>
unsigned dynamic_modint<id>::umod;
template <int id>
int dynamic_modint<id>::mod;
typedef dynamic_modint<-1> mint;
mint solve(mint c) {
  static const auto check = [](mint c) {
    return qpow(c, (mint::mod - 1) >> 1) == 1;
  };
  if (c == 0) return 0;
  if (!check(c)) throw "No solution!";
  static mt19937 rng{random_device{}()};
  mint a = rng();
  while (check(a * a - c)) a = rng();
  typedef pair<mint, mint> number;
  const auto mul = [=](number x, number y) {
    return make_pair(x.first * y.first + x.second * y.second * (a * a - c),
                     x.first * y.second + x.second * y.first);
  };
  const auto qpow = [=](number a, int b) {
    number r = {1, 0};
    for (; b; b >>= 1, a = mul(a, a))
      if (b & 1) r = mul(r, a);
    return r;
  };
  return qpow({a, 1}, (mint::mod + 1) >> 1).first;
}
int main() {
#ifndef NF
  cin.tie(nullptr)->sync_with_stdio(0);
#endif
  int T;
  cin >> T;
  while (T--) {
    int n, p;
    cin >> n >> p;
    mint::set_mod(p);
    try {
      mint ret = solve(n);
      if (ret == 0) {
        cout << 0 << endl;
      } else {
        int x = raw(ret), y = raw(-ret);
        cout << min(x, y) << " " << max(x, y) << endl;
      }
    } catch (const char *) {
      cout << "Hola!" << endl;
    }
  }
  return 0;
}