Cipolla算法学习小记

Posted 大奕哥&VANE

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了Cipolla算法学习小记相关的知识,希望对你有一定的参考价值。

转自:http://blog.csdn.net/doyouseeman/article/details/52033204

简介

Cipolla算法是解决二次剩余强有力的工具,一个脑洞大开的算法。 
认真看懂了,其实是一个很简单的算法,不过会感觉得出这个算法的数学家十分的机智。

基础数论储备

二次剩余

首先来看一个式子x2n(modp),我们现在给出n,要求求得x的值。如果可以求得,n为mod p的二次剩余,其实就是n在mod p意义下开的尽方。Cipolla就是一个用来求得上式的x的一个算法。

勒让德符号

勒让德符号是判断一个数是否为p的二次剩余的一个有力工具,p一定要为奇质数。(n,p)表示n为关于p的勒让德符号。其实就是判断n是否为p的二次剩余。 

(np)=???1pn,np1pnnp0pn

 

看起来好像是一大段废话,勒让德只是站在巨人的肩膀上总结了一下而已。 
其实勒让德还总结了一些性质,不过一般只用得到欧拉判别准则,不够勒让德符号是个积性函数可能还有点用。 
但还是不知道如何判断n是否为p的二次剩余,往下看欧拉判别准则

 ll Legendre(ll a, ll p)  
  {  
       return qsm(a, (p-1)/2, p);  
  } 

欧拉判别准则

先来回顾一下欧拉定理xφ(p)1(modp) 
因为这里的p是奇质数,所以xp11(modp) 
xp1进行开方操作,在虚数域中xp12±1(modp),如果等于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,然后对a2n进行开方操作(就是计算他勒让德符号的值),直到他们的勒让德符号为-1为止(就是开不了方为止)。 
就是找到一个a满足(a2n)p12=1 
先不要管为什么,后面会讲,我们现在就默默的去膜拜Cipolla的脑洞很大。

while(1){
    a=rand()%p;
    w=(a*a-n+p)%p;
    if(qsm(w,(p-1)/2)==p-1)break;
}

因为是随机找a,那么会不会要找很久。 
答案:不会! 
?1x2n(modp)p12n使 
1x2n(modp),如果存在不同的数u,v,使他们带入x后可以使方程有解,那么很显然满足u2v2|p,所以满足(u+v)(uv)|p,因为 
u2v2|p所以p不可能是(u-v)的倍数,所以(u+v)|p,那么这样的数对在p中存在的数量为p12 
根据定理1,随机找a的期望为2。

建立一个复数域

说是建立,其实根本不用程序去打,说是建立复数域只是跟方便理解。 
在平常学的复数域中,有一个i,满足i2=1。 
我们也建立一个类似的域,因为我们要对于a2n开方,而a2n有不是p的二次剩余,所以我们定义ω=a2n−−−−−√。那么现在的ω也像i一样,满足ω2=a2n,这样我们就定义了一个新的域。

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=a2n,而不像正常复数的i2=1。在这个域中的表示方式类似正常的复数:a+bω

得出答案

我们要求的是x2n(modp),x的值 
我们现在知道了aω之后,就能得出答案了。 
答案=(a+ω)p+12 
真的吗?真的!但是这个答案不是由实数和虚数组成的吗? 
根据拉格朗日定理,可以得出虚数处的系数一定为0。

 

为什么会有两个答案,比如4√=±2x2(px)2(modp)很显然,因为把后面拆开x2x22px+p2(modp),把p都约去,所以x2(px)2(modp)

证明原理

再搞出一些定理方便说明。

定理

?2(a+b)pap+bp(modp) 
2: 
(a+b)ppi=0Cipapibi(modp) 
pi=0p!(pi)!i!apibi(modp) 
可以发现,只有当i=0或i=p的时候式子没有p的因子,所以中间有p的因子的都可以省略,所以(a+b)pap+bp(modp) 
?3ωpω(modp) 
3ωp 
=ωp1ω 
=(ω2)p12ω 
=(a2n)p12ω——因为ω2=a2n 
=ω——因为满足(a2n)p12=1 
?4apa(modp) 
3ap 
ap11(modp) 
所以apaap1a(modp)

推导

问题:x2n(modp)求解x 
Cipolla:x≡(a+ω)p+12(modp)的实数 
直接把式子转化: 
(a+ω)p+12(modp) 
((a+ω)p+1)12(modp) 
((a+ω)p(a+ω))12(modp) 
((ap+ωp)(a+ω))12(modp)根据定理2 
((aω)(a+ω))12(modp)根据定理3和定理4 
(a2ω2)12(modp)根据定理3和定理4 
(a2(a2n))12(modp)满足ω2=a2n 
n12(modp) 
所以(a+ω)p+12n12n√(modp) 
这脑洞太大了!!!!!!!!!!!!!!

代码

 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 }

 

以上是关于Cipolla算法学习小记的主要内容,如果未能解决你的问题,请参考以下文章

二次剩余 Cipolla算法

2018.07.29(搜索)学习DFS算法小记

二次剩余-Cipolla

二次剩余(Cipolla)

二次剩余(Cipolla)

KM算法学习小记: