Scientific Calculator

silentmatt/miller-rabin

By Matthew Crumley

millerRabin and isProbablePrime functions

Tagged: number-theory primes

function() {

var lowprimes = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997];
var lplim = (1 * 2^26) / lowprimes[len lowprimes - 1];

// test primality with certainty >= 1-.5^t
var global:isProbablePrime = function(this, t) {
  var i;
  var x = |this|;
  var res = false;

  if (x <= lowprimes[len lowprimes - 1]) -> (
    i = 0;
    while (not res and i < len lowprimes) -> (
      if (x == lowprimes[i]) -> (
        res = true;
      );
      i += 1
    );
    res
  )
  else if (isEven(x)) -> (
    false
  )
  else (
    i = 1;
    res = true;
    while (res and i < len lowprimes) -> (
      var m = lowprimes[i];
      var j = i + 1;
      while (j < len lowprimes and m < lplim) -> (
        m *= lowprimes[j];
        j += 1
      );
      m = x mod m;
      while (res and i < j) -> (
        if (m mod lowprimes[i] == 0) -> (
          res = false
        );
        i += 1
      )
    );
    if (res) -> (
      millerRabin(x, t)
    )
    else (
      false
    )
  )
};

var shiftToOdd = function(n) {
  var k = 0;
  while (isEven(n)) -> (
    n /= 2;
    k += 1
  );
  [n, k]
};

// (protected) true if probably prime (HAC 4.24, Miller-Rabin)
var global:millerRabin = function(this, t) {
  var n1 = this - 1;
  var n1_k = shiftToOdd(n1);
  var k = n1_k[1];
  var r = n1_k[0];

  if (k <= 0) -> (
    false
  )
  else (
    var undefined;
    t == undefined ? t = 14 : 0;
    t = min((t + (isOdd(t) ? 1 : 2)) / 2, len lowprimes);
    var a;
    var i = 0;
    var res = true;
    while (res and i < t) -> (
      a = lowprimes[exact floor(random() * len lowprimes)];
      var y = modPow(a, r, this);
      if (y != 1 and y != n1) -> (
        var j = 1;
        while (res and j < k and y != n1) -> (
          j += 1;
          y = modPow(y, 2, this);
          if (y == 1) -> (
            res = false
          )
        );
        if (y != n1) -> (
          res = false
        )
      );
      i += 1
    );
    res
  )
};

}();

spam? | offensive?