Cipolla算法学习小记
来源:互联网 发布:秘密知乎 编辑:程序博客网 时间:2024/06/06 17:27
简介
Cipolla算法是解决二次剩余强有力的工具,一个脑洞大开的算法。
认真看懂了,其实是一个很简单的算法,不过会感觉得出这个算法的数学家十分的机智。
基础数论储备
二次剩余
首先来看一个式子
勒让德符号
勒让德符号是判断一个数是否为p的二次剩余的一个有力工具,p一定要为奇质数。(n,p)表示n为关于p的勒让德符号。其实就是判断n是否为p的二次剩余。
看起来好像是一大段废话,勒让德只是站在巨人的肩膀上总结了一下而已。
其实勒让德还总结了一些性质,不过一般只用得到欧拉判别准则,不够勒让德符号是个积性函数可能还有点用。
但还是不知道如何判断n是否为p的二次剩余,往下看欧拉判别准则
ll Legendre(ll a, ll p) { return qsm(a, (p-1)/2, p); }
欧拉判别准则
先来回顾一下欧拉定理
因为这里的p是奇质数,所以
对
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,然后对
就是找到一个a满足
先不要管为什么,后面会讲,我们现在就默默的去膜拜Cipolla的脑洞很大。
while(1){ a=rand()%p; w=(a*a-n+p)%p; if(qsm(w,(p-1)/2)==p-1)break;}
因为是随机找a,那么会不会要找很久。
答案:不会!
根据定理1,随机找a的期望为2。
建立一个复数域
说是建立,其实根本不用程序去打,说是建立复数域只是跟方便理解。
在平常学的复数域中,有一个
我们也建立一个类似的域,因为我们要对于
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,因为现在的域单位复数为
得出答案
我们要求的是
我们现在知道了
答案=
真的吗?真的!但是这个答案不是由实数和虚数组成的吗?
根据拉格朗日定理,可以得出虚数处的系数一定为0。
CN Cqsm(CN x,ll y){\\复数的快速幂 CN z;z.x=1,z.y=0;\\注意初始化 while(y){ if(y&1)z=z*x; x=x*x; y/=2; } return z;}
w=(a*a-n+p)%p; u.x=a,u.y=1;//为什么u.y是1——在复数的统计中只用统计系数就好了 u=Cqsm(u,(p+1)/2); ll yi=u.x,er=p-u.x; if(yi>er)swap(yi,er); if(yi==er){ printf("%lld\n",yi); } else{ printf("%lld %lld\n",yi,er); }
为什么会有两个答案,比如
证明原理
再搞出一些定理方便说明。
定理
可以发现,只有当i=0或i=p的时候式子没有p的因子,所以中间有p的因子的都可以省略,所以
所以
推导
问题:
Cipolla:x≡
直接把式子转化:
所以
这脑洞太大了!!!!!!!!!!!!!!
由于本人是一个蒟蒻
对于这个算法知道的也只有这么多了。
推荐题目
现在只找到了一个较能入手的题目。
Timus Online Judge 1132 Square Root
Code
上面那道题的代码。
#include<iostream>#include<cstdio>#include<cstring>#include<cmath>#include<algorithm>#define fo(i,a,b) for(i=a;i<=b;i++)using namespace std;typedef long long ll;ll i,j,k,l,t,m,ans,w,a;ll n,p;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;ll qsm(ll x,ll y){ ll z=1; while(y){ if(y&1)z=z*x%p; x=x*x%p; y/=2; } return z;}CN Cqsm(CN x,ll y){ CN z;z.x=1,z.y=0; while(y){ if(y&1)z=z*x; x=x*x; y/=2; } return z;}int main(){ for(scanf("%lld",&t);t;t--){ scanf("%lld%lld",&n,&p); n%=p; if(p==2){ printf("1\n"); continue; } if(qsm(n,(p-1)/2)==p-1){ printf("No root\n"); continue; } while(1){ a=rand()%p; w=(a*a-n+p)%p; if(qsm(w,(p-1)/2)==p-1)break; } u.x=a,u.y=1; u=Cqsm(u,(p+1)/2); ll yi=u.x,er=p-u.x; if(yi>er)swap(yi,er); if(yi==er){ printf("%lld\n",yi); } else{ printf("%lld %lld\n",yi,er); } }}
- Cipolla算法学习小记
- 二次剩余Cipolla算法学习小记
- 二次剩余Cipolla算法 【转载a_crazy_czy】
- SIFT算法学习小记
- SIFT算法学习小记
- Sift算法学习小记
- SIFT算法学习小记
- SIFT算法学习小记
- 贪心算法 学习小记
- SIFT算法学习小记
- BM算法学习小记
- KM算法学习小记
- km算法学习小记
- 【转贴】SIFT算法学习小记
- 莫队算法学习小记
- 莫队算法学习小记
- 莫队算法学习小记
- 莫队算法学习小记
- Matlab函数大全
- 备份恢复、多租户,样样都不能少 -- 谈谈BigInsights企业管理模块的作用
- Git版本控制使用方法入门教程
- 进阶级 - Git Hub 常用指南
- C# ArrayList 的使用
- Cipolla算法学习小记
- 杭电oj 2001 计算两点间的距离
- 为什么要将数据源注入会话工厂
- AJAX学习笔记
- 关于图片分辨率
- 【Codeforces Round 362 (Div 2)C】【STL-map 最近公共祖先思想】Lorenzo Von Matterhorn 数域二叉树的路径权值变更查询
- Java第二天
- Android dex分包方案以及热补丁修复
- Spring学习篇:AOP知识整理