Return to the planet of the JavaScript Primes

I finally broke down and decided to do the sensible thing: implement Eratosthenes’ Sieve in JavaScript.

Eratosthenes of Cyrene

Eratosthenes was an over‐achieving Greek polymath. Born in Cyrene in Greek north Africa around 276ʙᴄ, he, among other things, invented the discipline of geography, accurately measured the circumference of the Earth and gave us the leap year. He also came up with his eponymous sieve: a quick way of finding prime numbers.

Sieve it

The sieve works by removing composite (non‐prime) numbers from the set of the natural numbers by seiving out the multiples of each prime, one prime at a time.

Lets start with the numbers from two to twenty:

\[\{2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20\}\]

Two is a prime, so we remove any numbers greater than two that have it as a divisor:

\[ \{2, 3, {\color{red}4}, 5, {\color{red}6}, 7, {\color{red}8}, 9, {\color{red}10}, 11, {\color{red}12}, 13, {\color{red}14}, 15, {\color{red}16}, 17, {\color{red}18}, 19, {\color{red}20}\} \]

The next number that’s not already been sieved out, three, must be a prime, so we repeat the procedure for two and remove any multiples:

\[\{2, 3, 5, 7, {\color{red}9}, 11, 13, {\color{red}15}, 17, 19\}\]

Six, being divisible by two, is already removed, as are twelve and eighteen, so the only numbers we remove this time are nine and fifteen.

Five is the next number still in the list, but we don’t have any numbers left that are divisible by it. In fact, because \(5^2\) is greater than our biggest number, twenty, we know that we can stop looking for primes. If our numbers continued up further we would remove twenty‐five, thirty, thirty‐five, and so on, then move onto multiples of seven. As it stands, our list is now:

\[\{2, 3, 5, 7, 11, 13, 17, 19\}\]

i.e. all the prime numbers under twenty-one.

Solve it

Implementation is rather straight‐forward. We create an array of a set size and step through the algorithm outlined above.

function seive(primeList) {
  for (var p = 2 ; p < primeList.length ; p++) {
    if (primeList[p]) {
      // Only do this if p is marked as a prime

      for (var n = p*2 ; n < primeList.length ; n = n+p) {
        // increase n by p each time so we step through each integer
        // that can be divided by p, and, of course, if it's divisible
        // by p, then it's not prime
        primeList[n] = false;
      }
    }
  }
  return primeList;
}

// 0 and 1 are never prime
var primeList = [false, false, true, true, true, true, true, true,
                 true, true, true, true, true, true, true, true,
                 true, true, true, true];

primeList = seive(primeList);

// 0 → false, 1 → false, 2 → true, 3 → true, 4 → false, 5 → true, &c.

Now we can just do a simple, very fast, look‐up on the array to find out if a number is a prime or not.

The problem with this method is that we have to decide ahead of time what the largest number we can check is; if we want to increase our look‐up table range, we have to go back to the beginning and start again. It’s also rather inefficient if we just want to check one, very large, number as we’d potentially have to calculate every prime up to that number.

Large look‐up tables could also use a significant amount of memory and take a long time to create.

Save it

The solution is to use a hybrid approach. For small numbers we generate a look‐up table, and for larger numbers we use something like the Miller‐Rabin algorithm¹.

var P = (function(lookupSize) {
  var lookup = null;

  function millerRabin(n) {
    function isOdd(n) {
      return n%2===1;
    }

    function rnd(floor, ceil) {
      return floor+Math.floor(Math.random()*(ceil-floor));
    }

    // A custom exponentiation function. Raising numbers to large
    // powers can very quickly take javascript's numbers out of their
    // linear range, but we only need the modulus of the result, and
    // modular arithmetic allows us to cheat.
    function expmod(base, exponent, mod) {
      var result=1;

      while (exponent > 0) {
        if (isOdd(exponent)) {
          result=result*base%mod;
          exponent--;
        }

        base=base*base%mod;
        exponent=exponent/2;
      }
      return result;
    }

    // This check finds whether the number is composite. If it fails
    // to find the number composite it doesn't mean it's definitely
    // prime. It gets run several times for each number we want to
    // check so as to reduce the chances we give a false positive.
    function test(q, s) {
      var a=rnd(1, n-1);
      var apowq=expmod(a, q, n);

      if (apowq===1 || apowq===n-1) return true;

      for (var i=1 ; i < s ; i++) {
        var t=expmod(apowq, 1 << i, n);

        if (t===n-1 || t===1) {
          return true;
        }
      }

      return false;
    }

    for (var q=n-1, s=0 ; !isOdd(q) ; q/=2, s++);

    // Check up to twenty times
    for (var i=0, t=true ; i<20 && t ; i++)
      t=t && test(q, s);

    return t;
  }

  // initialise the array
  function initLookup(size) {
    l = Array(size);
    l[0] = false;
    l[1] = false;

    for (var i = 2 ; i < size ; i++) {
      l[i] = true;
    }

    return l;
  }

  // Seive out the composites
  function seive(primeList) {
    for (var p = 2 ; p*p < primeList.length ; p++) {
      if (primeList[p]) {
        for (var n = p*2 ; n < primeList.length ; n = n+p) {
          primeList[n] = false;
        }
      }
    }
    return primeList;
  }

  // Test divide against some low numbers
  function divisionTest(n) {
    var primes = [2, 3, 5, 7];
    for (var i = 0 ; i < primes.length ; i++) {
      if (n%primes[i] === 0) {
        return false;
      }
    }
    return true;
  }

  function isPrime(n) {
    // Generate lookup when it's first required
    if (lookup === null) {
      lookup = seive(initLookup(lookupSize));
    }

    // For our purposes negative primes are the same as positive
    // primes
    if (n < 0) {
      n = -n;
    }

    if (n < lookupSize) {
      return lookup[n];
    }

    // If the division test doesn't detect a composite, use
    // miller-rabin
    if (divisionTest(n)) {
      return millerRabin(n);
    }

    return false;
  }

  return {
    isPrime: isPrime
  };
})(500000); // Lookup size is half a million

Try it out:

Footnotes:

Because the Miller‐Rabin function squares numbers, and javascript uses sixty‐four bit floats (fifty‐three bit precision), you can’t safely check numbers larger than \(2^{26}\) (67,108,864).