# Greatest prime number smaller than N where N can be as large as 10^18

This is code for finding largest prime number smaller than $N$, where $N$ can be as large as $10^{18}$. But, it takes 1 minute for $10^{18}$. I need to pass it in 2-3 sec.

What changes should I make to pass it? My compiler used is g++ 4.3.2

The program here runs for multiple test cases.

#include <iostream>
#include <cmath>

#define c2 341550071728321
#define c1 4759123141

using namespace std;

unsigned long long int mul(unsigned long long int a, unsigned long long int b, unsigned      long long int mod)
{
int i;
unsigned long long int now = 0;
for (i = 63; i >= 0; i--) if (((a >> i) & 1) == 1) break;
for (; i >= 0; i--)
{
now <<= 1;
while (now > mod) now -= mod;
if (((a >> i) & 1) == 1) now += b;
while (now > mod) now -= mod;
}
return now;
}

unsigned long long int pow(unsigned long long int a, unsigned long long int p, unsigned long long int mod)
{
if (p == 0) return 1;
if (p % 2 == 0) return pow(mul(a, a, mod), p / 2, mod);
return mul(pow(a, p - 1, mod), a, mod);
}

bool MillerRabin(unsigned long long int n)
{
int l;
unsigned long long int ar = { 2, 3, 5, 7, 11, 13, 17, 19, 23};
if (n < c1) {

l = 3;
}
else if (n < c2) {
l = 7;
}
else {

l = 9;
}
//l = 9;
unsigned long long d = n - 1;
int s = 0;
while ((d & 1) == 0) { d >>= 1; s++; }
int i, j;
for (i = 0; i < l; i++)
{
unsigned long long int a  = min(n - 2, ar[i]);
unsigned long long int now = pow(a, d, n);
if (now == 1) continue;
if (now == n - 1) continue;
for (j = 1; j < s; j++)
{
now = mul(now, now, n);
if (now == n - 1) break;
}
if (j == s) return false;
}
return true;
}

