summaryrefslogtreecommitdiffstats
path: root/miller_rabin.c
diff options
context:
space:
mode:
Diffstat (limited to 'miller_rabin.c')
-rw-r--r--miller_rabin.c80
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;
+}
+