title: Miller-Rabin Primality test FIPS 186-3 implementation
body: <p>Im trying to implement the Miller-Rabin primality test according to the description in <a href="http://csrc.nist.gov/publications/fips/fips186-3/fips_186-3.pdf" rel="noreferrer">FIPS 186-3</a> C.3.1. No matter what I do, I cannot get it to work. The instructions are pretty specific, and I dont think I missed anything, and yet Im getting <code>true</code> for non-prime values.</p>
  
  <p>What did I do wrong?</p>
  
  <pre><code>template &lt;typename R, typename S, typename T&gt;
  T POW(R base, S exponent, const T mod){
      T result = 1;
      while (exponent){
          if (exponent &amp; 1)
              result = (result * base) % mod;
          exponent &gt;&gt;= 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 &amp; 1)){
          m &gt;&gt;= 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 &lt; iterations; i++){
          uint64_t b = (rand() % (W - 3)) + 2; // 2 &lt;= b &lt;= w - 2
          uint64_t z = POW(b, m, w);
          if ((z == 1) || (z == W))
              continue;
          else
              for(unsigned int j = 1; j &lt; a; j++){
                  z = POW(z, 2, w);
                  if (z == W)
                      continue;
                  if (z == 1)
                      return 0;// Composite
              }
      }
      return 1;// Probably Prime
  }
  </code></pre>
  
  <p>this:</p>
  
  <pre><code>std::cout &lt;&lt; MillerRabin_FIPS186(33) &lt;&lt; std::endl;
  std::cout &lt;&lt; MillerRabin_FIPS186(35) &lt;&lt; std::endl;
  std::cout &lt;&lt; MillerRabin_FIPS186(37) &lt;&lt; std::endl;
  std::cout &lt;&lt; MillerRabin_FIPS186(39) &lt;&lt; std::endl;
  std::cout &lt;&lt; MillerRabin_FIPS186(45) &lt;&lt; std::endl;
  std::cout &lt;&lt; MillerRabin_FIPS186(49) &lt;&lt; std::endl;
  </code></pre>
  
  <p>is giving me: </p>
  
  <pre><code>0
  1
  1
  1
  0
  1
  </code></pre>