C++

Forum

 
How to find whether a very long integer(upto 10^18) is a prime number or not in C++?

ayush5588 (1)
While solving a problem in C++,the range of numbers(integers only) to be used are from 4 to 10^18.So i am confused as how to calculate whether a number is a prime number or not upto such a long range?
Thank you in advance for any help..
Last edited on
Duthomhas (13208)
The Boost multiprecision library has a Miller-Rabin test built-in:

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
 
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/miller_rabin.hpp>
#include <chrono>
#include <iostream>
#include <random>

typedef boost::multiprecision::cpp_int number;

bool is_prime( const number& n )
{
  static std::ranlux48 rng( std::chrono::system_clock::now().time_since_epoch().count() );
  return miller_rabin_test( n, 40, rng );
}

int main()
{
  std::cout << "Test big numbers for primality.\n";
  
  number n;
  {
    std::cout << "n? ";
    std::cin >> n;
  }
  if (is_prime( n )) std::cout << "prime!\n";
  else               std::cout << "not prime.\n";
}

Running it:
$ a
Test big numbers for primality.
n? 123456789123456789000000000000000000000000000000000000000000023
prime!

This, of course, works for smaller numbers too.

Enjoy!
Last edited on
closed account (1vf9z8AR)
alternatively, you could take the number as a string and do the required divisions with strings but its a huge headache.
MikeStgt (466)
prime!
Probably not composite with an uncertainness of 8,3*10-25 or less (much less). How to push the MRT to a deterministic one is described here: https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants
in place of algorithm as pseudocode German Wikipedia shows a C++ function: https://de.wikipedia.org/wiki/Miller-Rabin-Test#Implementierung
alas for numbers up to 232 only, to go also up to 1018 some adjustments are necessary.

ayush5588, now decide if you do need 100% certainty or if a iota less is reasonable.

Edit: For the a. m. minor adjustments to make that MRT work also up to 1018 see here: http://www.cplusplus.com/forum/general/253330/#msg1115850
Type uint32_t is now uint64_t and what was uint64_t now __uint128_t
Last edited on
Duthomhas (13208)
LOL. Welcome to the “I can look stuff up on Wikipedia” section of our show.


Probably not composite with an uncertainness of 8,3*10-25 or less (much less)

You don’t consider 2-80 to be a sufficiently low probability of failure?

Because that’s significantly less probable than a passing neutrino changing the result of the function as it returns. That is, according to the experts. You know, the same experts who gave you that number in your wiki browsing. Yeah, I read them too.

If you need a better guarantee, use AKS or Baillie-PSW; both of which, BTW, are mentioned in that wiki link (and the second of which, while still “probabilistic”, turns non-failure probability up to 11).


See, “probable” is a “strictly-speaking” term used by mathematicians to say (in this case) that: no, we did not exhaustively compute against every possible divisor of n to prove it is prime.

But we can still be pretty darn smug about the result after 40 rounds of randomized tests.


[suggestions involving reading Wikipedia and vague FUD about “adjustments are necessary” when dealing with OP’s n compared to smaller n]

I admit, that’s some pretty impressive reasoning there... Introduce non-issues and scare the OP. Good job!

So... if you really wanna help a brother/sister out, then post some code, dude.

’Cause, you know, doing this stuff is really easy for beginners. Maybe only a couple lines of code, right?


MikeStgt has a past habit of visiting threads in an effort to derail my answers. It seems to irk him that there do exist a few small areas of computer science in which my knowledge exceeds his. He also enjoys piña coladas and long walks in the rain, and on cold nights will often find himself snuggling up with old BASIC programs so that he can untangle their spaghetti code. I have my ideas about why, but I could be wrong so I’ll keep them to myself.

lastchance (6980)
ayush5588 wrote:
While solving a problem in C++,the range of numbers(integers only) to be used are from 4 to 10^18


@ayush5588, that is well within the range of unsigned long long and standard direct prime-testing routines:
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
 
#include <iostream>
#include <limits>
#include <cmath>
using namespace std;

const unsigned long long MAXNUM = numeric_limits<unsigned long long>::max();


