二次剩余入门

什么是二次剩余

对于一个奇素数(p)，和一个整数(nin [0,p))，如果同余方程(x^2equiv npmod p)有解，那么称(n)是(p)的一个二次剩余。

关于二次剩余，专门有一个关于它的“勒让德记号”：

[left(frac a p ight)=egin{cases}0&a=0\1&exists xin[0,p)cap mathbb{Z},x^2equiv apmod p\-1& extrm{otherwise}end{cases} ]

并且关于勒让德记号还有一个欧拉判定法：

[left(frac a p ight)equiv a^{frac{p-1}2}pmod p ]

勒让德记号与欧拉判定法

给出证明如下：

(a^{frac{p-1}2}equiv 1pmod p)
1. 充分性：如果 (a) 是 (p) 的二次剩余，则 (a^{frac{p-1}2}equiv 1pmod p) 。
  
  显然，根据定义，存在整数 (xin[0,p)) ，使得 (x^2equiv apmod p) 。代入得到：
  
  [x^{p-1}equiv 1pmod p ]
  根据费马小定理显然成立。
2. 必要性：如果 (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) 的二次剩余的情况。
1. 充分性：如果 (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) ：

[(a+omega)^p=sum_{i=0}^{p}inom p ia^iomega^{p-i} ]

考虑放在模意义下，由于 (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) 。
于是有：

[(a+omega)^pequiv a^p+omega^ppmod p ]

根据费马小定理可以继续简化得到：

[(a+omega)^pequiv a+omega^p pmod p ]

由于 (left(frac b p ight)=-1) ，所以 (b^{frac{p-1}2}equiv -1pmod p) ，即 (omega^{p-1}equiv -1pmod p) ，即 (omega^pequiv -omegapmod p) 。于是可以继续简化：

[egin{aligned} (a+omega)^p&equiv a-omega&pmod p\ (a+omega)^{p+1}&equiv (a-omega)(a+omega)&pmod p\ &equiv a^2-b&pmod p\ &equiv n&pmod p end{aligned} ]

由于 (2|p+1) ，所以必然存在 ((a+omega)^{frac{p+1}2}) 。所以这就是其中的一个解。

现在来解释一下 (omega) 是怎么算出来的。

既然我们不能在模意义下直接开方，我们就直接暴力计算。

简单来说——就像虚数单位 (i^2=-1) 一样，我们定义一个单位 (omega^2=a^2-b) ，然后就可以构造一个“复数” (a+bomega) 。

相当于是将 (sqrt b) 扩入了模 (p) 整数域中。

我们可以得到这一种数的运算规律：

[(p+qomega)+(r+somega)=(p+q)+(r+s)omega ]

[(p+qomega)(r+somega)=(pr+qsb)+(ps+qr)omega ]

可以参考来自博客二次剩余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;
}