#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const long long int MAX = 1e18;
// witness for Miller-Rabin primality test
const int a[9] = {2, 3, 5, 7, 11, 13, 17, 19, 23};

// calculate x * y module p for big number (avoid overflow)
long long int mod_mul(long long int x, long long int y, 
                      long long int p) {
    long long int res = 0;
    x %= p, y %= p;
    while (y > 0) {
        if(y & 1) {
            res += x;
            if(res >= p) res -= p;
        }
        x <<= 1, y >>= 1;
        if(x >= p) x -= p;
    }
    return res;
}

// calculate x^y module p
long long int mod_pow(long long int x, long long int y, 
                      long long int p) {
    long long int res = 1;
    while (y > 0) {
        if (y & 1) res = mod_mul(res, x, p);
        x = mod_mul(x, x, p);
        y >>= 1;
    }
    return res;
}

// determine whether n is prime number
int is_prime(long long int n) {
    int i, j, s, pass;
    long long int d, x;
    if (n < 2) return 0;
    for (i = 0; i < 9; ++i) {
       if (n % a[i] == 0){
            if (n == a[i]) return 1;
            else return 0;
        }
    }
    // n - 1 = 2^s * d (d is odd number)
    d = n - 1, s = 0;
    while (d & 1 == 0) {
        d >>= 1;
        s++;
    }
    for (i = 0; i < 9; ++i) {
        x = mod_pow(a[i], d, n);
        if (x == 1 || x == n - 1) continue;
        pass = 0;
        for (j = 1; j < s; ++j) {
            x = mod_mul(x, x, n);
            if (x == 1) return 0;
            if (x == n - 1) {
                pass = 1;
         printf("3:x = %d \n",x);
                break;
            }
        }
        if (pass == 1) continue;
        return 0;
    }
    return 1;
}

int old_test(long long int n) {
    long long int i, j;
    j = (int) sqrt((double)n);    
    for (i = 2; i <= j; ++i) 
       if (n % i == 0) return 0;
    return 1;
}

/* 実験回数の上限 */
#define MAX_NUM 100000


int main(int argc, char *argv[]) {
    int i, j, first_prime;
    long long int n, ct=0;
    
    n = 0;
    srand((unsigned)time(NULL));

    while (ct++ <MAX_NUM){
        n = rand();
        if (is_prime(n) !=  old_test(n)) {
            printf("%lld　はエラー\n",n);
//            printf("is_prime = %d\n", is_prime(n));
//            printf("old_test = %d\n", old_test(n));
        }
        if (ct % 100 == 0) printf("%d\n", ct);
    }
    return 0;
}