Sunday 6 October 2013

Computing the Collatz sequence

Problem 14 from ProjectEuler concerns the Collatz sequence. From the problem's description:

The following iterative sequence is defined for the set of positive integers:
n → n/2 (n is even)
n → 3n + 1 (n is odd)

This shows how to compute the Collatz sequence for n. As an example, the Collatz sequence for 13 is:
13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1

The problem is to find which number, under one million, has the longest Collatz sequence. I took the following approach:

From the above example, the length of the Collatz for 13 is equal to the length for 40's  sequence plus one. In general, the length of the Colaltz sequence for n is equal to the length for second term of that sequence plus one. This can be implemented neatly with the recursive function:
uint getCollatzSeriesLength(uint value)
{
  uint seriesLength;

  if (value == 1) {
    seriesLength = 1;
  }
  else {
    if (value%2 == 0) {
      seriesLength getCollatzSeriesLength(value/2) + 1;
    }
    else {
      seriesLength = getCollatzSeriesLength(value*3+1) +1;
    }
  }
  
  return seriesLength;
}

The recursion is terminated by returning a length of 1 for a value 1. Using the above function for it takes about 580 milliseconds to calculate the Collatz series length for every integer from 1 to 1,000,000 on my machine.

I then made the following optimisation. Consider again the Collatz series of 13:
13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1
If the program starts at 1 then when it comes to 13 it has already calculated the series length for 10. If this value was cached then it could be returned without recursing all the way to the end of the series. Therefore when when the algorithm calculates a series length it should cache that, and the function should check that cache before calculating the series length for the next term. Note an upshot of this is that the series length's for 40 and 20 will be calculated in the course of calculating the series length of 13. The updated function looks like this (changes highlighted):
uint getCollatzSeriesLength(uint value)
{
  uint seriesLength;

  if ((value <= noTerms) && (lengthsCache[value] != 0)) {
    seriesLength = lengthsCache[value];
  }
  else if (value == 1) {
    seriesLength = 1;
  }
  else {
    if (value%2 == 0) {
      seriesLength = getCollatzSeriesLength(value/2) + 1;
    }
    else {
      seriesLength = getCollatzSeriesLength(value*3+1) +1;
    }
    if (value <= noTerms) {
      lengthsCache[value] = seriesLength;
    }
  }

  return seriesLength;
}

Here noTerms is the number of terms to be checked (one million) and lengthsCached is a vector initialised with one million and one zeros (so it's effectively indexed from one). This optimisation makes around a 21x difference in terms of speed, as it now takes about 28 milliseconds to run the program.

I made a final optimisation by considering that when the series length has been calculated for some value, n, then the series length for 2n will be one greater. Using this the following ammendment was made to the code:
uint getCollatzSeriesLength(uint value)
{
  uint seriesLength;

  if ((value <= noTerms) && (lengthsCache[value] != 0)) {
    seriesLength = lengthsCache[value];
  }
  else if (value == 1) {
    seriesLength = 1;
  }
  else {
    if (value%2 == 0) {
      seriesLength = getCollatzSeriesLength(value/2) + 1;
    }
    else {
      seriesLength = getCollatzSeriesLength(value*3+1) +1;
    }
    auto additionalSteps = 0;
    for (auto i=value; i<=noTerms; i*=2) {
      lengthsCache[value] = seriesLength additionalSteps++;
    }
  }

  return seriesLength;
}

This brings the program's execution time down by a further factor of ~2.3x to around 12 milliseconds on my machine.

No comments: