Milan-Rabin FIPS 186-3 Primary Test

Question

Im trying to implement the Miller-Rabin primitiveness test as described in FIPS 186-3 C.3.1. No matter what I do, I cannot get it to work. The instructions are pretty specific, and I don't think I missed something, and yet Im getting true for non-simple values.

What have I done wrong?

 template <typename R, typename S, typename T> T POW(R base, S exponent, const T mod){ T result = 1; while (exponent){ if (exponent & 1) result = (result * base) % mod; exponent >>= 1; base = (base * base) % mod; } return result; } // used uint64_t to prevent overflow, but only testing with small numbers for now bool MillerRabin_FIPS186(uint64_t w, unsigned int iterations = 50){ srand(time(0)); unsigned int a = 0; uint64_t W = w - 1; // dont want to keep calculating w - 1 uint64_t m = W; while (!(m & 1)){ m >>= 1; a++; } // skipped getting wlen // when i had this function using my custom arbitrary precision integer class, // and could get len(w), getting it and using it in an actual RBG // made no difference for(unsigned int i = 0; i < iterations; i++){ uint64_t b = (rand() % (W - 3)) + 2; // 2 <= b <= w - 2 uint64_t z = POW(b, m, w); if ((z == 1) || (z == W)) continue; else for(unsigned int j = 1; j < a; j++){ z = POW(z, 2, w); if (z == W) continue; if (z == 1) return 0;// Composite } } return 1;// Probably Prime } 

 std::cout << MillerRabin_FIPS186(33) << std::endl; std::cout << MillerRabin_FIPS186(35) << std::endl; std::cout << MillerRabin_FIPS186(37) << std::endl; std::cout << MillerRabin_FIPS186(39) << std::endl; std::cout << MillerRabin_FIPS186(45) << std::endl; std::cout << MillerRabin_FIPS186(49) << std::endl; 

gives me:

 0 1 1 1 0 1 

Answer 1

The only difference between your implementation and Wikipedia is that you forgot the second returned compound statement. You should get the opposite 0 at the end of the loop.

Edit: As Daniel noted, there is a second difference. The continuation continues the inner loop, not the outer loop, as intended.

 for(unsigned int i = 0; i < iterations; i++){ uint64_t b = (rand() % (W - 3)) + 2; // 2 <= b <= w - 2 uint64_t z = POW(b, m, w); if ((z == 1) || (z == W)) continue; else{ int continueOuter = 0; for(unsigned int j = 1; j < a; j++){ z = POW(z, 2, w); if (z == W) continueOuter = 1; break; if (z == 1) return 0;// Composite } if (continueOuter) {continue;} } return 0; //This is the line you're missing. } return 1;// Probably Prime 

In addition, if the input is even, it will always return as simple as possible, since a is 0. For this you need to add an additional check.

Answer 2

In the inner loop

  for(unsigned int j = 1; j < a; j++){ z = POW(z, 2, w); if (z == W) continue; if (z == 1) return 0;// Composite } 

you must break; instead of continue; when z == W In continue ing, in the next iteration of this loop, if there is one, z will become 1, and the candidate may have been mistakenly declared compound. Here it happens for 17, 41, 73, 89 and 97 among primes less than 100.