如何用分布式Pollard-Rho法对椭圆曲线离散对数问题（ECDLP）进行攻击（上）

最新推荐文章于 2022-09-22 01:40:01 发布

BUAA_Alchemist

最新推荐文章于 2022-09-22 01:40:01 发布

阅读量4.6k

点赞数 6

分类专栏：密码学文章标签： ECDLP Pollard Rho 分布式

本文链接：https://blog.csdn.net/BUAA_Alchemist/article/details/94437596

版权

密码学专栏收录该内容

4 篇文章 3 订阅

订阅专栏

上个学期，竞赛因为种种原因没有继续进行而搁浅。这里仅仅分享一下临时成果。
椭圆曲线离散对数问题赛题介绍如下
在这里插入图片描述

任何一个熟悉密码学的人基本上都会了解椭圆曲线离散对数问题。它比较普通素数域离散对数问题更具有难解性。因此相当多的公钥密码体制都构建在这个问题的基础上，如EC-Elgamal系统，以及数字签名、密钥交换协议等。再次不多介绍。

接下来给出上述数据文档的内容

		椭圆曲线参数形式为
		E/F_p: y^2=x^3+ax+b,n:=#E(F_p ),P∈E(F_p),r:=order(P),R∈<P>.
		求解k,1≤k≤r,s.t.R=[k]P.

第一类曲线参数：
(1)  p:= 211108170305887; a:=0;b:=7;n:=p;r:=p;
P:=( 47815642535808, 116240163507508);
R:=( 77983503452527, 143728424564583);

(2) p:= 13835058061724614657; a:=0;b:=20;n:=p;r:=p;
P:=( 616859655854051956, 12065166484289278801);
R:=( 5170466145333976578, 7139090565738339416);

(3) p:=906694364778591846139117; a:=0; b:=2;n:=p;r:=p;
P:=(475325122433864976165476 ,666857317692667708141096 );
R:=(519814743429987512024682, 392414632044199857512746);

(4) p:= 59421121885714719481295537269; a:=0;b:=10;n:=p;r:=p;
P:=(17547290159953212409742744311, 23276625757393135830641872446);
R:=(22782721588122522786532109807, 29566916200346945584248955766);

(5) p:=3894222643901127098494944603540019; a:=0; b:=2;n:=p;r:=p;
P:=(3474281736844926688615305014567004, 3343311742974261537268420184037101);
R:=(46006664812763786791056435590121, 1631023347800240287678172773820495);

(6) p:= 255211775190703851000955237173238443091; 
a:=0;b:=32;n:=p;r:=p;
P:=(82054120567654459070422632716611948091,
208358019881453692582450632924824868211);
R:=(54625405255845450735684923869953661538,
82751760415032967427013207973543496169);

(7) p:=16725558898897967356385545845388318567564081; 
a:=0; b:=39;n:=p;r:=p;
P:=(15866255640827385375149316462253403735143979, 13637900016555731147592334085022810288787804);
R:=(1641844846280313158776293444944029290477363, 1916339757694211104728115095781816602417918);

(8) p:= 1096126227998177188652856107362412783873814431647; 
a:=0;b:=5;n:=p;r:=p;
P:= (83173790821618364013911485269418053634749973470, 331912486564013055335956381322549180967694844964);
R:=(163008382252281273947629676562612579052560281605, 50645003441262331054159546612247657430098333396);

第二类曲线参数：
(9) p:= 140737488356467; a:=1;b:=0;
n:= 140737488356468; r:= 35184372089117;
P:=( 59477512413747, 79851403980273);
R:=( 125962305399026, 124644987166940);

(10) p:= 9223372036854782251; a:=1;b:=0;
n:= 9223372036854782252; r:= 2305843009213695563;
P:=(8930887567779448763 , 890632237231967440);
R:=( 4707122633993752935 , 3224323188778636920);

(11) p:=604462909807314587364667;a:=1;b:=0;
n:=604462909807314587364668; r:=151115727451828646841167;
P:=(587411173575122535454189, 184243119926212298785598);
R:=(539012794797204313781763, 513627008874848913042314);

(12) p:= 39614081257132168796771986051; a:=1;b:=0;
n:= 39614081257132168796771986052;
r:= 9903520314283042199192996513;
P:=( 17135037968192446511660916347 , 15756828316532436197954290560);
R:=( 14307140364976860571933505517 , 31289859728052046761658770776);

