转自：http://blog.csdn.net/doyouseeman/article/details/52033204
简介
Cipolla算法是解决二次剩余强有力的工具，一个脑洞大开的算法。
认真看懂了，其实是一个很简单的算法，不过会感觉得出这个算法的数学家十分的机智。
基础数论储备
二次剩余
首先来看一个式子),我们现在给出n，要求求得x的值。如果可以求得，n为mod p的二次剩余，其实就是n在mod p意义下开的尽方。Cipolla就是一个用来求得上式的x的一个算法。
勒让德符号
勒让德符号是判断一个数是否为p的二次剩余的一个有力工具，p一定要为奇质数。（n,p）表示n为关于p的勒让德符号。其实就是判断n是否为p的二次剩余。
看起来好像是一大段废话，勒让德只是站在巨人的肩膀上总结了一下而已。
其实勒让德还总结了一些性质，不过一般只用得到欧拉判别准则，不够勒让德符号是个积性函数可能还有点用。
但还是不知道如何判断n是否为p的二次剩余，往下看欧拉判别准则
ll Legendre(ll a, ll p)
{
return qsm(a, (p-1)/2, p);
} 欧拉判别准则
先来回顾一下欧拉定理)，如果等于1就肯定开的了方，为-1一定开不了。所以x是否为n的二次剩余就用这个欧拉判别准则。
if(qsm(n,(p-1)/2)==p-1){
printf("No root\n");
continue;
}-1在mod p意义下为p-1。
算法流程
给出n和p
就算我们n关于p的勒让德符号为1，那么要怎样取开n的方呢？
现在是脑洞大开的时候。
找一个数a
我们随机一个数a，然后对1
先不要管为什么，后面会讲，我们现在就默默的去膜拜Cipolla的脑洞很大。
while(1){
a=rand()%p;
w=(a*a-n+p)%p;
if(qsm(w,(p-1)/2)==p-1)break;
}因为是随机找a，那么会不会要找很久。
答案：不会！
2
根据定理1，随机找a的期望为2。
建立一个复数域
说是建立，其实根本不用程序去打，说是建立复数域只是跟方便理解。
在平常学的复数域中，有一个n，这样我们就定义了一个新的域。
struct CN{
ll x,y;
CN friend operator *(CN x,CN y){
CN z;
z.x=(x.x*y.x%p+x.y*y.y%p*w%p)%p;
z.y=(x.x*y.y%p+x.y*y.x%p)%p;
return z;
}
}u,v;像正常打复数运算一样我们定义一下运算符。可以发现z.x那个地方后面是*w而不是*1，因为现在的域单位复数为ω
得出答案
我们要求的是2
真的吗？真的！但是这个答案不是由实数和虚数组成的吗？
根据拉格朗日定理，可以得出虚数处的系数一定为0。
为什么会有两个答案，比如)。
证明原理
再搞出一些定理方便说明。
定理
)
推导
问题：)
这脑洞太大了！！！！！！！！！！！！！！
代码
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 ll w,a,n,p; 5 struct CN 6 { 7 ll x,y; 8 CN operator * (const CN &a)const 9 { 10 CN z; 11 z.x=(x*a.x%p+y*a.y%p*w%p)%p; 12 z.y=(x*a.y%p+y*a.x%p)%p; 13 return z; 14 } 15 }u; 16 ll qsm(ll x,ll y) 17 { 18 ll ans=1; 19 while(y) 20 { 21 if(y&1)ans=ans*x%p; 22 x=x*x%p;y>>=1; 23 } 24 return ans; 25 } 26 CN Cqsm(CN x,ll y) 27 { 28 CN z;z.x=1;z.y=0; 29 while(y) 30 { 31 if(y&1)z=z*x; 32 x=x*x;y>>=1; 33 } 34 return z; 35 } 36 int main() 37 { 38 int T; 39 scanf("%d",&T); 40 while(T--) 41 { 42 scanf("%lld%lld",&n,&p); 43 n%=p; 44 if(p==2) 45 { 46 printf("1\n"); 47 continue; 48 } 49 if(qsm(n,(p-1)/2)==p-1){ 50 puts("No root"); 51 continue; 52 } 53 while(1) 54 { 55 a=rand()%p; 56 w=(a*a%p-n+p)%p; 57 if(qsm(w,(p-1)/2)==p-1)break; 58 } 59 u.x=a;u.y=1; 60 u=Cqsm(u,(p+1)/2); 61 ll ans1=u.x,ans2=p-ans1; 62 if(ans1>ans2)swap(ans1,ans2); 63 if(ans1==ans2){ 64 printf("%lld\n",ans1); 65 } 66 else{ 67 printf("%lld %lld\n",ans1,ans2); 68 } 69 } 70 return 0; 71 }