fork download
copy 
  1. #include <iostream>
  2. #include <math.h>
  3. typedef long long int i64;
  4.  
  5. int jacobi(i64 k, i64 n)
  6. {
  7. 	int	res = 1;
  8. 	i64	kmn;
  9.  
  10. 	if( k < 0 ) { k = -k; if( n & 2 ) { res = -res; } }
  11.  
  12. 	while( (kmn = k % n) != 0 ) {
  13. 		k = n; n = kmn;
  14. 		while( !(n & 3) ) { n >>= 2; }
  15. 		if( !(n & 1) ) {
  16. 			n >>= 1;
  17. 			if( (k ^ (k >> 1)) & 2 ) { res = -res; }
  18. 		}
  19. 		if( (n & k & 2) ) { res = -res; }
  20. 	}
  21. 	return ( n == 1 ? res : 0 );
  22. }
  23.  
  24. i64 modprod(i64 a, i64 b, i64 mod)	// mod < 2^43
  25. {
  26. 	if( !(a & 0xffffffff80000000) && !(b & 0xffffffff80000000) ) {
  27. 		return (a * b) % mod;
  28. 	}
  29.  
  30. 	i64 a2 =  a / (1LL << 34);
  31. 	i64 b2 =  b / (1LL << 34);
  32. 	i64 a1 = (a / (1LL << 17)) % (1LL << 17);
  33. 	i64 b1 = (b / (1LL << 17)) % (1LL << 17);
  34. 	i64 a0 =  a % (1LL << 17);
  35. 	i64 b0 =  b % (1LL << 17);
  36.  
  37. 	i64 rem =			(a2 * b2) % mod;
  38. 	rem = ((rem << 17) + a2 * b1 + a1 * b2) % mod;
  39. 	rem = ((rem << 17) + a2 * b0 + a1 * b1 + a0 * b2) % mod;
  40. 	rem = ((rem << 17)			 + a1 * b0 + a0 * b1) % mod;
  41. 	rem = ((rem << 17)					   + a0 * b0) % mod;
  42.  
  43. 	return rem;
  44. }
  45.  
  46. i64 modpow(i64 base, i64 power, i64 mod)
  47. {
  48. 	i64 result = 1;
  49. 	while( power > 0 ) {
  50. 		if( power & 1 ) { result = modprod(result, base, mod); }
  51. 		base = modprod(base, base, mod);
  52. 		power >>= 1;
  53. 	}
  54. 	return result;
  55. }
  56.  
  57. bool is_prime(i64 n)
  58. {
  59. 	// Miller-Rabin (base 2)
  60. 	{	i64 d = (n - 1) >> 1;
  61. 		int s = 1;
  62. 		while( (d & 1) == 0 ) { d >>= 1; s ++; }
  63. 		i64 y = modpow(2, d, n);
  64. 		if( y != 1 ) {
  65. 			for(int t = s; y != n - 1 ;y = modprod(y, y, n)) {
  66. 				if( -- t == 0 ) { return false; }
  67. 			}
  68. 		}
  69. 	}
  70.  
  71. 	// square number check
  72. 	{	i64 s = (i64)(sqrt((double)n) + .5);
  73. 		if( s * s == n ) { return false; }
  74. 	}
  75.  
  76. 	// Lucas
  77. 	{	int d = 5, s = 1;
  78. 		while( jacobi(d * s, n) >= 0 ) { d += 2; s = -s; }
  79. 		d *= s;
  80. 		int p = 1, q = (1 - d) / 4;
  81. 		i64 k = n + 1, u = 1, v = p, m = 1, qm = q;
  82.  
  83. 		for(int b = 32; b > 0 ;b >>= 1) { if( k >= (m << b) ) { m <<= b; } }
  84.  
  85. 		while( (m >>= 1) > 0 ) {
  86. 			u = modprod(u, v, n);
  87. 			v = (modprod(v, v, n) - qm * 2) % n;
  88. 			qm = modprod(qm, qm, n);
  89.  
  90. 			if( (k & m) ) {
  91. 				i64 u2 = modprod(p, u, n * 2) + v;
  92. 				i64 v2 = modprod(d, u, n * 2) + modprod(p, v, n * 2);
  93. 				u = ( (u2 & 1) ? ((u2 + n) / 2) % n : (u2 / 2) % n );
  94. 				v = ( (v2 & 1) ? ((v2 + n) / 2) % n : (v2 / 2) % n );
  95. 				qm = modprod(qm, q, n);
  96. 			}
  97. 		}
  98. 		if( u != 0 ) { return false; }
  99. 	}
  100.  
  101. 	return true;
  102. }
  103.  
  104. int main()
  105. {
  106. 	int found = 0, s = 0, d = 4;
  107. 	for(i64 p = 5; p < (1LL << 43) ;p += (d = 6 - d)) {
  108. 		if( is_prime(p) ) {
  109. 			if( (s += d - 3) == 0 ) { std::cout << p << std::endl; }
  110. 		}
  111. 	}
  112. 	return 0;
  113. }
  114.  
Time limit exceeded #stdin #stdout 5s 4432KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
7
13
19
37
43
79
163
223
229
https://ideone.com/dJhsj8
language:
C++ (gcc 8.3)
created:
3 years ago
visibility:
public

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!