(13) p:=2596148429267413814265248164627363;a:=1;b:=0;
n:=2596148429267413814265248164627364;
r:=649037107316853453566312041156841;
P:=(622497953272208887929891881018762, 578415484635633376407094129167605);
R:=(1641315995010304174750922058961802, 498062680065145119027669812682594);

(14) p:= 170141183460469231731687303715884123283; a:=1;b:=0;
n:= 170141183460469231731687303715884123284;
r:= 42535295865117307932921825928971030821;
P:=(67561795560749592594257348250381095281 ,132967032564783335773794402620839168397);
R:=( 81767392472072972289932811718039895097 , 74483734086357842747707740743879456218);

(15) p:=11150372599265311570767859136324180753031283;a:=1;b:=0;
n:=11150372599265311570767859136324180753031284;
r:=2787593149816327892691964784081045188257821;
P:=(10280081406291279076050813261615407334758360, 10524832460733629896529554181939223304289725);
R:=(1653939556554266819965991031125172966638801, 8370614898071115688249737076845551200138218);

(16) p:= 730750818665451459101842416358141509827966272147; a:=1;b:=0;
n:= 730750818665451459101842416358141509827966272148;
r:= 182687704666362864775460604089535377456991568037;
P:=(221903508784709687556506235376523499020990986034 , 287854209568598432667871196372312232614552000893);
R:=( 219270464789952726868617287704232288912801490505, 578348937880270477205296252062048509994911317888);

第三类曲线参数：
(17) p:= 140737488355333; a:=-3;b:=234;
n:= 140737484527007;r:=n;
P:= (118344265104828, 25754556069705);
R:=( 8594518695631, 14966619423525);

(18) p:= 9223372036854775837; a:=-3;b:=37;
n:= 9223372038068412403;r:=n;
P:=(7220661838117791356 ,6477969291505257777);
R:=( 7819726923954549567 ,6266868167000835108);

(19) p:=604462909807314587353111;a:=-3;b:=95;
n:= 604462909807750541909849;r:=n;
P:=(428432040100075198254744, 95025782588400118295756);
R:=(472138605558837378507194, 138148835226095728614736);

(20) p:= 39614081257132168796771975177; a:=-3;b:=562;
n:= 39614081257132147751873462213;r:=n;
P:=(39220344157117715096559716859 , 3087260393566610895498596606);
R:=( 13160703755766149325136222951 , 6533567029273948676370251736);

(21) p:= 2596148429267413814265248164610099; a:=-3;b:=171;
n:= 2596148429267413788194246433592681;r:=n;
P:=(177627746966506201527866528484813 , 1735330667419561032038559795836927);
R:=( 1575178132453901115471501570466697, 538238983595050290992251454098891);

(22) p:= 170141183460469231731687303715884105757; a:=-3;b:=33;
n:= 170141183460469231728561996679834270597; r:=n;
P:=( 96152442714630401692673834213524521597 , 58626340067602219522071225434282008266);
R:=( 47795192678737921133501161460114957836 , 166078919686417913546029950371468891352);

在比赛文件刚刚出来时，我便开始着手尝试，作为一名ICPC选手，首先会想到接触过的大步小步（BSGS）法，这个可谓是解决离散对数问题的经典算法。对于三个文件的首个数据都在C++的long long范围之内，因此用大步小步法是可以尝试的。

注：以下关于椭圆曲线的运算，我会用乘法代表普遍而言的加法，幂次代表普遍而言的数乘，目的在于和“群”的定义更为符合

BSGS法的内容很简单，要求解 $P^x=R,\ x \in(0,n)$ ，我们可以设定一个数m，将 $n$ 均分成 $\lceil\frac{n}{m}\rceil$ 份，先求出 $P^{1},P^{2}…P^{m}$ ，将这m个值，储存于一个映射表（map）中，键为点，值为幂次，接着我们开始不断求出 $Q=R*P^{-\lceil\frac{n}{m}\rceil},R*P^{-\lceil\frac{n}{m}\rceil*2}…$ 每一次尝试将点Q在映射表中查找，倘若查找到，那么有 $R*P^{-{\lceil\frac{n}{m}\rceil}*i}=P^{j}$ ，于是有 $R=P^{{\lceil\frac{n}{m}\rceil}*i+j}$ ，也就是得到解。
BSGS法的本质是分块，或者比起单纯的暴力枚举，是在用空间换时间，理想状态下，空间复杂度和时间复杂度均可以做到 $O (s q r t (n))$ 级别（这里忽略映射表的修改与查找所需的复杂度）。但是注意到一般的存储设备是无法做到储存如此量级的数据的，因此真正实现BSGS法不现实，我们只可能把m设置的远小于 $s q r t (n)$ 用少量空间来换取时间，而这对于n很大的情况，效率依旧是难以被我们所接受的。
这里给出我的第一份代码，可以解出每一组的首个数据。