bool isPrime( unsigned long long n )
{
        if ( n     <  2 ) return false;
   else if ( n     <  4 ) return true;
   else if ( n % 2 == 0 ) return false;

   unsigned long long imax = sqrt( n + 0.5 );
   for ( unsigned long long i = 3; i <= imax; i += 2 )
   {
      if ( n % i == 0 ) return false;
   }
   return true;
}


int main()
{
   unsigned long long N;
   cout << MAXNUM << " is the largest allowed. Enter a number:\n";
   cin >> N;
   cout << N << ( isPrime( N ) ? " is " : " is not " ) << "prime\n";

   for ( N--; N > 1 && !isPrime( N ); N-- );
   if ( N > 1 ) cout << "The next lower prime is " << N;
}


18446744073709551615 is the largest allowed. Enter a number:
1234567890123456789
1234567890123456789 is not prime
The next lower prime is 1234567890123456737


18446744073709551615 is the largest allowed. Enter a number:
1000000000000000000
1000000000000000000 is not prime
The next lower prime is 999999999999999989
Last edited on
anup30 (968)
lol, i was also about to post as lastchance.
unsigned long long will support all 19 digit numbers.
so it's enough for op's problem.
MikeStgt (466)
after 40 rounds of randomized tests ...

... you assert "Prime!" and I conclude "very probable not a composite." This is the beginners forum and for that I prefere to be precise. As in education, trueness instead of half-truth.

BTW, it is not the first time you complain, in the past twice even very impertinently. Before you derail once more and forget your language, with C++ I am a novice too and try to learn by the problems of others -- problems with C++, not with "my knowledge exceeds his"-kindergarten competition.

Back to the subject. I am just about to compare MRT (code from a. m. German Wikipedia) against the very basic primality check from JLBorges http://www.cplusplus.com/forum/beginner/253520/#msg1114784 (modified: unsigned instead of int)
Result:
MRT counts 9591 primes from 1 to 100001 within 0.033 seconds.
SPC counts 9592 primes from 1 to 100001 within 0.006 seconds.
MRT counts 6247 primes from 8000000 to 8100000 within 0.051 seconds.
SPC counts 6247 primes from 8000000 to 8100000 within 0.042 seconds.
MRT counts 447 primes from 4294957295 to 4294967295 within 0.008 seconds.
SPC counts 446 primes from 4294957295 to 4294967295 within 8.764 seconds.

Summary: i) The claim in Wikipedia "Diese C++-Implementierung kann alle Zahlen kleiner 232 behandeln" is wrong, it fails for 1, _2,_ and 3.
ii) the primality check derivative (unsigned instead of int) from JLBorges fails for 4294967291.
iii) Wikipedia claims "if n < 9,080,191, it is enough to test a = 31 and 73"; with the given C++ code it fails.

Conclusion: Wikipedia needs some maintenance.

Edit: 2 is not a failure as MRT is for testing odd only -- https://de.wikipedia.org/wiki/Miller-Rabin-Test#Implementierung
Edit: see underscored
Last edited on
MikeStgt (466)
@lastchance
If you wait 4 seconds or several more makes a little difference. I hope you do not mind if I add a little wheel to your primality check. A small step for me ... suffices.
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
 
#include <iostream>
#include <cmath>
using namespace std;

bool isPrime( unsigned long long n )
{
   if ( n     <  2 ) return false;
   if ( n     <  4 ) return true;
   if ( n % 2 == 0 ) return false;
   if (!(n % 3)) return false;
   if (!(n % 5)) return false;
   if ( n < 7*7) return true;
   int wheel[] = {7, 11, 13, 17, 19, 23, 29, 31};
   unsigned long long imax = sqrt( n + 0.5 );    // + .5 as safety strategy
   for ( unsigned long long i = 0; i <= imax; i += 30 )
   
      for (int k = 0; k < 8; k++ )
      {
         if ( n % (i + wheel[k]) == 0 )
            return false;
      }
  }
  return true;
}


int main()
{
   unsigned long long N;
   cout << ~0ull << " is the largest allowed. Enter a number:\n";
   cin >> N;
   cout << N << ( isPrime( N ) ? " is " : " is not " ) << "prime\n";

   for ( N--; N > 1 && !isPrime( N ); N-- );
   if ( N > 1 ) cout << "The next lower prime is " << N;
}

