Optimized Longest Collatz Sequence

Question

I'm writing a program to calculate the longest collatz sequence for positive integers below some N. I have optimized it as far as I can - is there anything more I could do to optimize it? Any constructive criticism is welcome :)

A few interesting notes: It is easy to prove (by contradiction or otherwise) that the number that produces the longest chain exists in the second half of the numbers from 1-N (that is, it must be greater than N/2). If we try to run the loop from N/2 to N, the running time becomes substantially slower. Interesting... Also, the lower while loop is an interesting idea, but provides no tangible performance benefits (or losses, for that matter - which is why I kept it. It's an interesting idea).

Thank you!

#include <iostream>
#include <chrono>

using namespace std::chrono;

const int N = 1000000;
long long iCpy;
short ct;
short maxLength = 0;
short vals[N] = {0};

int main() {
  // Time the program:
  auto startTime = steady_clock::now();

  int i, maxIndex = 0;

  // calculate collatz lengths for all indices.
  // keep track of max value:
  for (i = 2; i < N; i++) {
    iCpy = i;
    ct = 0;

    while (iCpy != 1) {
      if (iCpy < N && vals[iCpy]) {
        ct += vals[iCpy];
        break;
      }

      if (iCpy & 1) {
        iCpy *= 3;
        iCpy++;
      } else {
        iCpy >>= 1; // Equivalent to iCpy /= 2
      }

      ct++;
    }

    iCpy = i;
    vals[iCpy] = ct;
    while (iCpy < N/2) {
      iCpy <<= 1; // Equivalent to iCpy *= 2
      vals[iCpy] = ++ct;
    }
    if (ct > maxLength) {
      maxLength = ct;
      maxIndex = iCpy;
    }
  }

  // output the max index and length:
  std::cout << "Max number: " << maxIndex << ", with length " << maxLength << std::endl;

  auto elapsed = duration_cast<milliseconds>(steady_clock::now() - startTime);
  std::cout << "Program execution took " << elapsed.count() << "ms\n";
  return 0;
}

Answer 1

Here are some suggestions:

• Don't use global variables.

• Instead of using a single big function, break it down to logical parts.

• const int N = 1000000; should be constexpr instead of const.

• Instead of = {0}, use {}. The latter makes it clear that you are value initializing all elements.

• i should be declared in the loop.

• Post-increment operators i++ should be changed to pre-increment ++i in a discarded-value expression.

• iCpy & 1 should be iCpy % 2 == 0. iCpy <<= 1; // Equivalent to iCpy *= 2 should be iCpy *= 2. Optimization doesn't happen in this way — the compiler knows such things much better.

• Use << '\n' instead of << std::endl when you don't need the flushing behavior.

Answer 2

The lower loop (which fills all power-of-two multiples of the outer index) is doing unnecessary work when i is even. If we find that vals[i] was already filled in by a previous iteration, then so was vals[2*i] and so on.

I improved run-time by about 5% simply by adding at the start of the outer loop:

    if (vals[i]) {
        continue;
    }

---

It's usual to report the lowest starting value of those that give rise to the same length chain. The current code doesn't do that, because of the second loop. To give the same results as seen on OEIS, then we need to change the logic to update maxIndex and maxLength only after computing Collatz(i), and before the *=2 loop.

---

Finally, we don't want to be measuring printing time as well as computation, so move the setting of elapsed before the printing of max number.

---

Improved version

I upped the range to 100 million so that I could get meaningful numbers of milliseconds on my machine. This consistently gives under 1750 ms, compared to slightly over 2000 ms with the original. It also gives the lowest index for a given count, consistent with the results shown on the Collatz entry in Wikipedia.

I've also applied changes suggested in other reviews

#include <chrono>
#include <iostream>

const unsigned int N = 100'000'000;
unsigned short vals[N] = {};

int main() {
    using namespace std::chrono;
    // Time the program:
    auto const startTime = steady_clock::now();

    unsigned int maxIndex = 0;
    unsigned int maxLength = 0;

    // calculate Collatz length for each index, and memoize
    for (unsigned int i = 2;  i < N;  ++i) {
        unsigned long iCpy = i;
        unsigned int ct = 0;

        while (iCpy > 1) {
            if (iCpy < N && vals[iCpy]) {
                ct += vals[iCpy];
                break;
            }

            iCpy = iCpy % 2
                ? 3 * iCpy + 1
                : iCpy / 2;
            ++ct;
        }
        if (ct > maxLength) {
            maxLength = ct;
            maxIndex = i;
        }

        // Pre-fill power-of-two multiples of this result
        for (iCpy = i;  iCpy < N && !vals[i];  iCpy *= 2) {
            vals[iCpy] = ct++;
        }
    }

    auto elapsed = duration_cast<milliseconds>(steady_clock::now() - startTime);

    // output the max index and length:
    std::cout << "Max number: " << maxIndex << ", with length " << maxLength <<  '\n';
    std::cout << "Program execution took " << elapsed.count() << "ms\n";
    return 0;
}

In case you wonder what happened to the continue mentioned in my first paragraph, it has to move to after the maxLength check, and it seems to optimise better when incorporated into the existing for loop condition, perhaps due to branch prediction.

Answer 3

You are more interested in the performance, and on style there are sufficient mentions.

The Collatz sequence can be reduced:

  if (iCpy & 1) {
    iCpy *= 3;
    ++iCpy;
    ++ct;
  }
  iCpy >>= 1;
  ++ct;

while (iCpy != 1) {
   if (iCpy & 1) {
    iCpy *= 3;
    ++iCpy;
    ++ct;
  }
  iCpy >>= 1;
  ++ct;
}

As a loop's intermediate step is not cached, instead first using a recursive function and on coming back update the result cache could be faster. Especially when the compiler eliminates the recursion.

int collatz(int n) {
    if (n <= 1) {
        return 0;
    }
    int st = n & 1;
    if (st) {
        n = 3*n + 1;
    }
    n /= 2;
    return 1 + st + collatz(n);
}

int collatz(int n) {
    if (n <= 1) {
        return 0;
    }
    int res = from_cache(n);
    if (res) {
        return res;
    }
    int n0 = n;

    int st = n & 1;
    if (st) {
        n = 3*n + 1;
    }
    n /= 2;
    res = collatz(n);
    res += 1 + st;
    to_cache(n0, res);
    return res;
}

Answer 4

I reduced the time of the algorithm by ~36% by reducing the algorithm and improving the Cache. Tested for N=500,000,000.

First lets look at the properties of the Collatz sequence

1) If the number is even, then we will perform a right-shift until the value is odd.