bool check_prime(unsigned long long int n) {

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

int main()
{
int t;

cin >> t;

while(t--) {
unsigned long long int n;
cin >> n;
while(1) {
if(check_prime(n)) {
cout << n << endl;
break;
}
else {
n--;
}
}
}

}

• Basically you go in the right direction: Miller-Rabin test, fast modular exponentiation. The calculation should be very fast (a lot less than 2-3 sec). So check where you loose time (profiling, or simple measure the time with with appropriate statements in the code) and find the error (I am not sure if your pow function is alright but I want not think about it). Oct 14 '14 at 16:35

My first suggestion is that you upgrade your compiler to gcc 4.8 or later, or clang 3.3 or later, so that you can compile with -std=c++11; C++11 really improves the language. I'll try to keep things to C++03, but I'll point out some opportunities for improvement with C++11.

As Jamal mentioned, unsigned long long is a mouthful. Using a typedef can allow you to change it more easily later, and allows you to document (at the point of the typedef) the requirements on the type.

// Integer must be able to represent at least the natural numbers up to 10^18.
typedef unsigned long long Integer;


I also recommend against using namespace std;. The top two answers to this SO question explain why this line is a bad practice and can even cause hard-to-diagnose bugs in your code.

Now let's see if we can clarify your Miller-Rabin code and identify areas for improvement. Let's start from the bottom and work up. You need to factor n - 1 into 2^s * d:

// n - 1 == d * 2^s
struct Factorization {
Integer n;
Integer s;
Integer d;
}

Factorization factor(Integer n) {
Factorization f;
f.n = n;
f.d = n - 1;
f.s = 0;
while (!(f.d % 2)) {
f.d /= 2;
++f.s;
}
return f;
}


Notice two things. I wrote the code in the way that is most straightforward for humans to read, not by trying to twiddle bits (e.g. d & 1, d >>= 1). It rarely pays to try to outsmart the compiler: write the operation you intend. Some time try writing both versions and benchmarking them with an optimization level of 2 or 3; I strongly doubt that the shifty version will be faster.

Second, I pulled out types and functions that make it clearer to readers what's going on. I had to look up a reference to Miller-Rabin to understand what in the world you were doing in your function, because it's all bit-shift this and if n < c1 that. (PS comments don't slow down your code :)

OK, now let's pull out a function to test a single base. A function called MillerRabin that returns bool is not very self documenting; let's see if we can fix that:

enum Primality {
COMPOSITE,
PROBABLY_PRIME,
};

Primality CheckMillerRabinWitness(const Factorization& f, Integer witness) {
Integer x = pow_mod(witness, f.d, f.n);
if (x == 1 || x == f.n - 1) return PROBABLY_PRIME;

for (int i = 0; i < s - 1; ++i) {
x = square_mod(x, f.n);
if (x == 1) return COMPOSITE;
if (x == f.n - 1) return PROBABLY_PRIME;
}
return COMPOSITE;
}


This function really lacks comments; here's where all the math is, and it could use a link to a reference, or maybe a couple of lemmas from number theory to explain to the reader what's going on (I'll leave this as an exercise :)

OK, we see we need two operations: square_mod, which squares a number modulo n, and pow_mod, which raises a number to a power modulo n. These can be implemented in terms of your pow and mul. We'll come back to these.

Now let's implement the whole algorithm:

Primality MillerRabinTest(Integer n) {
if (n <= 3) return COMPOSITE;
if (!(n % 2)) return COMPOSITE;
const Factorization f = factor(n);
const std::vector<Integer>& witnesses = SelectWitnesses(n);

for (unsigned i = 0; i < witnesses.size(); ++i) {
const Integer witness = witnesses[i];  // See note (1) below
if (CheckMillerRabinWitness(f, witness) == COMPOSITE) return COMPOSITE;
}
return PROBABLY_PRIME;
}


(1) with C++11, this line and the line previous can be replaced by

for (const Integer witness : witnesses) {


One thing to notice: The high-level algorithm here is absolutely clear: we factorize, we get a list of witnesses, and we see if any of the witnesses twig that n is composite. If there are no such witnesses, we admit that this number is probably prime.

Another thing we can notice: when n is even, we skip all the work. Your code performs pow(a, d, n) for each witness a whenever n is even. I wouldn't be surprised to see a 10-33% performance improvement from that change alone.

OK, now let's implement SelectWitnesses:

const std::vector<Integer>& SelectWitnesses(Integer n) {
// Note: http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
// lists minimal sets of witnesses that are sufficient for deterministic
// testing of primes for certain ranges of inputs. We've selected only
// a few of those ranges; it may be worthwhile to add more ranges here.
static const std::vector<Integer> w2k = { 2 };
static const std::vector<Integer> w25M = { 2, 3, 5 };
static const std::vector<Integer> w3T = { 2, 3, 5, 7, 11, 13 };
static const std::vector<Integer> w341T = { 2, 3, 5, 7, 11, 13, 17 };
static const std::vector<Integer> wmax = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };

if (n < 2047ull) { return w2k; }
else if (n < 25326001ull) { return w25M; }
else if (n < 3474749660383ull) { return w3T; }
else if (n < 341550071728321ull) { return w341T; }
else if (n < 3825123056546413051ull) { return wmax; }
else { throw std::out_of_range("Input is too large."); }
}


Here, we have a potential efficiency gain over your old function. In your implementation, the integer array had to be written onto the stack for each call to MillerRabin. Here, each vector is created exactly once per invocation of your program (NOTE: as with all performance-related issues, PROFILE BOTH VERSIONS. You may find that the memory locality provided by your implementation is more important than the cost of writing 9 values to the stack on every call to MillerRabin.). You'll also be somewhat faster for smaller numbers, but that is probably a small effect.

Note: I've written the C++11 version of this function (using braced initialization for the vectors) because the C++03 version is too annoying to write. It's doable though.

OK, check_prime:

bool is_prime(Integer n) {
return MillerRabinTest(n) == PROBABLY_PRIME;
}


First of all, I've renamed the function -- a name like check_prime doesn't tell you what true means -- is_prime clearly returns true if its argument is prime and false otherwise. Also, I've deleted all the dead code in there -- it's bad practice to keep code around that can't be executed.

OK, this review is already getting on pretty far, so I don't have time to go into pow and mul, but I suspect that there's some room for improvement there. Again, don't try to fool the compiler -- if you're trying to figure out whether something is even, check against x % 2 instead of x & 1, and if you're trying to divide by 2, divide by 2 instead of shifting right.

You'll definitely want to use a profiler. I'm partial to Google's gperftools, but there are other options. Focus your efforts on the functions that you spend the most time on. You can also use your profiler to do line-by-line annotations on your source to see where the hot spots are.

You'll also want to write some unit tests sooner than later. As you're monkeying around with pow and mul (and other functions, of course), it'll be very easy to introduce a bug, and unit tests will help prevent that.

Disclaimer: I haven't tried to compile, much less tested, any of the code in this review.

• These macros are not needed in C++:

#define c2 341550071728321
#define c1 4759123141


They should be constants, using the const keyword (or constexpr in C++11):

const int c2 = 341550071728321;
const int c1 = 4759123141;


I've only kept the names the same here because I have no idea what they mean. Try to avoid using single-char variable names, unless their meaning is already obvious.

• You use unsigned long long in many places, which could look verbose. You could use a typedef so that you can use another (shorter) name in place of them:

typedef unsigned long long ull;


Better yet, you can use std::uint64_t from <cstdint>. It is shorter, and it should be guaranteed to hold 64 bits on any system.

• This is quite a bit of whitespace:

if (n < c1) {

l = 3;
}
else if (n < c2) {
l = 7;
}
else {

l = 9;
}


You already use four spaces for indentation, so it should be consistent:

if (n < c1) {
l = 3;
}
else if (n < c2) {
l = 7;
}
else {
l = 9;
}


You're also lacking indentation within main() and check_prime(), but it's okay in your other functions. I don't see the reason for all this inconsistency if you're not in a hurry to write the code itself. Make it clean now so that fine-tuning will be easier, both for you and for others.

There's a very simple improvement: You check all the integers using the expensive Miller-Rabin test. But if you take for example n = 999,999,999,999,999,998 then I can tell you that it's not a prime!

Before you start an expensive test, check whether n is divisible by 2, 3, 5, 7, 11 and so on. How far do you go? To decide whether to test for divisibility by 17, the chance is one in 17 that this test detects a composite number, 16/17 that it doesn't. If Miller-Rabin takes 17 times longer than division by 17, or worse, divide by 17. That alone should get you close to your time limit.

Now the function "mul" is very, very, slow. You want to calculate (a * b) modulo n, where a, b and n are up to 10^18. You use a loop that runs 63 times. To calculate (a * b) modulo n, you calculate div = (a * b) / n) and mod = (a * b) - n * div. Here's a way to calculate this quickly, assuming an Intel processor where "long double" has a 64 bit mantissa.

Now you know (hopefully) that unsigned long long multiplication gives the correct result, modulo 2^64. So you calculate

long double one_over_n = 1 / (long double) n;
unsigned long long div = (unsigned long long) (one_over_n * a * b);
unsigned long long mod = (a * b) - n * div;


What happens? a * b / n is calculated in long double floating point. Since the result is < 10^18 and there are 64 bits of mantissa, this is quite close to the correct result, even though you will have rounding errors. Converting to (unsigned long long) gives the correct result, or possibly one too large or too small. So mod is either correct, or too small by n, or too large by n. No problem. Add these two lines:

if (mod > (1ull << 63)) mod += n;
if (mod >= n) mod -= n;


This will run a lot faster if you change the recursion in pow to a loop, and calculate one_over_n once outside the loop. These two changes should get you under the time limit easily.