Documentation autogenerated by HWEB

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

typedef long long int_t;

int_t pow_(int_t b, int_t e);
void do_error(const char *fmat, ...);

int main(int argc, char **argv)
{
    int_t n, r, s, a, j;
    if (argc <= 1) do_error("Usage: rabinmiller <possible-prime>");

     

Initialize n. Determine r and s such that

n = 2	r	s+1

    n = strtol(argv[1], NULL, 10);
    if (n % 2 == 0) do_error("n can't be even!");
    s = n-1;
    for (r=0; (s % 2) == 0; ++r) s >>= 1;

     

Now, for each a with 1 ≤ a ≤ n-1, we say that n passes the test for that a if and only if a^s ≡ 1(modn) or a^2js ≡ -1(modn) for some 0 ≤ j ≤ r-1. A prime will pass the tests for every a.

    for (a = 1; a < n; ++a) {
        if (pow_(a, s) % n == 1) goto pass_test;
        for (j = 0; j < r; ++j)
          if (pow_(a, (1 << j)*s) % n == n-1) goto pass_test;
        goto fail_test;
      pass_test:
        continue;
    }

    printf("%ld passes the Rabin-Miller primality test!\n"
           "It may be a prime.\n", (long)n);
    return 0;

  fail_test:
    printf("%ld fails the Rabin-Miller test.\n"
           "It is not a prime.\n", (long)n);
    return 0;
}


int_t pow_(int_t b, int_t e)
{
    if (e < 0) return 0;
    else if (e == 0) return 1;
    else if (e == 1) return b;
    else if (e == 2) return b*b;
    else return pow_(b, e/2) * pow_(b, e - e/2);
}


void do_error(const char *fmat, ...)
{
    va_list ap;
    va_start(ap, fmat);
    vfprintf(stderr, fmat, ap);
    fprintf(stderr, "\n");
    va_end(ap);
    exit(EXIT_FAILURE);
}