2) If the number is odd, then it will take one iteration of n=3*n+1 to become even.

Proof

Let's rewrite our equation n = 3*n + 1 = (2*n) + (n+1) = (n<<1) + (n+1).

The first term (n<<1) is guaranteed to be even, because a left shift will pad a 0 to the lsb.

The second term (n+1) is guaranteed to be even, because an odd number plus one is even.

The sum of two even terms is even.

As a result, n = 3*n + 1 is guaranteed to create an even n (when n is odd).

3) n = 3*n+1 for odd n will never be 1.

Proof

See #2. Result must be even.

Now let's rewrite our Collatz sequence

// if even, right shift until odd
while ((n & 1) == 0) {
    n >>= 1;
    ++len;
}
while ( n != 1 ) {
    // odd (guaranteed)
    n = 3*n + 1;
    // even (guaranteed)
    n >>= 1;
    // if even, right shift until odd
    while ((n & 1) == 0) {
        n >>= 1;
        ++len;
    }
    // account for 2 guaranteed steps (see above)
    len += 2;
}

Now let's improve the Cache

From our properties (see above), we notice that the Collatz sequence will frequently alternate between positive and negative numbers. The maximum number of odd numbers that can be seen before seeing an even number is one. The maximum number of even numbers that can be seen before seeing an odd number depends on N, but it is relatively small due to right shifting (max 31 for an unsigned int). For this reason, I have decided to only store odd numbers in the Cache, which will reduce our Cache to half the size.

Additionally, I have decided to store values in the cache as I encounter them--rather than set only the initial values at the end of each sequence. In order to achieve this, we will need to store the values in a preliminary buffer, and then, at the end of the Collatz sequence, add them to the cache and calculate the length by taking the difference between the length that we first encountered the number, and the final length of the sequence.

Since we are storing the values we encounter rather than the initial values. We can now improve our algorithm by starting at N/2 as you mentioned, because we will still cache values for the range [2, N/2).

Here is the resulting function

void collatz(unsigned N, unsigned& r_index, unsigned& r_length) {
    unsigned *cache, n0, len, i=0, lenMax=0;
    unsigned long long n;

    // cache size
    const unsigned cacheMaxN = N;
    const unsigned cacheSize = (cacheMaxN/2) + 1; // store all odd from (1,cacheMaxN) exclusive

    // cache
    cache = new unsigned[cacheSize];
    memset(cache, 0, cacheSize*sizeof(unsigned));

    // cache buffer (for temporarily holding cache values)
    struct CacheObject {
        unsigned n;
        unsigned n0_to_n_len;
    };
    std::vector<CacheObject> buffer;
    buffer.reserve(128);
    CacheObject cObj;

    for (n0 = N/2; n0 < N; ++n0 ) {
        // initial
        n = n0;
        len = 0;
        // collatz
        // if even, right shift until odd
        while ((n & 1) == 0) {
            n >>= 1;
            ++len;
        }
        while ( n != 1 ) {
            // odd (guaranteed), check cache
            if ( n < cacheMaxN ) {
                if ( cache[ n>>1 ] ) {
                    // use cache value
                    len += cache[ n>>1 ];
                    break;
                } else {
                    // store cache value
                    // store to cache as difference between current length and final length
                    // needs to resolve final length first, so store in queue
                    cObj.n = n;
                    cObj.n0_to_n_len = len;
                    buffer.push_back( cObj );
                }
            }
            // odd  (guaranteed)
            n = 3*n + 1;
            // even (guaranteed)
            n >>= 1;
            // if even, right shift until odd
            while ((n & 1) == 0) {
                n >>= 1;
                ++len;
            }
            // account for 2 guaranteed steps (see above)
            len += 2;
        }
        // buffer to cache
        if ( buffer.size() != 0 ) {
            for (CacheObject o : buffer) {
                cache[ o.n>>1 ] = len - o.n0_to_n_len;
            }
            buffer.clear();
        }
        // set max len
        if ( len > lenMax ) {
            lenMax = len;
            i = n0;
        }
    }

    r_index = i;
    r_length = lenMax;

    delete [] cache;
}