Timus Online Judge传送门
Vjudge传送门
Description
The number x is called a square root of a modulo n (root( a, n)) if x * x = a (mod n). Write the program to find the square root of number a by given modulo n.
Input
One number K in the first line is an amount of tests ( K ≤ 100000). Each next line represents separate test, which contains integers a and n (1 ≤ a, n ≤ 32767, n is prime, a and n are relatively prime).
Output
For each input test the program must evaluate all possible values root( a, n) in the range from 1 to n − 1 and output them in increasing order in one separate line using spaces. If there is no square root for current test, the program must print in separate line: ‘No root’.
Example
Input
5
4 17
3 7
2 7
14 31
10007 20011
Output
2 15
No root
3 4
13 18
5382 14629
解题分析
二次剩余模板题, 这里用 C i p o l l a Cipolla Cipolla算法做, 只讨论奇质数的情况。
首先我们有欧拉准则:
∃
x
,
x
2
≡
a
(
m
o
d
p
)
⇐
⇒
a
p
−
1
2
≡
1
(
m
o
d
p
)
\exists x, x^2\equiv a(mod\ p)\Leftarrow\Rightarrow a^{\frac{p-1}{2}}\equiv 1(mod\ p)
∃x,x2≡a(mod p)⇐⇒a2p−1≡1(mod p)
证明:
-
充分性:
a p − 1 2 ≡ ( x 2 ) p − 1 2 ≡ x p − 1 ≡ 1 ( m o d p ) a^{\frac{p-1}{2}} \\ \equiv (x^2)^{\frac{p-1}{2}} \\ \equiv x^{p-1} \\ \equiv1(mod\ p) a2p−1≡(x2)2p−1≡xp−1≡1(mod p) -
必要性:
设 g g g为 p p p的原根, g i ≡ a ( m o d p ) g^{i}\equiv a(mod\ p) gi≡a(mod p)。
那么有 g i ( p − 1 ) 2 ≡ 1 ( m o d p ) g^{\frac{i(p-1)}{2}}\equiv 1(mod\ p) g2i(p−1)≡1(mod p)
就有 ( p − 1 ) ∣ i ( p − 1 ) 2 (p-1)|\frac{i(p-1)}{2} (p−1)∣2i(p−1), 即 i i i为偶数。
那么 x x x有一组解为 g i 2 g^{\frac{i}{2}} g2i。
对于一组 a , n a,n a,n, 我们已经可以用欧拉准则判定是否有解(特判2)。
现在我们先需要找到一组满足 b 2 − a = w b^2-a=w b2−a=w, 且w为非二次剩余的 b b b和 w w w。这个也很好做, 只需要瞎 r a n d ( ) rand() rand()b, 然后判 w p − 1 2 w^{\frac{p-1}{2}} w2p−1是否在 m o d p mod\ p mod p意义下等于 − 1 -1 −1就好了。
不妨扩域, 设 i = w i=\sqrt{w} i=w, 将一个数表示为 a + b i a+bi a+bi, 类似于复数。
定义代数系统
<
G
,
+
,
×
>
<G,+,\times>
<G,+,×>满足:
(
a
1
,
b
1
i
)
+
(
a
2
,
b
2
i
)
=
(
a
1
+
a
2
,
(
b
1
+
b
2
)
i
)
(
a
1
,
b
1
i
)
∗
(
a
2
,
b
2
i
)
=
(
a
1
a
2
+
b
1
b
2
w
,
(
a
1
b
2
+
a
2
b
1
)
i
)
(a_1,b_1i)+(a_2,b_2i)=(a_1+a_2,(b_1+b_2)i) \\ (a_1,b_1i)*(a_2,b_2i)=(a_1a_2+b_1b_2w,(a_1b_2+a_2b_1)i)
(a1,b1i)+(a2,b2i)=(a1+a2,(b1+b2)i)(a1,b1i)∗(a2,b2i)=(a1a2+b1b2w,(a1b2+a2b1)i)
那么可以发现,
G
G
G实际上是一个环, 满足结合律、交换律。
然后又有一个神仙结论: x = ( b + i ) p + 1 2 x=(b+i)^{\frac{p+1}{2}} x=(b+i)2p+1。
证明(以下都在模
p
p
p意义下讨论):
x
2
=
(
b
+
i
)
p
+
1
=
(
b
+
i
)
p
(
b
+
i
)
=
∑
k
=
0
p
(
p
k
)
b
k
i
p
−
k
(
b
+
i
)
x^2= (b+i)^{p+1} \\ =(b+i)^p(b+i) \\ =\sum_{k=0}^p\binom{p}{k}b^ki^{p-k}(b+i)
x2=(b+i)p+1=(b+i)p(b+i)=k=0∑p(kp)bkip−k(b+i)
发现
(
p
k
)
\binom{p}{k}
(kp)当且仅当
k
=
0
k=0
k=0或
k
=
p
k=p
k=p的时候后不被
p
p
p整除, 那么有
∑
k
=
0
p
(
p
k
)
b
k
i
p
−
k
(
b
+
i
)
=
(
b
p
+
i
p
)
(
b
+
i
)
=
(
b
+
i
p
)
(
b
+
i
)
=
(
b
+
i
×
(
i
2
)
p
−
1
2
)
(
b
+
i
)
=
(
b
+
i
×
w
p
−
1
2
)
(
b
+
i
)
=
(
b
−
i
)
(
b
+
i
)
=
b
2
−
w
\sum_{k=0}^p\binom{p}{k}b^ki^{p-k}(b+i) \\ =(b^p+i^p)(b+i) \\ =(b+i^p)(b+i) \\ =(b+i\times (i^2)^{\frac{p-1}{2}})(b+i) \\ =(b+i\times w^{\frac{p-1}{2}})(b+i) \\ =(b-i)(b+i) \\ =b^2-w
k=0∑p(kp)bkip−k(b+i)=(bp+ip)(b+i)=(b+ip)(b+i)=(b+i×(i2)2p−1)(b+i)=(b+i×w2p−1)(b+i)=(b−i)(b+i)=b2−w
得证。
代码如下:
#include <cstdio>
#include <cmath>
#include <cctype>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define R register
#define IN inline
#define W while
#define gc getchar()
#define ll long long
struct QF {ll re, im, w;};
struct INFO {ll b, w;};
IN QF operator * (const QF &x, const QF &y) {return {x.re * y.re + x.w * y.im * x.im, x.im * y.re + y.im * x.re, x.w};}
IN QF operator + (const QF &x, const QF &y) {return {x.re + y.re, x.im + y.im, x.w};}
IN QF fpow(QF x, R ll tim, R ll mod)
{
QF ret = {1, 0, x.w};
W (tim)
{
if (tim & 1) ret = ret * x;
x = x * x, tim >>= 1;
ret.re %= mod, ret.im %= mod;
x.re %= mod, x.im %= mod;
}
return ret;
}
IN ll fpow(R ll base, R ll tim, R ll mod)
{
ll ret = 1;
W (tim)
{
if (tim & 1) ret = ret * base % mod;
base = base * base % mod, tim >>= 1;
}
return ret;
}
IN bool check(R int val, R int n) {return fpow(val, (n - 1) / 2, n) == 1;}
IN INFO find(R int a, R int n)
{
ll b;
W (233)
{
b = rand() % n;
if (!check((b * b % n - a + n) % n, n)) return {b, (b * b % n - a + n) % n};
}
}
int main(void)
{
using namespace std;
int T, n, a; INFO dat; ll res;
ios_base::sync_with_stdio(false);
cin >> T;
W (T--)
{
cin >> a >> n; a %= n;
if (n == 2) {cout << "1" << endl; continue;}
if (!check(a, n)) {cout << "No root" << endl; continue;}
dat = find(a, n);
res = fpow({dat.b, 1, dat.w} ,((n + 1) >> 1), n).re;
if (res * 2 == n) cout << res << endl;
else
{
if (res * 2 > n) res = n - res;
cout << res << " " << n - res << endl;
}
}
}