fork download 
  1. #include <bits/stdc++.h>
  2. #define PRIME 1
  3. #define COMPOSITE 0
  4. #define MAXNUMFACT 64
  5.  
  6. using namespace std;
  7. typedef unsigned long long LLU;
  8.  
  9. LLU Fa[MAXNUMFACT];
  10. int nfact = -1;
  11. map <LLU, LLU> mp;
  12.  
  13. int miller_rabin_testing(int acc, LLU N);
  14. void decompose(LLU p, int *k, LLU *m);
  15. LLU llu_fsqrt(LLU N);
  16. LLU mod_mul(LLU a, LLU b, LLU p);
  17. LLU random64();
  18. LLU next_prng(LLU x, LLU N);
  19. LLU gcd(LLU a, LLU b);
  20. LLU find_factor(LLU N);
  21. LLU mod_power(LLU a, LLU b, LLU p);
  22. void fast_native_fact(LLU N);
  23. int llu_cmp(const void *a, const void *b);
  24. void pollard_rho_fact(LLU N);
  25. void print_factors();
  26.  
  27. int main () {
  28.     ios::sync_with_stdio(false); cin.tie(0);
  29.     LLU n;
  30.     cin >> n;
  31. 	pollard_rho_fact(n);
  32. 	print_factors();
  33. 	return 0;
  34. }
  35.  
  36. void pollard_rho_fact(LLU(N)){
  37. 	if(N == LLU(1)) return;
  38. 	else if(N < (LLU(1) << 10)) fast_native_fact(N);
  39. 	else if(miller_rabin_testing(10, N) == PRIME) Fa[++nfact] = N;
  40. 	else {
  41. 		LLU p = (LLU)0;
  42. 		while(p == 0 || p == N) p = find_factor(N);
  43. 		pollard_rho_fact(p);
  44. 		pollard_rho_fact(N/p);
  45. 	}
  46. }
  47.  
  48. LLU find_factor(LLU N) {
  49. 	LLU x = random64()% N;
  50. 	LLU y = x;
  51. 	LLU p = 1;
  52. 	do{
  53. 		x = next_prng(x, N);
  54. 		y = next_prng(next_prng(y, N), N);
  55. 		if (x > y) p = gcd(x-y,N);
  56. 		else p = gcd(y-x,N);
  57. 	} while(p == 1);
  58. 	return p;
  59. }
  60.  
  61. LLU next_prng(LLU x, LLU N){
  62. 	LLU xs = mod_mul(x,x,N);
  63. 	return xs + LLU(1);
  64. }
  65.  
  66. void fast_native_fact(LLU N){
  67. 	while (N % 2 == 0){
  68. 		Fa[++nfact] = LLU(2);
  69. 		N = N/2;
  70. 	}
  71. 	while (N % 3 == 0){
  72. 		Fa[++nfact] = LLU(3);
  73. 		N = N/3;
  74. 	}
  75. 	if (N == 1) return;
  76. 	if (miller_rabin_testing(10, N) == PRIME ) {
  77. 			Fa[++nfact] = N;
  78. 			return;
  79. 	}
  80. 	LLU Q = 5;
  81. 	int diff = 2;
  82. 	int sqrtN = llu_fsqrt(N);
  83. 	while (Q <= sqrtN) {
  84. 		if (N %Q == 0) {
  85. 			Fa[++nfact] = Q;
  86. 			N = N/Q;
  87. 			sqrtN = llu_fsqrt(N);
  88. 			if (miller_rabin_testing(10, N) == PRIME ) {
  89. 				Fa[++nfact] = N;
  90. 				return;
  91. 			}
  92. 		}
  93. 		else {
  94. 			Q = Q + LLU(diff);
  95. 			diff = (diff&2) + 2;
  96. 		}
  97. 	}
  98. }
  99.  
  100. void print_factors() {
  101. 	for (int i = 0; i <= nfact; i++){
  102. 		mp[Fa[i]]++;
  103. 	}
  104. 	for (auto it : mp) {
  105.         cout << it.first << " " << it.second << "\n";
  106. 	}
  107. }
  108.  
  109. int llu_cmp(const void *a, const void *b){
  110. 	return *(LLU*)a - *(LLU*)b;
  111. }
  112.  
  113. int witness(LLU a,  LLU N){
  114.     LLU m;
  115. 	int k;
  116. 	decompose(N-1, &k, &m);
  117. 	LLU B[k+1];
  118. 	B[0] = mod_power(a, m, N);
  119. 	if(B[0] == 1) return PRIME;
  120. 	else {
  121. 		int i = 1;
  122. 		while(i <= k) {
  123. 			B[i] =  mod_mul(B[i-1], B[i-1], N);
  124. 			if(B[i] == 1) {
  125. 				if(B[i-1] == N-1) return PRIME;
  126.                 else return COMPOSITE;
  127. 			}
  128. 			i++;
  129. 		}
  130. 	}
  131. 	return COMPOSITE;
  132. }
  133.  
  134. int miller_rabin_testing(int acc, LLU N) {
  135. 	LLU a = 0;
  136. 	for (int i = 0; i <= acc; i++) {
  137. 		a = random64()%(N-2) + 2;
  138. 		if (gcd(a,N) != 1) return COMPOSITE;
  139. 		else if (mod_power(a, N-1, N) != 1) return COMPOSITE;
  140. 		else if (witness(a, N) == COMPOSITE) return COMPOSITE;
  141. 	}
  142. 	return PRIME;
  143. }
  144.  
  145. void decompose(LLU p, int *k, LLU *m){
  146. 	int i = 0;
  147. 	while ((p & 1) == 0){
  148. 		i++;
  149. 		p = p >> 1;
  150. 	}
  151. 	*k = i;
  152. 	*m = p;
  153. }
  154.  
  155. LLU mod_mul(LLU a, LLU b, LLU p){
  156. 	if (b == 1) return a % p;
  157. 	else {
  158. 		LLU x = mod_mul(a, b>>1, p);
  159. 		if ((b & 1) == 0) return (x<<1)%p;
  160. 		else return (((x<<1)%p)+a)%p;
  161. 	}
  162. 	return 0;
  163. }
  164.  
  165. LLU mod_power(LLU a, LLU b, LLU p){
  166. 	if (b == 1) return a % p;
  167. 	else {
  168. 		LLU x = mod_power(a, b/2, p);
  169. 		if((b & 1) == 0) return mod_mul(x, x, p);
  170. 		else return (mod_mul(mod_mul(x, x, p), a, p));
  171. 	}
  172. }
  173.  
  174. LLU gcd(LLU a, LLU b){
  175. 	LLU tmp;
  176. 	while  (b != 0){
  177. 		tmp = b;
  178. 		b = a % b;
  179. 		a = tmp;
  180. 	}
  181. 	return a;
  182. }
  183.  
  184. LLU random64() {
  185. 	LLU n  = 0;
  186. 	int l = rand()%(RAND_MAX-2)+2;
  187. 	int r = rand()%(RAND_MAX) + rand()%2;
  188. 	n = n | LLU(r);
  189. 	n = n << 32;
  190. 	n = n | LLU(l);
  191. 	return n;
  192. }
  193.  
  194. LLU llu_fsqrt(LLU  N){
  195. 	LLU X, A2, R;
  196. 	X = N;
  197. 	A2 = LLU(1) << 62;
  198. 	R = 0;
  199. 	while (N < A2) A2 >>= 2;
  200. 	while (A2 > 0) {
  201. 		if (X >= R + A2) {
  202. 			X -= R + A2;
  203. 			R = (R >> 1) + A2;
  204. 			A2 = A2>> 2;
  205. 		}
  206. 		else {
  207. 			A2 = A2 >> 2;
  208. 			R = R >> 1;
  209. 		}
  210. 	}
  211. 	return R;
  212. }
Success #stdin #stdout 0.01s 5532KB
comments (?)
stdin 
12
stdout 
2 2
3 1
https://ideone.com/xfeeYF
language:
C++14 (gcc 8.3)
created:
1 year ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!