/* Source: * http://en.literateprograms.org/Miller-Rabin_primality_test_%28C%29 * */ #include #include #define COMPOSITE 0 #define PRIME 1 unsigned short modular_exponent_16(unsigned short base, unsigned short power, unsigned short modulus) { unsigned long result = 1; int i; for (i = 15; i >= 0; i--) { result = (result*result) % modulus; if (power & (1 << i)) result = (result*base) % modulus; } return (unsigned short) result; /* will not truncate since modulus is a unsigned short */ } int miller_rabin_pass_16(unsigned short a, unsigned short n) { unsigned short a_to_power, s, d, i; s = 0; d = n - 1; while ((d % 2) == 0) { d /= 2; s++; } a_to_power = modular_exponent_16(a, d, n); if (a_to_power == 1) return PRIME; for (i = 0; i < s-1; i++) { if (a_to_power == n-1) return PRIME; a_to_power = modular_exponent_16(a_to_power, 2, n); } if (a_to_power == n-1) return PRIME; return COMPOSITE; } int miller_rabin_16(unsigned short n) { if (n <= 1) return COMPOSITE; if (n == 2) return PRIME; if (miller_rabin_pass_16( 2, n) == PRIME && (n <= 7 || miller_rabin_pass_16( 7, n) == PRIME) && (n <= 61 || miller_rabin_pass_16(61, n) == PRIME)) return PRIME; else return COMPOSITE; } int main(int argc, char* argv[]) { unsigned short n = (unsigned short) atoi(argv[1]); puts(miller_rabin_16(n) == PRIME ? "PRIME" : "COMPOSITE"); return 0; }