title: Implementing Miller-Rabin in C
body: <p>I'm trying to implement the <a href="https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test" rel="nofollow noreferrer">Miller-Rabin primality test</a> in C99, but I'm coming across some problems getting it to work. I crafted a small test-set to verify whether or not the implementation works, here's how I'm checking for primes </p>
  
  <pre><code>int main() {
      int foo[11] = {0, 1, 2, 3, 4, 7, 28, 73, 125, 991, 1000};
      for (int i = 0; i &lt; 11; i++) {
          printf("%s; ", isprime(foo[i], 5000) ? "Yes" : "No");
      }
      return 0;
  }
  </code></pre>
  
  <p>From the numbers listed, the expected output would be</p>
  
  <blockquote>
    <p>No; No; Yes; Yes; No; Yes; No; Yes; No; Yes; No;</p>
  </blockquote>
  
  <p>However, as implemented , the output I get is the following:</p>
  
  <blockquote>
    <p>No; No; Yes; Yes; No; Yes; No; No; No; No; No; </p>
  </blockquote>
  
  <p>Here's how I wrote the algorithm</p>
  
  <pre><code>int randint (int low, int up){
      return rand() % (++up - low)+low;
  }
  
  int modpow(int a, int b, int m) {
      int c = 1;
      while (b) {
          if (b &amp; 1) {
              c *= a;
          }
          b &gt;&gt;= 1;
          a *= a;
      }
      return c % m;
  }
  
  bool witness(int a, int s, int d, int n) {
      int x = modpow(a,d,n);
      if(x == 1) return true;
      for(int i = 0; i&lt; s-1; i++){
          if(x == n-1) return true;
          x = modpow(x,2,n);
      }
      return (x == n- 1);
  }
  
  bool isprime(int x, int j) {
      if (x == 2) {
          return true;
      }
      if (!(x &amp; 1) || x &lt;= 1) {
          return false;
      }
      int a = 0;
      int s = 0;
      int d = x - 1;
  
      while (!d&amp;1){
          d &gt;&gt;=1;
          s+=1;
      }
      for(int i = 0; i &lt; j; i++){
          a = randint(2, x-1);
          if(!witness(a,s,d,x)){
              return false;
          }
      }
  
      return true;
  }
  </code></pre>
  
  <p>What am I doing wrong? Why is the test failing for "large" primes, but working for very small ones? How may I fix this?</p>