# Optimized Longest Collatz Sequence

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:

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;
}

• using unsigned types should speed it up a bit Oct 11, 2019 at 5:28
• It's interesting that you report that the lower loop has no effect on your timings. I measured about 20× difference with it (more visible when I set N to 100'000'000, as the result was too small at just 1'000'000). Oct 15, 2019 at 8:27

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.

• The compiler will optimize i++ to ++i in discarded-value expressions. Oct 11, 2019 at 12:22
• @Kyy13, that's true for integers, but with preincrement, then readers don't need to double-check what type i is. (i.e. ++i is always right, and i++ is sometimes right, depending on the type of i. So let's consistently use the one that's always right.) Oct 11, 2019 at 12:57
• @Kyy13 In addition to what Toby Speight said, ++i just means “increment i” whereas i++ means “copy the value, increment, and then give me the original value” (just like std::exchange). In the same way you don’t want to write int a = static_cast<int>('\008') * 42LL / 21.0 - true; in place of int a = 15; even though compiler optimizations should compile them the same way. Oct 11, 2019 at 13:19
• @TobySpeight I get the same disassembly with i++ and ++i regardless of using floating or integral types when the return value is not needed. Can you explain what you mean? Oct 11, 2019 at 13:59
• @Kyy13, I'm thinking of class types with overloaded ++ operators. When the operator isn't inlineable, your compiler has much less ability to optimise. Oct 11, 2019 at 14:02

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:

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.

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;
}


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;
}