Edit: see comment in bold
Last edited on
Duthomhas (13208)
MikeStgt also likes slow dances with strawman and listens to music by XY. A thing of evidence means nothing with but a liberal sprinkling of big words and tangenterines.

Fortunately, wholly missing the point always manages to keep him from considering himself speared in any concrete context, and he lives happily thinking airy and important thoughts designed to redirect the unsavvy beginner through deep waters that would otherwise have been blithely left uncrossed.

How very noble!
anup30 (968)
MikeStgt wrote:
"
MRT counts 9591 primes from 1 to 100001 within 0.033 seconds.
SPC counts 9592 primes from 1 to 100001 within 0.006 seconds.
MRT counts 6247 primes from 8000000 to 8100000 within 0.051 seconds.
SPC counts 6247 primes from 8000000 to 8100000 within 0.042 seconds.
MRT counts 447 primes from 4294957295 to 4294967295 within 0.008 seconds.
SPC counts 446 primes from 4294957295 to 4294967295 within 8.764 seconds.
"
Miller–Rabin Test is a Probabilistic test. it actually checks if a number is composite. but there is probablity that it can fail to find compositeness. so, MRT can give more number for range test, than SPC(simple primality check?). how "MRT counts 9591" which is less than "SPC counts 9592" ? can you share the MRT code you used.
MikeStgt (466)
Miller–Rabin Test is a Probabilistic test.

Correct. But Wikipedia also shows a chapter 'Deterministic variants', see
https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants

The C++ code I use is from here: https://de.wikipedia.org/wiki/Miller-Rabin-Test#Implementierung
Alas I had to add few filters, see following code. BTW, I changed it meanwhile, added a simple Primality check With Wheel (PWW) besides the Simple Primality Check (SPC), changed the numbers range, so the output is now:

MRT counts 78498 primes from 1 to 1000001 within 0.909 seconds.
SPC counts 78498 primes from 1 to 1000001 within 0.273 seconds.
PWW counts 78498 primes from 1 to 1000001 within 0.196 seconds.
MRT counts 45166 primes from 4000000000 to 4001000000 within 1.728 Seconds.
SPC counts 45166 primes from 4000000000 to 4001000000 within 13.172 seconds.
PWW counts 45166 primes from 4000000000 to 4001000000 within 9.588 seconds.
MRT counts 21 primes from 18446744073709550615 to 18446744073709551615 within 0.007 seconds.
SPC skipped.
PWW counts 21 primes from 18446744073709550615 to 18446744073709551615 within 310.68 seconds.

--------------------------------
Process exited after 337.7 seconds with return value 0


The code is "work in progress", not for use, much of it will not make it to a final version.
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
 
/*
 *  Teste Laufzeit verschiedener Tests auf prim.
 *  Dazu werden jeweis eine messbare Anzahl von Zahlen
 *  verschiedener Grösse getestet:
 *  a) 1..100'000,
 *  b) 8'000'000..8'100'000,
 *  c) 1'000 vor ~0 bis ~0.
 *
 *	0. Version: 190519, MF
 *
 */

#include <iostream>		/* cout, cin, endl */
#include <ctime>		/* clock */
#include <cmath>		/* sqrt */
#include <cstdint>		/* __uint128_t */
#include <bitset>		/* for bool array of prime Y/N */
using namespace std;

#define tytok unsigned long long
 
/* from https://de.wikipedia.org/wiki/Miller-Rabin-Test#Implementierung */
#include <cstdint>
// using std::uint32_t;
// bool mrt (const uint32_t n, const uint32_t a) { // n ungerade, 1 < a < n-1
bool mrt (const tytok n, const tytok a) {		 // n ungerade, 1 < a < n-1
   if( n < 2 ) return false;	// MF
   if( n == 2 ) return true;	// MF
   if( n == 3 ) return true;	// MF
   const tytok m = n - 1;
   tytok d = m >> 1, e = 1;
   while (!(d & 1)) d >>= 1, ++e;
   __uint128_t p = a, q = a;
//   unsigned long long p = a, q = a;
   while (d >>= 1) { // exponentiere modular: p = a^d mod n
      q *= q, q %= n; // quadriere modular: q = q^2 mod n
      if (d & 1) p *= q, p %= n; // multipliziere modular: p = (p * q) mod n
   }
   if (p == 1 || p == m) return true; // n ist wahrscheinlich prim
   while (--e) {
      p *= p, p %= n;
      if (p == m) return true;
      if (p <= 1) break;
   }
   return false; // n ist nicht prim
}

