【Coel.做题笔记】【旁观者…】二次剩余- Cipolla 算法
题前闲语
这周末就是省选了，甚至考场就在这个机房，可惜我并没有参加的机会。
唉，今年得好好努力了！
题目简介
洛谷传送门
给出 \(N,p\)，求解方程
多组数据。
保证 \(p\) 是奇素数。
输入输出格式
输入格式
接下来 \(T\) 行，每行两个整数，分别是 \(N\) 和 \(p\) 。
输出格式
对于每一行输出：
若有解，则按 \(\bmod ~p\) 后递增的顺序输出在 \(\bmod~ p\) 意义下的全部解。
若两解相同，只输出其中一个。
若无解，则输出 Hola! 。
前置知识
之前学习了很多模意义下的操作，比如说离散对数（\(BSGS\)），同余方程（扩展欧几里得），乘法逆元（费马小定理），光速快速幂（扩展欧拉定理），同余方程组（中国剩余定理）等等。接下来看看模意义下的开方运算：二次剩余。
什么是二次剩余？
当关于 \(x\) 的同余方程
存在解时，则称 \(x\) 为模 \(p\) 意义下的二次剩余；反之则称为非二次剩余。
以下是几个二次剩余里用到的定理与符号：
勒让德符号
这是一个关于二次剩余的符号（以下的方程特指上文提到的同余方程）。
欧拉判别准则
欧拉：没想到吧又是我！
证明略过不提其实是懒得写。
这样我们就可以用快速幂计算出勒让德符号，进而判定二次剩余是否存在了。
二次剩余的个数
在模意义之下的二次剩余数与非二次剩余数相等，即为 \(\frac{p-1}{2}\)。
比较专业的证明需要用到拉格朗日定理，这里采用一种更为简单的证明方式。
假设原方程在模意义下的两个不同解分别为 \(x_1\) 和 \(x_2\) ，由平方差公式可以知道 \(p|(x_1+x_2)(x_1-x_2)\)。
当 \(p\) 为奇素数时，会有\((x_1-x_2) \bmod p \ne 0\)，那么\((x_1+x_2) \bmod p = 0\)，也就是说 \(x_1,x_2\)互为模意义下相反数，故共有 \(\frac{p-1}{2}\) 个二次剩余。
与此同时我们也可以知道，对于每个二次剩余而言，它对应的解的个数恰好为两个，且它们的和 \(x_1 +x_2\equiv 0\)。
Cipolla 算法
讲了这么多，终于到算法核心内容了。
首先我们用 rand() 求出一个 \(a\) 使得 \(a^2-n\) 为非二次剩余。由上面二次剩余个数可以知道，这个过程的期望复杂度为 \(O(2)\)。
接下来在模意义下暴力开根。定义 \(\omega ^2 \equiv a^2-n\), 类似虚数的操作，我们可以用 \(a+\omega b\) 表示这个模意义下的所有数。
此时原方程的解便是\((a+\omega)^{\frac{p-1}{2}}\) 我们来证明一下。
在模意义对解平方，则有可以表示为\(x^2 \equiv (a+\omega)^{p-1}\)。
暴力推柿子时间！
当然在推柿子之前还得来个引理：为什么不在前面全部说完
Lemma： \((a+b)^p \equiv a^p +b^p\)。用二项式定理展开以后，系数不为 \(1\) 的所有部分都会被模 \(p\) 消去，得证。
最后利用这个引理推柿子：
\(x^2 \equiv (a+\omega)^{p-1}\)
\(\equiv (a+\omega)^p(a+\omega)\)
\(\equiv (a-\omega)(a+\omega)\) （运用Lemma）
\(\equiv a^2 - \omega ^2\)
\(\equiv a^2 - (a^2 - n)\) （运用定义）
\(\equiv n\).
此外还可以保证 \(\omega\) 的系数必然为零，证明略。
\(Cipolla\) 算法的瓶颈在于快速幂，所以复杂度为 \(O(log n)\)。
代码实现
如果能理解透 \(Cipolla\) 算法的基本内容，实现就没什么难度了。
首先模仿虚数来定义一个数域，它的乘法定义和虚数是一样的：
struct Complex {
int x, y;
inline Complex(int x = 0, int y = 0): x(x), y(y) {}
};
inline Complex operator*(Complex a, Complex b) {
//乘法，注意取摸和溢出规避
int x = ((a.x * b.x % p + a.y * b.y % p * mul % p) % p + p) % p;
int y = ((a.x * b.y % p + a.y * b.x % p) % p + p) % p;
return Complex(x, y);
}
然后是整数和新数域下的快速幂：
int qpow_int(int a, int b, int p) {
int ans = 1;
while (b) {
if (b & 1) ans = ans * a % p;
a = a * a % p;
b >>= 1;
}
return ans;
}
int qpow_comp(Complex a, int b, int p) {
Complex ans = Complex(1, 0);
while (b) {
if (b & 1) ans = ans * a;
a = a * a;
b >>= 1;
}
return ans.x % p;
}
最后是算法主体的实现，还是有几个注意点的。
int Cipolla(int n, int p) {
n %= p;
if (qpow_int(n, (p - 1) / 2, p) == p - 1)//欧拉判别准则
return -1;
int a = rand() % p;
mul = ((a * a % p - n) % p + p) % p; //mul为a^2-n
while (qpow_int(mul, (p - 1) / 2, p) != p - 1){ //随机数找合法的a
a = rand() % p;
mul = ((a * a % p - n) % p + p) % p;
}
return qpow_comp(Complex(a, 1), (p + 1) / 2, p);
}
完整代码：
#include <cstdio>
#include <cctype>
#include <iostream>
#include <cstdlib>
#include <ctime>
#define int long long
#define WYSI 727727727 //奇怪的随机数种子
namespace FastIO {
#ifdef ONLINE_JUDGE
#define _getchar_nolock getchar_unlocked
#define _putchar_nolock putchar_unlocked
#endif
inline int read() {
int x = 0, f = 1;
char ch = _getchar_nolock();
while (!isdigit(ch)) {
if (ch == '-')
f = -1;
ch = _getchar_nolock();
}
while (isdigit(ch)) {
x = x * 10 + ch - '0';
ch = _getchar_nolock();
}
return x * f;
}
inline void write(int x) {
static int buf[35];
int top = 0;
if (x < 0) {
x = -x;
_putchar_nolock('-');
}
do {
buf[top++] = x % 10;
x /= 10;
} while (x);
while (top)
_putchar_nolock(buf[--top] + '0');
_putchar_nolock(' ');
}
} // namespace FastIO
using namespace FastIO;
using namespace std;
int T, n, p;
int mul;
struct Complex {
int x, y;
inline Complex(int x = 0, int y = 0): x(x), y(y) {}
};
inline Complex operator*(Complex a, Complex b) {
int x = ((a.x * b.x % p + a.y * b.y % p * mul % p) % p + p) % p;
int y = ((a.x * b.y % p + a.y * b.x % p) % p + p) % p;
return Complex(x, y);
}
int qpow_int(int a, int b, int p) {
int ans = 1;
while (b) {
if (b & 1) ans = ans * a % p;
a = a * a % p;
b >>= 1;
}
return ans;
}
int qpow_comp(Complex a, int b, int p) {
Complex ans = Complex(1, 0);
while (b) {
if (b & 1) ans = ans * a;
a = a * a;
b >>= 1;
}
return ans.x % p;
}
int Cipolla(int n, int p) {
n %= p;
if (qpow_int(n, (p - 1) / 2, p) == p - 1)
return -1;
int a = rand() % p;
mul = ((a * a % p - n) % p + p) % p;
while (qpow_int(mul, (p - 1) / 2, p) != p - 1){
a = rand() % p;
mul = ((a * a % p - n) % p + p) % p;
}
return qpow_comp(Complex(a, 1), (p + 1) / 2, p);
}
signed main(void) {
srand(WYSI);
T = read();
while (T--) {
n = read(), p = read();
if (n == 0) { //n=0时答案自然为0
write(n);
puts("");
continue;
}
int ans1 = Cipolla(n, p), ans2 = p - ans1;
if (ans1 == -1) {
puts("Hola!");
} else {
if (ans1 == ans2) {
write(ans1);
puts("");
} else {
write(min(ans1, ans2));
write(max(ans1, ans2));
puts("");
}
}
}
return 0;
}
题后闲话
最近还是感觉有点摸鱼了，得好好加油才行！
浙公网安备 33010602011771号