Prime generator in Python and C I have written the Rabin-Miller test in Python and ported it to C. But somehow, the C version is ~20 times slower on my system, and I am not sure how. Anyway, here is my code.

Python

import random
import time

def rabin_miller(n):
    if n in [2, 3, 5]:
        return True
    if n &lt; 2 or n % 2 == 0:
        return False

    s = 0
    d = n - 1
    while d % 2 == 0:
        s += 1
        d &gt;&gt;= 1

    for k in range(2 ** 4):
        a = random.randrange(2, n - 1)
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        c = False
        for r in range(s - 1):
            x = pow(x, 2, n)
            if x == 1:
                return False
            if x == n - 1:
                break
        else:
            return False
    return True

for i in range(10000):
    rabin_miller(i)


C

#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;
#include &lt;time.h&gt;

long mod_pow(long b, long e, long m){
    if(m == 1)
        return 0;
    long c = 1;
    for(long i = 0; i &lt; e; i++)
        c = (c * b) % m;
    return c;
}

int rabin_miller(long n){
    if(n == 2 || n == 3)
        return 1;
    if(n &lt; 2 || n % 2 == 0)
        return 0;

    long s = 0, d = n - 1, k = 10, a, x, stop;
    while(d % 2 == 0)
        s += 1, d &gt;&gt;= 1;

    srand(time(NULL));
    for(long i = 0; i &lt; k; i++){
        a = (long) (((float) rand() / RAND_MAX) * (n - 3)) + 2;
        x = mod_pow(a, d, n);
        if(x == 1 || x == n - 1)
            continue;
        stop = 1;
        for(long r = 0; r &lt; s - 1; r++){
            x = mod_pow(x, 2, n);
            if(x == 1)
                return 0;
            if(x == n - 1){
                stop = 0;
                break;
            }
        }
        if(stop){
            return 0;
        }
    }
    return 1;
}

int main(int argc, char *argv[]){
    for(long i = 1; i &lt; 10000; i++){
        rabin_miller(i);
    }
    return 0;
}
I believe the problem is that long mod_pow is very inefficient, as it takes e steps to complete, as compared to Python's approach which likely only takes log(e) steps. To fix this, you should use an algorithm similar to the following

function modular_pow(base, exponent, modulus)
if modulus = 1 then return 0
Assert :: (modulus - 1) * (modulus - 1) does not overflow base
result := 1
base := base mod modulus
while exponent &gt; 0
    if (exponent mod 2 == 1):
       result := (result * base) mod modulus
    exponent := exponent &gt;&gt; 1
    base := (base * base) mod modulus
return result