• 【学习笔记】Pollard's rho学习


    开坑时间:2019.10.18周五

    完善时间:2019.11.16周六

    学习原因及其他

    没什么原因,就是想学。有可能是因为今天在机房,csl到处问Pollard's rho怎么写,我随即发现自己不会,决定去学习。

    2019-10-18

    入门

    入门,初步学习:xyx的博客

    初步了解Pollard's rho的过程。认识到它的本质以及大致过程。

    • 要分解(N)

    • 伪随机生成数列(A={a_1,a_2,a_3,...}),每一项由前一项用固定方式推得。常用方式为(a_i=f(a_{i-1})=(a_{i-1}^2+C)mod\,N)其中(C)为一不变常数。

    • 每次查看(gcd(a_{2i}-a_i,N)),若不为1,说明此时找到了一个约数(gcd(a_{2i}-a_i,N))。若此时(gcd eq N),说明我们确实找到了一个非平凡约数。但如果(gcd=N),就没有找下去的必要了。

    • 这实际上是找到了循环节。即,如果最终我们有一个因数(P),那么(a_i\,mod\,P)实际是有循环节的。当(gcd eq 1)是,我们发现(a_{2i})(a_i)实际上就到了循环节上一个相同的位置。

    • 如果此时(gcd=N),就说明(a_{2i})(a_i)也到了(mod\,N)循环节上的同一位置。这就说明不行了,我们要换一个(C)

    • 我们假设要经过(Q)步进入循环节,循环节长(R)步,那么(a_{2i})(a_{i})要对应上循环节同一点,所需时间是(O(Q+R))的。

    • 所以这个算法的关键就是期望多少次,(a_i\,mod\,P)会冲突。生日悖论告诉我们是(O(sqrt{P}))(O(N^{frac{1}{4}}))

    生日悖论

    进一步学习:关于生日悖论发现一篇文章:一篇文章

    这篇文章关于生日悖论讲了一些。可以帮助理解为什么会在很快的时间内成功,但是没有证明期望。

    关于名字

    发现了为什么要叫Pollard's rho:因为循环节这东西长得很像( ho)

    关于(a_{2i}-a_i)

    这东西实际上就是Floyd的找环算法。

    2019-10-21

    代码1

    LL PR_zuo(LL cn,LL a)
    {
    	LL lin1 = 2, lin2 = get_ne(lin1,a,cn), lin3;
    	while(1)
    	{
    		lin3 = gcd((lin2-lin1+cn)%cn, cn);
    		if(lin3 != 1) return lin3 == cn ? 0 : lin3;
    		lin1 = get_ne(lin1,a,cn); lin2 = get_ne(get_ne(lin2,a,cn),a,cn);
    	}
    }
    LL Pollard_rho(LL cn)
    {
    	if(cn == 1) return 1;
    	if(!(cn&1)) return 2;
    	if(Miller_Rabin(cn)) return cn;
    	LL a = 4, lin;
    	while(1) {if(lin = PR_zuo(cn,a)) return lin; a += 3; }
    }
    

    这份代码时间复杂度(O(N^{frac{1}{4}}alpha(N))) ,其中(alpha(N))表示求(O(N))大小的两个数的gcd的时间复杂度。

    代码2

    LL PR_zuo(LL cn,LL a)
    {
    	LL lin1 = 2, lin2 = get_ne(lin1,a,cn), lin3;
    	int B = 100;
    	while(1)
    	{
    		LL lin = 1, lin4 = lin1, lin5 = lin2;
    		for(int i = 1;i<=B;i++) 
    		{
    			lin = ksc(lin,(lin5-lin4+cn)%cn,cn);
    			lin4 = get_ne(lin4,a,cn); lin5 = get_ne(get_ne(lin5,a,cn),a,cn);
    		}
    		lin3 = gcd(lin,cn);
    		if(lin3 != 1) {if(lin3 != cn) return lin3; break; }
    		lin1 = lin4; lin2 = lin5;
    	}
    	while(1)
    	{
    		lin3 = gcd((lin2-lin1+cn)%cn, cn);
    		if(lin3 != 1) return lin3 == cn ? 0 : lin3;
    		lin1 = get_ne(lin1,a,cn); lin2 = get_ne(get_ne(lin2,a,cn),a,cn);
    	}
    }
    LL Pollard_rho(LL cn)
    {
    	if(cn == 1) return 1;
    	if(!(cn&1)) return 2;
    	if(Miller_Rabin(cn)) return cn;
    	LL a = 4, lin;
    	while(1) {if(lin = PR_zuo(cn,a)) return lin; a += 3; }
    }
    

    这份代码是(O(N^{frac{1}{4}}+frac{N^{frac{1}{4}}}{B}alpha(N)+Balpha(N)))(B)是一个选定的常数。

    但是我这一份慢得很,不知道为什么。正在努力卡常。

    2019-10-27

    一个最新进展:我发现(a)初值选得好,会让程序效率有一定的提升(大概20%)

    13是一个不错的选择。

    2019-11-16

    今天CSP复赛Day1,下午不太想写题,所以就开始卡常。

    1. 快速乘不要用倍增,直接用__int128(没有c++11,手写也比倍增快。)
    2. 取模要科学,不要动不动取模。那份优化版(就是加了(B))的,就是因为这两件事才特别慢的。
    3. 快速乘就是(哔——)(哔——),一开始取模写错了,就很难受。直接让我死循环。

    代码——终极版

    #include<bits/stdc++.h>
    #define LL long long
    using namespace std;
    template<typename T>void Read(T &cn)
    {
    	char c;int sig = 1;
    	while(!isdigit(c = getchar()))if(c == '-')sig = -1; cn = c-48;
    	while(isdigit(c = getchar()))cn = cn*10+c-48; cn*=sig;
    }
    template<typename T>void Write(T cn)
    {
    	if(cn < 0) {putchar('-'); cn = 0-cn; }
    	int wei = 0; T cm = 0; int cx = cn%10; cn/=10;
    	while(cn)cm = cm*10+cn%10,cn/=10,wei++;
    	while(wei--)putchar(cm%10+48),cm/=10;
    	putchar(cx+48);
    }
    const int MAXPRI = 6;
    const int pri[] = {2,61,3,13,17,23};
    //inline LL ksc(LL cn,LL cm,LL MOD) {cn%=MOD; cm%=MOD; LL ans = 0; while(cm) ans = ans + ((cm&1) ? cn : 0), cn = cn+cn, cm>>=1, ans = ((ans>=MOD)?(ans-MOD):ans), cn = ((cn>=MOD)?(cn-MOD):cn); return ans; }
    inline LL ksc(__int128 cn,__int128 cm,LL MOD) {return cn*cm%MOD; }
    LL ksm(LL cn,LL cm,LL MOD) {LL ans = 1; while(cm) ans = ksc(ans,1+(cn-1)*(cm&1),MOD), cn = ksc(cn,cn,MOD), cm>>=1; return ans; }
    LL gcd(LL cn,LL cm) {return cm ? gcd(cm,cn%cm) : cn; }
    int MR_jian(LL cn, int cm, LL a, LL N)
    {
    	a = ksm(a, cn, N);
    	if(a == 1 || a == N-1) return 1;
    	while(cm--)
    	{
    		a = ksc(a,a,N);
    		if(a == 0 || a == 1) return 0;
    		if(a == N-1) return 1;
    	}
    	return 0;
    }
    int Miller_Rabin(LL cn)
    {
    	if(cn == 2) return 1;
    	if(cn == 0 || cn == 1 || (!(cn&1))) return 0;
    	LL wei = 0, cm = cn-1, cx;
    	while(!(cm&1)) wei++, cm>>=1;
    	for(int i = 0;i<MAXPRI;i++) if(pri[i] == cn) return 1;
    	for(int i = 0;i<MAXPRI;i++) if(!MR_jian(cm,wei,pri[i],cn)) return 0;
    	return 1;
    }
    inline LL get_ne(LL cn,LL N,LL a) {return (ksc(cn,cn,N) + a)%N; }
    LL Pollard_rho(LL cn)
    {
    	if(!(cn&1)) return 2;
    	LL a = 13;
    	int b = 100;
    	while(1)
    	{
    		LL lin1 = 2, lin2 = get_ne(lin1,cn,a);
    		while(1)
    		{
    			LL bx = lin1, by = lin2;
    			LL lei = 1;
    			for(int i = 1;i<=b;i++) lei = ksc(lei,bx-by+cn,cn), bx = get_ne(bx,cn,a), by = get_ne(get_ne(by,cn,a),cn,a);
    			LL lin = gcd(lei,cn);
    			if(lin != 1){
    				for(int i = 1;i<=b;i++)
    				{
    					lin = gcd(lin1-lin2+cn,cn);
    					if(lin == cn) break;
    					if(lin != 1) return lin;
    					lin1 = get_ne(lin1,cn,a); lin2 = get_ne(get_ne(lin2,cn,a),cn,a);
    				}
    				break;
    			}
    			lin1 = bx; lin2 = by;
    		}
    		a += 23;
    	}
    }
    LL n,t;
    int shi1,shi2,kai;
    LL get_zuida(LL cn)
    {
    	if(Miller_Rabin(cn)) return cn;
    	LL cm = Pollard_rho(cn);
    	return max(get_zuida(cm), get_zuida(cn/cm));
    }
    void zuo()
    {
    	Read(n); 
    	if(Miller_Rabin(n)) {puts("Prime"); return; }
    	Write(get_zuida(n)); puts("");
    }
    int main() {Read(t); while(t--) zuo(); return 0; }
    
    
  • 相关阅读:
    51nod 1122 机器人走方格 V4(矩阵乘法)
    51nod 1092 回文字符串
    51nod 1254 最大子段和 V2(递推)
    容斥原理求gcd为k的数对个数
    51nod 1115 最大M子段和 V3
    51nod 1053 最大M子段和 V2(贪心)
    洛谷P1792 [国家集训队]种树
    洛谷P1484 种树(反悔贪心,双向链表+堆)
    51 nod 1052 最大M子段和
    51 nod 1051 最大子矩阵和
  • 原文地址:https://www.cnblogs.com/czyarl/p/11701276.html
Copyright © 2020-2023  润新知