//第一个版本，采用分块+map的方式，作为尝试性求解
#include<cstdio>
#include<algorithm>
#include<map>
#define mo 211108170305887LL //模数
#define a 0LL
#define b 7LL
#define x first
#define y second
using namespace std;
using LL=long long;
using Point=pair<LL,LL>; //定义ECC上的点

Point P={47815642535808LL, 116240163507508LL};
Point R={77983503452527LL, 143728424564583LL};
Point inv;
LL ans,block,sz;
map<Point,LL> M;

LL mul(LL p,LL q) //icpc 常用的快速乘模技巧
{
	p=(p%mo+mo)%mo,q=(q%mo+mo)%mo;
	return ((p*q-(LL)((long double)p/mo*q+1e-6)*mo)%mo+mo)%mo;
}

LL quick_power(LL p, LL q)  //快速幂
{
	LL res=1LL,base=p;
	while(q)
	{
		if(q&1LL)
			res=mul(res,base);
		base=mul(base,base);
		q>>=1;
	}
	return res;
}

Point Point_square(Point Q)  //点的平方
{
	Point P;
	LL tmp=mul((3LL*mul(Q.x,Q.x)%mo+a)%mo,quick_power(2LL*Q.y%mo,mo-2LL));
	P.x=(mul(tmp,tmp)-2LL*Q.x%mo+mo)%mo;
	P.y=(mul(tmp,(Q.x-P.x+mo)%mo)-Q.y+mo)%mo;
	return P;
}

Point Point_mul(Point P, Point Q)  //点的乘积
{
	if(P==Q)
		return Point_square(P);
	Point res;
	LL tmp=mul((Q.y-P.y+mo)%mo,quick_power((Q.x-P.x+mo)%mo,mo-2LL));
	res.x=(mul(tmp,tmp)-Q.x+mo-P.x+mo)%mo;
	res.y=(mul(tmp,(P.x-res.x+mo)%mo)-P.y+mo)%mo;
	return res;
}

Point Point_quick_power(Point p, LL q)  //点的快速幂
{
	Point res=p,base=p;
	--q;
	while(q)
	{
		if(q&1LL)
			res=Point_mul(res,base);
		base=Point_mul(base,base);
		q>>=1LL;
	}
	return res;
}

int main()
{
	Point PP=P,QQ=R;
	block=LL(sqrtl(mo))*10LL;  //定义分块个数，乘以10是一个内存限制的体现
	sz=(mo-1)/block+1;  //求出块的大小
	inv=Point_quick_power(P,sz);  //P^-sz
	inv.y=mo-inv.y;
	for(int i=1;i<=sz;i++)
	{
		M[PP]=i;
		PP=Point_mul(PP,P);
	}
	printf("%d",M.size());
	for(int i=0;i<=block;i++)
	{
		if(i%(block/100)==0)
			printf("%d\n",i);
		if(M.count(QQ))  //找到碰撞，求出答案
		{
			printf("%lld",i*sz+M[QQ]);
			break;
		}
		QQ=Point_mul(QQ,inv);  //求出新的QQ
	}
	return 0;
}

运行效果

139189752582973

Process exited after 420.6 seconds with return value 0 请按任意键继续. . .

（另：一开始给的数据文件第三组的第一个数据是有问题的，还是我最早发现跑不出来，官方才修正过来）
通过BSGS法解决了小数据是一个很好的开始。但接下来就不那么容易了，原本想继续沿用BSGS法，看看能不能求解出第二组数据。当时因为long long不够用，套用了自己写的大数板子，发现由于大量的取模操作，显得无比低效。另外，根据效率估计，即使大数的运算速度和long long类型一致，可能BSGS法也需要以日计去跑第二组数据，这就有些尴尬，毕竟之后的天文数据还有很多呢。
先解决大数问题，大数是解决问题的一个关键，算法相同情况下，好的大数基础上只需要一天那么可能坏的大数运算就需要一个月。这不是一个ICPC上可以忽略的常数问题。
现有的大数运算库（基于C++）有很多，有号称最快的gmp，密码学相关的库miracl，以及openssl和boost等库中也有大数运算。最后我选择了miracl，原因在于：

带有椭圆曲线相关的基础运算，会比自己基于大数实现更高效
大量的已有函数，用起来比较方便
安装比较容易

