title: Miller-Rabin test: bug in my code
body: <p>I've written a <a href="http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test" rel="nofollow">Miller-Rabin primality test</a> based on the following pseudo code:</p>
  
  <pre><code>Input: n &gt; 2, an odd integer to be tested for primality;
         k, a parameter that determines the accuracy of the test
  Output: composite if n is composite, otherwise probably prime
  write n − 1 as 2s·d with d odd by factoring powers of 2 from n − 1
  LOOP: repeat k times:
     pick a randomly in the range [2, n − 1]
     x ← ad mod n
     if x = 1 or x = n − 1 then do next LOOP
     for r = 1 .. s − 1
        x ← x2 mod n
        if x = 1 then return composite
        if x = n − 1 then do next LOOP
     return composite
  return probably prime
  </code></pre>
  
  <p>The code I have rarely gets past 31 (if I put it in a loop to test numbers from 2 to 100). There must be something wrong but I can't see what it is.</p>
  
  <pre><code>bool isProbablePrime(ulong n, int k) {
      if (n &lt; 2 || n % 2 == 0) 
          return n == 2;
  
      ulong d = n - 1;
      ulong s = 0;
      while (d % 2 == 0) {
          d /= 2;
          s++;
      }
      assert(2 ^^ s * d == n - 1); 
      outer:
      foreach (_; 0 .. k) {
          ulong a = uniform(2, n);
          ulong x = (a ^^ d) % n;
          if (x == 1 || x == n - 1)
              continue;
          foreach (__; 1 .. s) {
              x = (x ^^ 2) % n;
              if (x == 1) return false;
              if (x == n - 1) continue outer;
          }
          return false;
      }
      return true;
  }
  </code></pre>
  
  <p>I've also tried the variant</p>
  
  <pre><code>    ...
  
      foreach (__; 1 .. s) {
          x = (x ^^ 2) % n;
          if (x == 1) return false;
          if (x == n - 1) continue outer;
      }
      if ( x !=  n - 1) return false;  // this is different
  
      ...
  </code></pre>
  
  <p>I have a different version of the test that works correctly but it uses modpow. I'd like to have a version that stays closer to the pseudo code that's part of the <a href="http://rosettacode.org/wiki/Miller-Rabin_primality_test" rel="nofollow">rossetta.org task description</a>.</p>
  
  <p><strong>Edit</strong>: Re: overflow problem. I suspected something like that. I'm still puzzled why the Ruby version doesn't have that problem. It probably handles it differently under the hood.
  If I use BigInt, the code does work, but becomes a lot slower than when I use modpow. So I guess I can't get away from that. It's a pity Phobos doesn't have a modpow built-in, or I must have overlooked it.</p>
  
  <pre><code>ulong x = ((BigInt(a) ^^ d) % BigInt(n)).toLong();
  </code></pre>