https://www.geeksforgeeks.org/pollards-rho-algorithm-prime-factorization/

Pollard’s Rho Algorithm for Prime Factorization

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
/* C++ program to find a prime factor of composite using 
Pollard's Rho algorithm */
#include<bits/stdc++.h> 
#include <sys/time.h>
using namespace std; 
typedef long long int lli;

//微妙级随机数种子（ t=1 时），t = 1000 为毫秒级
void micro_srand(int t){
	struct timeval tv;
    gettimeofday(&tv,NULL);
    srand((tv.tv_sec * t ) + (tv.tv_usec / t));
}

/* Function to calculate (base^exponent)%modulus */
lli power(lli base, int exponent,lli modulus) 
{ 
	/* initialize result */
	lli result = 1; 

	while (exponent > 0) 
	{ 
		/* if y is odd, multiply base with result */
		if (exponent & 1) 
			result = (result * base) % modulus; 

		/* exponent = exponent/2 */
		exponent = exponent >> 1; 

		/* base = base * base */
		base = (base * base) % modulus; 
	} 
	return result; 
} 


// This function is called for all k trials. It returns 
// false if n is composite and returns false if n is 
// probably prime. 
// d is an odd number such that d*2<sup>r</sup> = n-1 
// for some r >= 1 
bool miillerTest(lli d, lli n) 
{ 
	// Pick a random number in [2..n-2] 
	// Corner cases make sure that n > 4 
	lli a = 2 + rand() % (n - 4); 

	// Compute a^d % n 
	lli x = power(a, d, n); 

	if (x == 1 || x == n-1) 
	return true; 

	// Keep squaring x while one of the following doesn't 
	// happen 
	// (i) d does not reach n-1 
	// (ii) (x^2) % n is not 1 
	// (iii) (x^2) % n is not n-1 
	while (d != n-1) 
	{ 
		x = (x * x) % n; 
		d *= 2; 

		if (x == 1)	 return false; 
		if (x == n-1) return true; 
	} 

	// Return composite 
	return false; 
} 

// It returns false if n is composite and returns true if n 
// is probably prime. k is an input parameter that determines 
// accuracy level. Higher value of k indicates more accuracy. 
bool isPrime(lli n, int k) 
{ 
	// Corner cases 
	if (n <= 1 || n == 4) return false; 
	if (n <= 3) return true; 

	// Find r such that n = 2^d * r + 1 for some r >= 1 
	lli d = n - 1; 
	while (d % 2 == 0) 
		d /= 2; 

	// Iterate given nber of 'k' times 
	for (int i = 0; i < k; i++) 
		if (!miillerTest(d, n)) 
			return false; 

	return true; 
} 

/* method to return prime divisor for n */
lli  PollardRho(lli n) 
{ 
	/* initialize random seed */
	micro_srand(1);

	/* no prime divisor for 1 */
	if (n==1) return n; 

	/* even number means one of the divisors is 2 */
	if (n % 2 == 0) return 2; 
	
	//如果是素数返回 
	if(isPrime(n,1)) return n;
	

	/* we will pick from the range [2, N) */
	lli x = (rand()%(n-2))+2; 
	lli y = x; 

	/* the constant in f(x). 
	* Algorithm can be re-run with a different c 
	* if it throws failure for a composite. */
	lli c = (rand()%(n-1))+1; 

	/* Initialize candidate divisor (or result) */
	lli d = 1; 

	/* until the prime factor isn't obtained. 
	If n is prime, return n */
	while (d==1) 
	{ 
		/* Tortoise Move: x(i+1) = f(x(i)) */
		x = (power(x, 2, n) + c + n)%n; 

		/* Hare Move: y(i+1) = f(f(y(i))) */
		y = (power(y, 2, n) + c + n)%n; 
		y = (power(y, 2, n) + c + n)%n; 

		/* check gcd of |x-y| and n */
		d = __gcd(abs(x-y), n); 

		/* retry if the algorithm fails to find prime factor 
		* with chosen x and c */
		if (d==n) return PollardRho(n); 
	} 

	return d; 
} 

/* driver function */
int main() 
{ 
	freopen("in","r",stdin);
	lli n,ans[100],t,cnt=0; 
	scanf("%lld",&n);
	while(n!=1){
		t=PollardRho(n);
		if(isPrime(t,1)) {
			ans[cnt++]=t;
			n/=t;
		}
	}
	sort(ans,ans+cnt);
	for(int i=0;i<cnt;i++)
		printf("%d ",ans[i]);
	return 0; 
}