diff options
Diffstat (limited to 'miller_rabin.c')
-rw-r--r-- | miller_rabin.c | 80 |
1 files changed, 80 insertions, 0 deletions
diff --git a/miller_rabin.c b/miller_rabin.c new file mode 100644 index 0000000..3794ee6 --- /dev/null +++ b/miller_rabin.c @@ -0,0 +1,80 @@ +/* Source: + * http://en.literateprograms.org/Miller-Rabin_primality_test_%28C%29 + * + */ + +#include <stdio.h> +#include <stdlib.h> + +#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; +} + |