A Few Sequences More

| 0 Comments

In the last few posts[1][2] we took a look at sequences and series, from the arithmetic and geometric to the Fibonacci, defined by the rules
\[ \begin{align*} s_0 &= 0\\ s_1 &= 1\\ s_n &= s_{n-1} + s_{n-2} \end{align*} \]
At the risk of being tedious, there are two more sequences that I'd like to cover before we move on to a new topic.

The Primal Stream

We can trivially make a sequence of the prime numbers, those integers greater than one that cannot be wholly divided by any positive integers other than themselves and one, by arranging them in ascending order. For example, its first few terms are
\[ \begin{align*} s_1 &= 2\\ s_2 &= 3\\ s_3 &= 5\\ s_4 &= 7\\ s_5 &= 11 \end{align*} \]
Unfortunately, the terms are anything but trivial to find.

The first thing that we need to do is work out how to test whether or not a number is a prime. The simplest algorithm is to simply test whether or not it can be divided by any integer less than it and greater than one, as done by program 1.

Program 1: Testing For Primes

Whilst this certainly works, it's hardly a model of efficiency. For one thing, if \(a\) wholly divides \(i\) so that
\[ i \div a = b \]
for some integer \(b\), then trivially
\[ i \div b = a \]
It must be the case, therefore, that if \(a\) is greater than \(b\) and
\[ i = a \times b \]
then we would have found that \(b\) is a divisor before we tested \(a\). Consequently, we don't need to test divisors that are greater than the square root of \(i\), as demonstrated by program 2.

Program 2: Stopping The Tests Short

Note that we don't need to consider the case when \(c\) is equal to the square root of \(i\) either, since it'd trivially be a divisor. We can also avoid the cost of the square root operation by comparing the square of \(c\) to \(i\) instead, as shown by program 3.

Program 3: Cutting Out The Square Root

We're still doing far more work than is necessary though; there's no point in checking whether or not \(i\) can be wholly divided by four if we already know that it can't be wholly divided by two, for example. Much better would be to only check whether or not \(i\) can be divided by prime numbers.

Since we're aiming to generate the sequence of primes in ascending order, we can easily do this if we use memoization. By saving the terms that we have already found, we can speed up the code both by never looking for them again and by using them as trial divisors when looking for new ones.
The ak.primeSequence function given in listing 1 uses exactly this strategy.

Listing 1: ak.primeSequence
var primes = [ 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];
var n = primes.length;

ak.primeSequence = function(i) {
 var c, j, p;

 if(i!==ak.floor(i) || i<0) {
  throw new Error('invalid index in ak.primeSequence');
 }
 
 if(i===ak.INFINITY) return i;
 if(i<n) return primes[i];

 c = primes[n-1];
 while(n<=i) {
  c += 2;

  j = 0;
  do {p = primes[++j];}
  while(p*p<c && c%p!==0);

  if(p*p>c) {
   primes.push(c);
   ++n;
  }
 }
 return c;
};

Now, there are a few details of this function that I think are worthy of explanation.

Firstly, we are again following the programming convention of starting the sequence at an index of zero rather than the mathematical convention of starting at an index of one.
Secondly, my expectation is that we'll generally be more interested in earlier terms of the prime sequence than in latter terms and so the primes cache is pre-populated with all of the primes less than one hundred so that each of them will be returned in constant time.
Next, the only way that we can know that we have found the \(i^\mathrm{th}\) prime is if we have already found all but the \(i^\mathrm{th}\). It makes sense, therefore, to add each that we find to the cache whenever we need to iterate up to one that wasn't already in it.
When doing so we can exploit the fact that the only even prime is two by only checking the odd numbers after the last prime in the cache, which we do here by iterating in steps of two and not testing for divisibility by two.
Finally, we only test whether or not a number is divisible by the primes in the cache whose squares are less than it to avoid the redundant testing done by program 3.

Program 4 demonstrates ak.primeSequence by using it to enumerate the first one hundred primes.

Program 4: ak.primeSequence

I'm sure that you'll agree that this is a good deal faster than the algorithm implemented by program 3!

The Prime Number Theorem

If we define a function \(\pi(x)\) as the number of primes less than or equal to \(x\), then it has a limit as \(x\) tends to infinity of
\[ \lim_{x \to \infty} \pi(x) = \frac{x}{\ln x} \]
which we can rewrite in terms of our prime sequence as
\[ \lim_{n \to \infty} \frac{s_n}{n \times \ln s_n} = 1 \]
This is known as the prime number theorem and is demonstrated by program 5 which plots
\[ \frac{s_{2^i}}{2^i \times \ln s_{2^i}} \]
against \(i\).