个人认为，如果不是miracl能提供椭圆曲线相关运算的话，我更倾向于选择采用Java或者Python，而非费尽心思去安装一个大数库。事实上，单纯从大数加减乘除运算速度来讲，miracl甚至略逊于Python（所以不要单单讲Python是个低效的语言，同样的功能下Python可能远远超乎你想象）。
另外，“安装比较容易”其实也是相比较于其它库而言。事实上因为C++没有跨平台性，装任何一个扩展库都是极其麻烦的一件事。从github上下载完毕之后，我倒腾了也不知道有多久，才终于在VS2017上能成功使用这个库。具体过程不再详细叙述。
下面给出基于Miracl库的BSGS法的代码，注：这里只是练习使用大数big.h，没有使用ecn.h中的椭圆曲线内容

//和ecdlpsolver1和2的方法完全相同，但是采用了miracl库的big类进行大整数运算（主要是为了学习使用方法）
//第一个正式的求解器，能求解三类中的第一组。
//deathmask,2019.1.23

#include<iostream>
#include<algorithm>
#include<map>
//#include<bitset>
#include "big.h"
#define LIMIT 24
//#define hashmod 53000007
using namespace std;

Miracl precision(LIMIT);

using Point = pair<Big, Big>;

Point P = {"118344265104828", "25754556069705"};
Point R = {"8594518695631", "14966619423525"};
//Point R = { "45389087271562","48223143538350" };
Big mo = "140737488355333";
Big a = "140737488355330";
Big r = "140737484527007";
Point inv,PP,RR;
Big ans, block, per,cnt;
map<Point, int> M;

int sz = 4000000;

inline Big mul(Big x, Big y)
{
	return modmult(x, y, mo);
}

inline Big div(Big x, Big y)
{
	return moddiv(x, y, mo);
}

Point Pt_square(Point Q)
{
	Point P;
	Big tmp = div(mul(Q.first, Q.first) * 3  + a, Q.second <<1);
	P.first = (mul(tmp, tmp) + (mo-Q.first<<1) ) % mo;
	P.second = (mul(tmp, Q.first - P.first + mo) - Q.second + mo) % mo;
	return P;
}

Point Pt_mul(Point P, Point Q)
{
	if (P == Q)
		return Pt_square(P);
	Point res;
	Big tmp = div(Q.second - P.second + mo, Q.first - P.first + mo);
	res.first = (mul(tmp, tmp) - Q.first + mo - P.first + mo) % mo;
	res.second =( mul(tmp, P.first - res.first + mo) - P.second + mo)%mo;
	return res;
}

Point Pt_power(Point p, Big q)
{
	Point res = p, base = p;
	--q;
	while (!q.iszero())
	{
		if (q%2)
			res = Pt_mul(res, base);
		base = Pt_mul(base, base);
		q = q >> 1;
	}
	return res;
}

/*inline int Pt_hash(const Point &x)
{
	return (pow(x.first, 5, hashmod) + pow(x.second, 3, hashmod)) % hashmod;
}*/

int main()
{
	ios::sync_with_stdio(false);
//	cout << Pt_power(P, "58891866538906").second << endl;
	PP = P, RR = R;
	block = (r - 1) / sz + 1;
	inv = Pt_power(P, sz);
	inv.second = mo - inv.second;
	per = sz / 100;
	for (int i = 1; i <= sz; i++)
	{
		M[PP] = i;
		PP = Pt_mul(PP, P);
		if (++cnt == per)
		{
			cout << "loaded " << i << endl;
			cnt = 0;
		}
	}
	cnt = 0;
	cout << "load_complete " << endl;
	per = block / 100;
	for (Big i = 0; i < block; ++i)
	{
		if (per == ++cnt)
		{
			cout << "\tstep:" << i << endl;
			//cout << RR.first << ' ' << RR.second << endl;
			cnt = 0;
		}
		if (M.count(RR))
		{
			cout << "The answer is " << M[RR] + i*sz << endl;
			break;
		}
		RR = Pt_mul(RR, inv);
	}
	return 0;
}

大数的问题暂时解决了。接着，我学习了一种Pollard-Rho法。有必要说明的是，Pollard-Rho是一类方法，比如在大数分解质因数问题方面，也有一套Pollard-Rho法。细节不叙述，原理可以参考 https://www.cnblogs.com/dalt/p/8437119.html ，代码如下

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<vector>
#include<ctime>
using namespace std;
using LL=long long;

vector<LL> f;

inline LL B_rand()
{
	return (LL)rand()*rand()*rand()*rand();
}

