Source code in the bc programming language

seed = 1   /* seed of the random number generator */
scale = 0

/* Random number from 0 to 32767. */
define rand() {
  /* Cheap formula (from POSIX) for random numbers of low quality. */
  seed = (seed * 1103515245 + 12345) % 4294967296
  return ((seed / 65536) % 32768)
}

/* Random number in range [from, to]. */
define rangerand(from, to) {
  auto b, h, i, m, n, r

  m = to - from + 1
  h = length(m) / 2 + 1  /* want h iterations of rand() % 100 */
  b = 100 ^ h % m        /* want n >= b */
  while (1) {
    n = 0                /* pick n in range [b, 100 ^ h) */
    for (i = h; i > 0; i--) {
      r = rand()
      while (r < 68) { r = rand(); }  /* loop if the modulo bias */
      n = (n * 100) + (r % 100)       /* append 2 digits to n */
    }
    if (n >= b) { break; }  /* break unless the modulo bias */
  }
  return (from + (n % m))
}



/* n is probably prime? */
define miller_rabin_test(n, k) {
  auto d, r, a, x, s

  if (n <= 3) { return (1); }
  if ((n % 2) == 0) { return (0); }

  /* find s and d so that d * 2^s = n - 1 */
  d = n - 1
  s = 0
  while((d % 2) == 0) {
     d /= 2
     s += 1
  }

  while (k-- > 0) {
    a = rangerand(2, n - 2)
    x = (a ^ d) % n
    if (x != 1) {
      for (r = 0; r < s; r++) {
        if (x == (n - 1)) { break; }
        x = (x * x) % n
      }
      if (x != (n - 1)) {
        return (0)
      }
    }
  }
  return (1)
}

for (i = 1; i < 1000; i++) {
  if (miller_rabin_test(i, 10) == 1) {
    i
  }
}
quit