/* function to check if n is a prime number uses simple trial division */
bool is_prime( tytok n )
{
    if( n < 2 ) return false ;
    if( n == 2 ) return true ; // 2 is a prime number
    if( n%2 == 0 ) return false ; // even numbers greater than two are not prime

    // trial division by odd numbers up to the square root of n
    for( tytok div = 3 ; (div*div) <= n ; div += 2 ) if( n%div == 0 ) return false ;
    return true ; // not evenly divisible by any of the candidates
}

/* source: https://www.codeproject.com/articles/465041/primality-test */
bool IsPrimeWithWheel(tytok number)
{
	if ( number  < 2 ) return false;
	if ( number  < 4 ) return true;
	if ( number == 5 ) return true;
	if (!(number % 2)) return false;
	if (!(number % 3)) return false;
	if (!(number % 5)) return false;
	if ( number < 7*7) return true;
	int wheel[] = {7, 11, 13, 17, 19, 23, 29, 31};
	tytok rt = sqrt(number);
	for(tytok r = 0; r < rt; r += 30)
		for(int i = 0; i < 8; ++i)
                if (!(number % (r + wheel[i])))
                	return false;
	return true;
}

int main()
{
	tytok lo(1), max(~0), mid(4000000000), delta1(1000000), delta2(1000);
	tytok n, cnt(0);

	clock_t t0 = clock();		/* time(R) */
	for (n = lo; n <= lo+delta1; ++n)
		if (mrt(n, 2)
		 && mrt(n, 3))		cnt++;
	double te = (double)( clock() - t0 ) / CLOCKS_PER_SEC;  /* time(E) */
	cout << "MRT counts " << cnt << " primes from " << lo << " to " << lo + delta1 << " within " << te << " seconds.\n";
	
	cnt = 0; t0 = clock();		/* time(R) */
	for (n = lo; n <= lo+delta1; ++n)
		if (is_prime(n)) cnt++;
	te = (double)( clock() - t0 ) / CLOCKS_PER_SEC;  /* time(E) */
	cout << "SPC counts " << cnt << " primes from " << lo << " to " << lo + delta1 << " within " << te << " seconds.\n";

	cnt = 0; t0 = clock();		/* time(R) */
	for (n = lo; n <= lo+delta1; ++n)
		if (IsPrimeWithWheel(n)) cnt++;
	te = (double)( clock() - t0 ) / CLOCKS_PER_SEC;  /* time(E) */
	cout << "PWW counts " << cnt << " primes from " << lo << " to " << lo + delta1 << " within " << te << " seconds.\n";
	
	for (n = lo; n <= lo+delta1; ++n)
		if ((IsPrimeWithWheel(n)) != (mrt(n, 2) && mrt(n, 3)))
			cout << "\nDiff: " << n << ", mrt(n, 2) .., 3): " << mrt(n, 2) << mrt(n, 3) << "\n\n";

	cnt = 0; t0 = clock();		/* time(R) */
	for (n = mid; n <= mid+delta1; ++n)
//		if (mrt(n, 31) && mrt(n, 73))
		if (mrt(n, 2)
		 && mrt(n, 7)
		 && mrt(n, 61))		cnt++;
	te = (double)( clock() - t0 ) / CLOCKS_PER_SEC;  /* time(E) */
	cout << "MRT counts " << cnt << " primes from " << mid << " to " << mid + delta1 << " within " << te << " seconds.\n";
	
	cnt = 0; t0 = clock();		/* time(R) */
	for (n = mid; n <= mid+delta1; ++n)
		if (is_prime(n)) cnt++;
	te = (double)( clock() - t0 ) / CLOCKS_PER_SEC;  /* time(E) */
	cout << "SPC counts " << cnt << " primes from " << mid << " to " << mid + delta1 << " within " << te << " seconds.\n";

	cnt = 0; t0 = clock();		/* time(R) */
	for (n = mid; n <= mid+delta1; ++n)
		if (IsPrimeWithWheel(n)) cnt++;
	te = (double)( clock() - t0 ) / CLOCKS_PER_SEC;  /* time(E) */
	cout << "PWW counts " << cnt << " primes from " << mid << " to " << mid + delta1 << " within " << te << " seconds.\n";

	cnt = 0; t0 = clock();		/* time(R) */
	for (n = max - delta2; n < max; ++n)
		if (mrt(n,  2) && 
		    mrt(n,  3) &&
		    mrt(n,  5) &&
		    mrt(n,  7) &&
		    mrt(n, 11) &&
		    mrt(n, 13) &&
		    mrt(n, 17) &&
		    mrt(n, 19) &&
		    mrt(n, 23))		cnt++;
	te = (double)( clock() - t0 ) / CLOCKS_PER_SEC;  /* time(E) */
	cout << "MRT counts " << cnt << " primes from " << max - delta2 << " to " << max << " within " << te << " seconds.\n";
/*	
/*	cnt = 0; t0 = clock();		/* time(R) */
/*	for (n = max-delta2; n < max; ++n)
/*		if (is_prime(n)) cnt++;
/*	te = (double)( clock() - t0 ) / CLOCKS_PER_SEC;  /* time(E) */
/*	cout << "SPC counts " << cnt << " primes from " << max - delta2 << " to " << max << " within " << te << " seconds.\n";
*/
	cout << "SPC skipped.\n";

	cnt = 0; t0 = clock();		/* time(R) */
	for (n = max-delta2; n < max; ++n)
		if (IsPrimeWithWheel(n)) cnt++;
	te = (double)( clock() - t0 ) / CLOCKS_PER_SEC;  /* time(E) */
	cout << "PWW counts " << cnt << " primes from " << max - delta2 << " to " << max << " within " << te << " seconds.\n";
/*
	for (n = max-delta2; n < max; ++n)
		if ((is_prime(n)) != (mrt(n, 2) && mrt(n, 7) && mrt(n, 61)))
			cout << "\nDiff: " << n << " SPC mrt(n, 2) ..7 ..61: " << is_prime(n) << ' ' << mrt(n, 2) << mrt(n, 7) << mrt(n, 61) << "\n\n";
*/
}
helios (17574)
Correct. But Wikipedia also shows a chapter 'Deterministic variants', see
Right at the top of that same article it states that the deterministic algorithm assumes the extended Riemann hypothesis being true. So unless you've proven that, you aren't really increasing the certainty of the answer by using Miller's algorithm.
MikeStgt (466)
the deterministic algorithm assumes the extended Riemann hypothesis being true
Yes, I saw that hint too -- but my goal is to stay in the range up to 264 so to my ratiocination, performing the required tests as listed here
https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases
is sufficing to call it deterministic. Even worse than ignoring the nonproven extended Riemann hypothesis -- I hope the Wikipedia article could be correct. I confess, I did not read the original papers referenced by a. m. Wikipedia. Did you? ;)
anup30 (968)
last paragraph of wikipedia>Miller–Rabin primality test>Deterministic variants:
The Miller test is not used in practice. For most purposes, proper use of the probabilistic Miller–Rabin test or the Baillie–PSW primality test gives sufficient confidence while running much faster. It is also slower in practice than commonly used proof methods such as APR-CL and ECPP which give results that do not rely on unproven assumptions. For theoretical purposes requiring a deterministic polynomial time algorithm, it was superseded by the AKS primality test, which also does not rely on unproven assumptions.
MikeStgt (466)
Is there somewhere a C++ routine available for the AKS primality test?
Edit: no more questions ;)
Last edited on
anup30 (968)
>"Is there somewhere a C++ routine available for the AKS primality test?"
check out this: https://sourceforge.net/projects/aksprimality/
MikeStgt (466)
check out
7 MB of code + lib -- that's all? No, kidding ;)
Thanks for the link.
Topic archived. No new replies allowed.