Program 5: The Prime Number Theorem

Here the first tick on the \(y\) axis is at one and, whilst the convergence to the limit is admittedly rather slow, I think that it's clear that the points get closer and closer to it as \(i\) increases.

Now I'll happily confess that I haven't the foggiest idea how to prove the prime number theorem, but it does at least shed some light upon why ak.primeSequence is faster than the algorithm given in program 3; in the worst case that the number \(c\) that we're testing is in fact prime, program 3 performs approximately \(\sqrt{c}\) trial divisions whereas ak.primeSequence performs approximately \(\sqrt{c}\) divided by the logarithm of \(\sqrt{c}\).

Halton, Who Goes There?

The last type of sequence that I'd like to cover, at least for now, are Halton sequences[3]; a family of low discrepancy, or quasi random, sequences. Unlike pseudo random number generators[4], which are designed to approximate uniformly distributed discrete random variables that choose from a finite set of evenly spaced numbers with equal probability, quasi random sequences are designed to progressively fill intervals as evenly as possible. Whilst these might sound like identical goals, there is a distinction, albeit a subtle one.
For a uniformly distributed random variable, say over an interval of width \(W\), we expect every subinterval of width \(w\) to contain approximately \(\frac{w}{W} \times n\) out of any \(n\) observations, for sufficiently large \(n\).
For a quasi random sequence over an interval of width \(W\), however, we require that every subinterval of width \(w\) contains approximately \(\frac{w}{W} \times n\) out of the first \(n\) terms, for all \(n\).

Halton's trick was to define the terms of a sequence by reversing the digits of their indices and sticking a decimal point in front of them. For example
\[ \begin{align*} s_{1} &= 0.1\\ s_{12} &= 0.21\\ s_{123} &= 0.321 \end{align*} \]
So the sequence begins with the unit interval's tenths, followed by its hundredths, starting with the first after each tenth, then the second and so on, as demonstrated by program 6.

Program 6: One Hundred Terms

Now there's no particular reason why we should only do this with decimal numbers; for each base we can construct a Halton sequence using the same trick of digit reversal. Doing so with their text representations is horribly inefficient however, requiring conversions between numbers, strings and arrays. We can do much better by using arithmetic operations instead, as done by ak.haltonSequence in listing 2
.
Listing 2: ak.haltonSequence
ak.haltonSequence = function(base) {
 if(!isFinite(base) || base!==ak.floor(base) || base<=1) {
  throw new Error('invalid base in ak.haltonSequence');
 }

 return function(i) {
  var bn = 1;
  var x = 0;

  if(i!==ak.floor(i) || i<0) {
   throw new Error('invalid index in ak.haltonSequence');
  }

  if(i===ak.INFINITY) return 1;

  while(i>0) {
   bn *= base;
   x += (i%base)/bn;
   i = ak.floor(i/base);
  }
  return x;
 };
};

Here we're using repeated integer division by the base to iterate through the digits of the index, from the least to the most significant, reading them off using the remainder operation. To reverse the digits we simply divide them by increasing powers of the base, which we're keeping track of in bn, and accumulate the sum of the resulting numbers in x, so that if
\[ i = i_n \, i_{n-1} \, \dots \, i_1 \, i_0 = \sum_{k=0}^n i_k \times b^k \]
where \(b\) is the base and \(\sum\) is the summation sign, then
\[ x = \sum_{k=0}^n i_k \times b^{-(k+1)} = 0. \, i_0 \, i_1 \, \dots \, i_{n-1} \, i_n \]
as required.

Program 7 compares the distribution of the terms of a sequence created with ak.haltonSequence, which is defined in HaltonSequence.js by the way, with that of the output of Math.random by plotting a histogram of the former in green together with one of the latter in blue.

Program 7: Pseudo Versus Quasi Random

This clearly illustrates the distinction between the evenly spread out terms of a quasi random sequence and the equally likely results of a pseudo random number generator and is a fitting conclusion to this post!

References

[1] A Sequence Of Predictable Events, www.thusspakeak.com, 2016.

[2] The Fibonacci Code, www.thusspakeak.com, 2016.

[3] Halton, J., Algorithm 247: Radical-inverse quasi-random point sequence, Communications of the ACM, Vol. 7, No. 12, 1964.

[4] Pseudo Algorithmical, www.thusspakeak.com, 2016.

Leave a comment