LL quick_mul(LL x, LL y, LL mo)
{
	//x,y过大，long double可能丢失个位，使用前最好能保证x,y比mo小，否则可以先对mo取模 
	//这里能保证。 
	return (x*y-(LL)((long double)x/mo*y)*mo+mo)%mo;
}

LL quick_power(LL x, LL y, LL mo)
{
	LL res=1,base=x;
	while(y)
	{
		if(y&1)
			res=quick_mul(res,base,mo);
		base=quick_mul(base,base,mo);
		y>>=1;
	}
	return res;
}

bool miller_rabin(LL n, int s)  //miller_rabin判断大质数
{
	if(n==2)
		return true;
	if(n==1||!(n&1))
		return false;
	int t=0;
	LL a,x,y,u=n-1;
	while(!(u&1))
		t++,u>>=1;
	for(int i=1;i<=s;i++)
	{
		a=B_rand()%(n-1)+1;
		x=quick_power(a,u,n);
		for(int j=0;j<t;j++)
		{
			y=quick_mul(x,x,n);
			if(y==1&&x!=1&&x!=n-1)
				return false;
			x=y;
		}
		if(x!=1)
			return false;
	}
	return true;
}

LL Pollard_Rho(LL n, LL c)
{
	LL i=1,k=2,x=B_rand(),y=x;
	for(;;)
	{
		i++;
		x=(quick_mul(x,x,n)+c)%n;  //更新随机量x，迭代步
		LL p=__gcd((y-x+n)%n,n);
		if(p!=1&&p!=n)  //y-x mod m=0 and y!=x
			return p;
		if(y==x)  //这一步是一个建议步，如果y=x，属于没有作用的碰撞
		//继续找在这个环上找有效碰撞，有可能，但概率低，步数大，还有概率死循环，不如从头开始
			return n;
		if(i==k)  //到k步时更新k值
			y=x,k<<=1;
	}
}

void find(LL n, LL c)  //分解n，c是一个随机步进值
{
	if(n==1)
		return ;
	if(miller_rabin(n,30))
	{
		f.push_back(n);
		return ;
	}
	LL p=n,k=c;
	while(p>=n)  //尝试分解
		p=Pollard_Rho(p,c),c=(c+1)%n;  //修改c
	find(p,k);
	find(n/p,k);
}

LL n;

int main()
{
	srand(time(0));
	while(scanf("%lld",&n)==1)
	{
		f.clear();
		find(n,max(B_rand()%n,1LL));
		sort(f.begin(),f.end());
		for(LL i:f)
			printf("%lld ",i);
		puts("");
	}
	return 0;
}

Rho法的复杂度有多种证明。但最直观的理解还是上面网址中所给出，偏向于图论的方式。假设一个有m个点的图，每个点都有一个出度，这样随机形成的图中，无论是子基环图（Rho形）还是环本身，包含点数的期望都会是 $\sqrt{m}$ 级别的。
倘若要分解的数是n，那么将迭代产生的数列进行连接可以得到一个基环图，不妨称之为大Rho。考虑其中最小的那个因子m，必然不会超过 $\sqrt{n}$ ，把n个点中关于m同余的点合并成一个点，那么每次迭代步形成的图的大小就不会超过 $\sqrt{\sqrt{n}}$ ，可以称之为小Rho，我们也可以发现，大Rho的环形可以看作由多个小Rho的环剪开顺次连接形成，我们现在要做的，就是在大Rho上找出点x和y，使得x和y不同，但在小Rho图中对应的点是相同的（碰撞）。
那么实际上要解决的就是找碰撞值，或者说在小Rho图上做一个找环操作。有相当多的方法，比如最简单的记录表法（需要内存），Floyd转圈（O(1)内存），还有效率更优越的Brent法等。以上的程序，采用的便是典型的Brent法。
回到离散对数问题。Pollard-Rho法的解决方案如下。

