什么是二次剩余
对于一个奇素数(p),和一个整数(nin [0,p)),如果同余方程(x^2equiv npmod p)有解,那么称(n)是(p)的一个二次剩余。
关于二次剩余,专门有一个关于它的“勒让德记号”:
并且关于勒让德记号还有一个欧拉判定法:
勒让德记号与欧拉判定法
给出证明如下:
-
(a^{frac{p-1}2}equiv 1pmod p)
-
充分性:如果 (a) 是 (p) 的二次剩余,则 (a^{frac{p-1}2}equiv 1pmod p) 。
显然,根据定义,存在整数 (xin[0,p)) ,使得 (x^2equiv apmod p) 。代入得到:
[x^{p-1}equiv 1pmod p ]根据费马小定理显然成立。
-
必要性:如果 (a^{frac{p-1}2}equiv 1pmod p) ,则 (a) 是 (p) 的一个二次剩余。
由于 (p) 是奇素数,所以可以找到它的原根 (g) ,且必然存在一个 (k) ,使得 (g^kequiv apmod p) 。
代入得到:
[g^{frac{k(p-1)}2}equiv 1pmod p ]由于 (g^{p-1}equiv 1pmod p) ,所以 (p-1|frac{k(p-1)}2) 。
由于 (frac{k(p- 1)}2) 和 (p-1) 都是整数,所以 (frac k 2) 也是整数。所以 (2|k) 。
我们便可以找到一个 (xequiv g^{frac k2}pmod p) ,这样就找到了 (x^2equiv apmod p) 。
-
-
(a^{frac{p-1}2}equiv 0pmod p)
当然,如果 (p|a) ,这个判断法是显然成立的。
-
(a^{frac{p-1}2}equiv -1pmod p)
接着,证明 (a) 不是 (p) 的二次剩余的情况。
-
充分性:如果 (a) 不是 (p) 的二次剩余,则(a^{frac{p-1}2}equiv -1pmod p)。
由于 (p) 是一个奇素数,所以 (forall iin[0,p)cap mathbb{Z},exists jin[0,p)capmathbb Z,s.t.ijequiv apmod p) 。
于是我们可以将 (p-1) 个数分成 (frac{p-1}2) 组,每组的积为 (a) 。
将这 (frac{p-1}2) 个组全部乘起来,得到了 ((p-1)!equiv a^{frac{p-1}2}pmod p) ,而根据威尔逊定理可以得到 ((p-1)!equiv -1pmod p) ,因此 (a^{frac{p-1}2}equiv -1pmod p) 。
-
于是这个等式就被证明了。
有用的性质
还有一个有用的性质:
奇素数 (p) 的二次剩余恰好有 (frac{p-1}2) 个。该性质等价于, (p-1) 个数中可以按照平方分组,每组有且仅有两个元素。
首先,证明每一对和为 (p) 的数,它们的平方相等。
对于 (u,vin[1,p)capmathbb Z,u>v) ,假设 (u^2equiv v^2pmod p) ,所以 (p|(u+v)(u-v)) 。
由于 (u-vin[1,p-2]) ,所以 (p|u+v) 。
然后考虑不同组 ((u,v)) 和 ((w,x)) 。
根据上面有 (u+v=p,w+x=p) 。假设 (u^2equiv w^2pmod p) ,所以还有 (u+w=p) 。
所以 (v=w) ,所以 (u=x) ,所以 ((u,v)) 就是 ((w,x)) ,与假设不符。所以每一组数有且仅有两个数。
Cipolla 算法
这个算法可以用于求解 (x^2equiv npmod p) 的方程的解(前提当然是(n)是(p)的二次剩余)
流程如下:
-
随机一个 (ain[0,p)cap mathbb{Z}) ,计算 (b=a^2-n) 。如果 (left(frac b p ight)=-1) ,进行 2 ,否则重复 1 。
-
这是本算法最富想象力的一个步骤——我们要对 (b) 开方!具体过程我们待会儿慢慢港,现在先设这个根为 (omega=sqrt b) 。
-
得到了 (omega) 之后,答案就是 ((a+omega)^{frac{p+1}2}) ,以及 (p-(a+omega)^{frac{p+1}2}) (前文提到过)。
首先说明 (3) 这一步骤。
根据二项式定理我们可以展开 ((a+omega)^p) :
考虑放在模意义下,由于 (inom{p}{i}=frac{p!}{i!(p-i)!}) ,且 (iin[0,p]capmathbb Z) ,所以如果 (iin(0,p)capmathbb Z) ,那么 (inom p i) 中的 (p) 就不会被约掉。换言之, (forall iin(0,p)mathbb Z,p|inom p i) 。
于是有:
根据费马小定理可以继续简化得到:
由于 (left(frac b p ight)=-1) ,所以 (b^{frac{p-1}2}equiv -1pmod p) ,即 (omega^{p-1}equiv -1pmod p) ,即 (omega^pequiv -omegapmod p) 。于是可以继续简化:
由于 (2|p+1) ,所以必然存在 ((a+omega)^{frac{p+1}2}) 。所以这就是其中的一个解。
现在来解释一下 (omega) 是怎么算出来的。
既然我们不能在模意义下直接开方,我们就直接暴力计算。
简单来说——就像虚数单位 (i^2=-1) 一样,我们定义一个单位 (omega^2=a^2-b) ,然后就可以构造一个“复数” (a+bomega) 。
相当于是将 (sqrt b) 扩入了模 (p) 整数域中。
我们可以得到这一种数的运算规律:
可以参考来自博客二次剩余Cipolla算法学习小记(博主SFN1036 On CSDN)(如有侵权可以联系删除)的图片:
知道了这个之后就可以写一个结构体或者类正常运算了。
例题
来自URAL 1132的 Square Root,入门模板题。代码:
#include <cstdio>
#include <cstdlib>
typedef long long LL;
#define int LL
#define random myRandom
template<typename _T>
void read( _T &x )
{
x = 0;char s = getchar();int f = 1;
while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}
while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}
x *= f;
}
template<typename _T>
void write( _T x )
{
if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }
if( 9 < x ){ write( x / 10 ); }
putchar( x % 10 + '0' );
}
template<typename _T>
_T MIN( const _T a, const _T b )
{
return a < b ? a : b;
}
template<typename _T>
_T MAX( const _T a, const _T b )
{
return a > b ? a : b;
}
int N, mod, w;
struct comp
{
int x, y;
comp() { x = y = 0; }
comp( const int a ) { x = a, y = 0; }
comp( const int a, const int b ) { x = a, y = b; }
comp operator * ( const comp & b ) const
{ return comp( ( 1ll * x * b.x % mod + 1ll * y * b.y % mod * w % mod ) % mod,
( 1ll * x * b.y % mod + 1ll * y * b.x % mod ) % mod ); }
void operator *= ( const comp & b ) { *this = *this * b; }
};
int qkpow( int base, int indx )
{
int ret = 1;
while( indx )
{
if( indx & 1 ) ret = 1ll * ret * base % mod;
base = 1ll * base * base % mod, indx >>= 1;
}
return ret;
}
comp qkpow( comp base, int indx )
{
comp ret = 1;
while( indx )
{
if( indx & 1 ) ret *= base;
base *= base, indx >>= 1;
}
return ret;
}
int random() { return rand() % mod; }
int inv( const int a ) { return qkpow( a, mod - 2 ); }
int fix( const int a ) { return ( a % mod + mod ) % mod; }
int chk( const int a )
{
int t = qkpow( a, mod - 1 >> 1 );
if( t + 1 == mod ) return -1;
return t;
}
signed main()
{
comp t;
int T, a;
read( T );
while( T -- )
{
read( N ), read( mod );
if( mod == 2 ) { puts( "1" ); continue; }
N %= mod;
if( chk( N ) <= 0 ) { puts( "No root" ); continue; }
srand( mod );
while( true )
{
a = random();
w = fix( 1ll * a * a % mod - N );
if( chk( w ) < 0 ) break;
}
t = comp( a, 1 );
t = qkpow( t, mod + 1 >> 1 );
if( t.x * 2 == mod ) { write( t.x ), putchar( '
' ); continue; }
int ans1 = t.x, ans2 = mod - t.x;
write( MIN( ans1, ans2 ) ), putchar( ' ' ), write( MAX( ans1, ans2 ) ), putchar( '
' );
}
return 0;
}