我们令G表示椭圆曲线上的点集合，然后我们将G随机分成三个大小相近的部分 $S_1,S_2,S_3$ （可以用哈希函数的方法）
然后我们设定一个用于生成随机的步数的迭代函数f：
$T_{i+1}=f\left(T_{i}\right)=\left\{\begin{array}{l}{R*T_{i}, T_{i} \in S_{1}} \\ {T_{i}^2, T_{i} \in S_{2}} \\ {P*T_{i}, T_{i} \in S_{3}}\end{array}\right.$
令 $T_i= P^{A_i}*R^{B_i}$ ，则有：
$A_{i+1}=\left\{\begin{array}{l}{A_{i}, T_{i} \in S_{1}} \\ {2A_{i} \bmod r, T_{i} \in S_{2}} \\ {(A_{i}+1) \bmod r, T_{i} \in S_{3}}\end{array}\right.$
$B_{i+1}=\left\{\begin{array}{l}{(B_{i}+1) \bmod r, T_{i} \in S_{1}} \\ {2B_{i} \bmod r, T_{i} \in S_{2}} \\ {B_{i}, T_{i} \in S_{3}}\end{array}\right.$
然后初始化参数，取 $T_0=P，A_0=1，B_0=0$ ，然后不断生成配对 $T_i,T_{2i})$ ，直到我们找到了某个m，有 $T_m=T_{2m}$ ，则有：
$\left\{\begin{array}{l}{T_{m}=P^{A_{m}} *R^{B_{m}}} \\ {T_{2 m}=P^{A_{2 m}}*R^{B_{2 m}} }\end{array}\right.$
于是可以求出
$Ans=\frac{A_{2 m}-A_{m}}{B_{m}-B_{2 m}}(\bmod\ r)$

表面上看，和之前的分解大数的Rho法差别很大，但从图论的角度理解的话，确实是一回事。我们将用迭代函数F为G中的点连边（近乎随机），形成多个基环，期望长度为 $\sqrt r$ 级别。接着我们期望要找到一个有效对撞。
注：此处叙述采用的找环方法为Floyd法，也可以用其他方法。
下面给出一份我的首个用Rho法求解ECDLP的代码

//Aphrodite: pollard-rho法，效率低，待改进,能处理三个类的第一个数据
#include<cstdio>
#include<ctime>
#include "ecn.h"
using namespace std;

Miracl precision(50, 0);  //初始化大数系统
miracl *mip = &precision;
Big mo = "140737488355333";
Big a = "-3";
Big b = "234";
Big r = "140737484527007";
/*
Big mo = 47;
Big a = 34;
Big b = 10;
Big r = 41;*/
Big A1 = 1, B1 = 0, A2 = 1, B2 = 0;
ECn P, R, Fm, F2m;
int c0, c1, c2;

//char buf[55];
Big x;

inline int Block(ECn &num)  //判断在那个S区
{
	/*register unsigned int h=0;
	mip->IOBASE = 16;
	buf << x;
	for (register int i = 0; buf[i]; i++)
	h = h * 31 + buf[i];
	buf << y;
	for (register int i = 0; buf[i]; i++)
	h = h * 37 + buf[i];
	mip->IOBASE = 10;
	return h%3;*/
	num.get(x);
	return x[0] % 3;  //用x坐标最低位的值模3的余数来确定，速度比较快
}

inline void F(Big &A, Big &B, ECn &Tmp)
{
	switch (Block(Tmp))  //迭代函数
	{
	case 0:
		++B;
		//	c0++;
		Tmp += R;
		break;
	case 1:
		A <<= 1;
		B <<= 1;
		//	c1++;
		Tmp += Tmp;
		break;
	default:
		++A;
		//	c2++;
		Tmp += P;
	}
	if (A >= r)  //取模
		A -= r;
	if (B >= r)  //取模
		B -= r;
}

Big pollard_rho()
{
	register int cnt = 0;
	Fm = F2m = P;
	do
	{
		F(A1, B1, Fm);  //Floyd法
		F(A2, B2, F2m);
		F(A2, B2, F2m);
		if (++cnt % 50000 == 0)  //中间信息
			cout << A1 << ' ' << B1 << ' ' << A2 << ' ' << B2 << ' ' << Fm << ' ' << F2m << endl;
	} while (Fm != F2m);
	B2 %= r, B1 %= r;  //貌似可以去掉？
	return moddiv(A2 - A1 + r, B1 - B2 + r, r);  //求出答案，这里忘了判断B1==B2要重新开始的情况了，小疏忽
}

int main()
{
	ios::sync_with_stdio(false);
	long long sta = clock();
	ecurve(a, b, mo, MR_PROJECTIVE);
	P = ECn("118344265104828", "25754556069705");
	R = ECn("8594518695631", "14966619423525");
	cout << P << ' ' << R << endl;
	cout << "Ans:" << pollard_rho() << endl;
	cout << (double)(clock() - sta) / CLOCKS_PER_SEC << endl;
	return 0;
}

以及后来用Brent法的版本

//Aldebaran:pollard-rho法，有增进
#include<cstdio>
#include<ctime>
#include "ecn.h"
using namespace std;

const int sz = 15;
const int csz = 20;

Miracl precision(20, 0);  //初始化大数系统，20单元，每单元极限储存
miracl *mip = &precision;
Big mo = "140737488355333";
Big a = "-3";
Big b = "234";
Big r = "140737484527007";
Big A1, B1, A2, B2,x;
ECn P, R, Fm, F2m;
ECn stp[sz];
Big da[sz], db[sz];

inline int Block(ECn &num)
{
	num.get(x);
	return x[0] % csz;
}

//这里对F函数进行了修改：3/4的概率进行线性迭代，迭代步数是随机化的，1/4的概率进行平方迭代
//修改的启发部分来自于论文。虽然没有理解原因，但效果似乎确实微微胜于前者

inline void F(Big &A, Big &B, ECn &Tmp)
{
	register int i = Block(Tmp);
	if (i < sz)
	{
		A += da[i]; 
		B += db[i];
		Tmp += stp[i];
	}
	else
	{
		A <<= 1;
		B <<= 1;
		Tmp += Tmp;
	}
	if (A >= r)
		A -= r;
	if (B >= r)
		B -= r;
}

Big pollard_rho()  //改用Brent法，效率提高
{
	register long long cnt = 0,limit=2,taken=1;
	Fm = F2m = mul(A1,P,B1,R);
	F(A2, B2, F2m);
	while(Fm!=F2m)
	{
		if (taken == limit)
		{
			taken = 0;
			limit <<= 1;
			A1 = A2;
			B1 = B2;
			Fm = F2m;
		}
		++taken;
		F(A2, B2, F2m);
		if(++cnt%100000==0)
			cout <<cnt/100000<<':'<< A1 << ' ' << B1 << ' ' << A2 << ' ' << B2 << ' ' << Fm << ' ' << F2m << endl; //中间信息
	}
	B2 %= r, B1 %= r;  //这里同样忘了判断B1是不是等于B2，疏忽
	return moddiv(A2 - A1 + r, B1 - B2 + r, r);
}

int main()
{
	ios::sync_with_stdio(false);
	long long sta = clock();
	ecurve(a, b, mo, MR_PROJECTIVE);
	P = ECn("118344265104828", "25754556069705");
	R = ECn("8594518695631", "14966619423525");
	for (int i = 0; i < sz; i++)
	{
		da[i] = rand(r);
		db[i] = rand(r);
		stp[i] = mul(da[i], P, db[i], R);
		cout << stp[i] << endl;
	}
	A1 = A2 = rand(r);
	B1 = B2 = rand(r);
	cout << P << ' ' << R << endl;
	cout << "Ans:" << pollard_rho() << endl;
	cout << (double)(clock() - sta) / CLOCKS_PER_SEC << endl;
	return 0;
}

对比而言，Pollard-Rho法已经可以拥有和BSGS相当的时间效率 $O(\sqrt n)$ ，而不必需要很可怕的空间复杂度。
接下来的这个版本进一步优化。从Brent法改进到被认为是最有效的Quisquater-Delescaille法（以下简称QD法），它的效率会胜于Brent法，而且非常适合多机运行（分布式）。首先，设置一个“桩点判断函数”，这个函数必须高效，作用是筛选一小部分具备某些特征的点进行记录，这些点的选取会呈现均匀的特点。如果在迭代的过程中不断记录桩点，那么如果碰到一个新的桩点，每次在表中查询，如果表中能查到，那么可能就是一个有效的碰撞了。
无论哪种找环的方法，在迭代完一圈之后都要额外迭代一段，才能发现确实成环。然而QD法的额外这一段的期望长度，是比Floyd法和Brent法都要更优的，从图的角度也很好理解这一点。
对于桩点出现的频率需要注意，频率太低意味着需要迭代更多的期望长度，但频率太高的话对期望长度的缩小是非常有限的，而且额外的代价会很大（很大的映射表，不断地查询，修改）。
可以看出，QD法基本上是记录表法和其他 $O (1)$ 空间方法的折衷，是在用合适的空间来换取效率的进步。
另一个优化在于我对Miracl库理解的加深。用VS的profile功能给程序测试的时候，发现ECn的get函数开销不小。开始很困惑，觉得get只是返回一个已有坐标呀。后来改用getx函数，更奇怪的事情发生了，得到的结果会有时候不是正确的值？这是库的bug吗？仔细阅读文档我才发现，出于运算效率的原因，ECn中点的坐标，有时候并不是真实的坐标，getx得到的就是“假坐标”，而get函数会先调用一个normalise函数，对坐标进行“规则化”，之后才返回求出的坐标结果。得知这一点就很棒了，我们只要在合适的时间调用normalise函数，那么没有运算过的话，调用getx函数就可以返回x坐标了。这个改进可以说是比较有效的提升。

//Shura：pollard-rho法，继续增进从上一版本的brent法改为了QD法
#include<cstdio>
#include<ctime>
#include<map>
#include "ecn.h"
using namespace std;

const int sz = 15;
const int csz = 20;
const int check_mod = 20003;  //桩点的频率

Miracl precision(20, 0);
miracl *mip = &precision;
Big mo = "13835058061724614657";
Big a = "0";
Big b = "20";
Big r = mo;
Big A1, B1, A2, B2, x;
Big x2, y2, x3, y3;
ECn P, R, Fm;
ECn stp[sz];
Big da[sz], db[sz];

struct cmp_key  //比较函数，用于map
{
	inline bool operator () (const ECn &p, const ECn &q) const
	{
		p.getxy(x2, y2);
		q.getxy(x3, y3);
		return x2 < x3 || x2 == x3&&y2 < y3;
	}
};

map<ECn, pair<Big, Big>, cmp_key> M;

inline int Block(ECn &num)
{
	num.getx(x);
	return x[0] % csz;
}

inline void F(Big &A, Big &B, ECn &Tmp)
{
	register int i = Block(Tmp);
	if (i < sz)
	{
		A += da[i];
		B += db[i];
		Tmp += stp[i];
	}
	else
	{
		A <<= 1;
		B <<= 1;
		Tmp += Tmp;
	}
	if (A >= r)
		A -= r;
	if (B >= r)
		B -= r;
}

inline bool check(ECn &Fm)
{
	Fm.getx(x);
	return x[0] % check_mod == 0;
}

Big pollard_rho()
{
	register long long cnt = 0, msz=0;
	for (;; F(A1, B1, Fm))
	{
		normalise(Fm);
		if (check(Fm))  //检查是不是桩点
			if (M.count(Fm))  //如果已存在
			{
				A2 = M[Fm].first, B2 = M[Fm].second;
				if(B1!=B2)  //这里判断了B1是不是等于B2，如果等的话会选择从头来过
				//会不会有更好的解决方式呢？
					break;
				A1 = rand(r);
				B1 = rand(r);
				Fm = mul(A1, P, B1, R);
			}
			else
				M[Fm] = { A1,B1 },++msz;  //插入
		if (msz > 2000000)  //map的大小太大了
		{
			cout << "WARNING!!!" << endl;
			exit(0);
		}
		if (++cnt % 100000 == 0)
			cout << cnt / 100000 << ':' << A1 << ' ' << B1 << ' ' << Fm <<' '<< M.size()<< endl;
	}
	return moddiv(A2 - A1 + r, B1 - B2 + r, r);
}

int main()
{
	ios::sync_with_stdio(false);
	long long sta = clock();
	ecurve(a, b, mo, MR_PROJECTIVE);
	P = ECn("616859655854051956", "12065166484289278801");
	R = ECn("5170466145333976578", "7139090565738339416");
	for (int i = 0; i < sz; i++)
	{
		da[i] = rand(r);
		db[i] = rand(r);
		stp[i] = mul(da[i], P, db[i], R);
		cout << stp[i] << endl;
	}
	A1 = rand(r);
	B1 = rand(r);
	Fm = mul(A1, P, B1, R);
	cout << P << ' ' << R << endl;
	cout << "Ans:" << pollard_rho() << endl;
	cout << (double)(clock() - sta) / CLOCKS_PER_SEC << endl;
	return 0;
}

至此，单机的Pollard-Rho算法基本没有很大的优化空间了。然而即便如此，跑每组的第二组数据，运气不佳仍然需要大概五六个小时。一方面我找方法对底层运算进行优化，另一方面开始着手设计分布式Pollard-Rho法。篇幅留给（下）

BUAA_Alchemist

关注

6
点赞
踩
40

收藏

觉得还不错? 一键收藏
2
评论
如何用分布式Pollard-Rho法对椭圆曲线离散对数问题（ECDLP）进行攻击（上）

上个学期，竞赛因为种种原因没有继续进行而搁浅。这里仅仅分享一下临时成果。椭圆曲线离散对数问题赛题介绍如下任何一个熟悉密码学的人基本上都会了解椭圆曲线离散对数问题。它比较普通素数域离散对数问题更具有难解性。因此相当多的公钥密码体制都构建在这个问题的基础上，如EC-Elgamal系统，以及数字签名、密钥交换协议等...
复制链接

扫一扫