My code

… was written in C++11 and can be compiled with G++, Clang++, Visual C++. You can download it, too. Or just jump to my GitHub repository.

       
#include <iostream>
       
#include <vector>
       
#include <cmath>
       
 
       
// find the asymptotic value, quite often off by a few numbers (first four digits tend to be okay)
       
// note: unused code, but helped me a lot during debugging

       
double fast(unsigned long long n)
       
{
       
  // golden ratio
       
  const auto phi = (1 + sqrt(5.0)) / 2;
       
  const auto constant = pow(phi, 2 - phi);
       
  // see https://en.wikipedia.org/wiki/Golomb_sequence
       
  return constant * pow(n, phi - 1);
       
}
       
 
       
int main()
       
{
       
  unsigned int limit = 1000000;
       
  std::cin >> limit;
       
 
       
  // (10^6)^3 = 10^18
       
  auto cubicLimit = (unsigned long long) limit * limit * limit;
       
 
       
  // index 0 is not used, golomb(1) = 1
       
  std::vector<unsigned long long> golomb = { 0, 1 };
       
 
       
#define LOW_MEMORY
       
#ifdef LOW_MEMORY
       
  // precompute golomb[i] = G(i)
       
  // stop when 1*G(1) + 2*G(2) + 3*G(3) + ... + i*G(i) >= 10^18
       
  unsigned long long products = 1;
       
  for (unsigned long long i = 2; products < cubicLimit; i++)
       
  {
       
    // https://en.wikipedia.org/wiki/Golomb_sequence
       
    auto current = 1 + golomb[i - golomb[golomb[i - 1]]];
       
    golomb.push_back(current);
       
    products += current * i;
       
  }
       
 
       
#else
       
 
       
  // precompute golomb[i] = G(i)
       
  // and sums[i]     =   G(1) +   G(2) +   G(3) + ... +   G(i)
       
  // and products[i] = 1*G(1) + 2*G(2) + 3*G(3) + ... + i*G(i)
       
  std::vector<unsigned long long> sums     = { 0, 1 };
       
  std::vector<unsigned long long> products = { 0, 1 };
       
  // stop when products[i] >= 10^18
       
  for (unsigned long long i = 2; products.back() < cubicLimit; i++)
       
  {
       
    auto current = 1 + golomb[i - golomb[golomb[i - 1]]];
       
 
       
    golomb  .push_back(current);
       
    sums    .push_back(current     + sums.back());
       
    products.push_back(current * i + products.back());
       
  }
       
#endif
       
 
       
  // will contain the result
       
  unsigned long long sum = 0;
       
 
       
#ifdef LOW_MEMORY
       
  unsigned long long lastSums     = 0;
       
  unsigned long long sums         = 1;
       
 
       
  unsigned long long lastProducts = 0;
       
                     products     = 1;
       
#endif
       
 
       
  // find products[index - 1] < i <= products[index]
       
  auto index = 1;
       
  for (unsigned long long i = 1; i < limit; i++)
       
  {
       
    // n = i^3
       
    auto n = i * i * i;
       
 
       
    // find products[index - 1] < i <= products[index]
       
#ifdef LOW_MEMORY
       
    while (products < n)
       
    {
       
      index++;
       
      lastSums     = sums;
       
      sums        += golomb[index];
       
      lastProducts = products;
       
      products    += golomb[index] * index;
       
    }
       
#else
       
    while (products[index] < n)
       
      index++;
       
#endif
       
    // note: n will be in ascending order, therefore I re-use index from previous iterations
       
    //       in most cases it will be already the correct value
       
 
       
    // find linear interpolation between products[index - 1] and products[index]
       
#ifdef LOW_MEMORY
       
    auto from = lastProducts;
       
    auto to   = products;
       
#else
       
    auto from = products[index - 1];
       
    auto to   = products[index];
       
#endif
       
    auto ratio = (n - from) / double(to - from);
       
 
       
    // and apply it to sums[index - 1] and sums[index]
       
#ifdef LOW_MEMORY
       
    auto low  = lastSums;
       
    auto high = sums;
       
#else
       
    auto low  = sums[index - 1];
       
    auto high = sums[index];
       
#endif
       
 
       
    // round up
       
    auto offset = ceil((high - low) * ratio);
       
    // note: convert to integer as soon as possible to avoid losing digits due to double's limited precision
       
    auto result = (unsigned long long)offset + low;
       
 
       
    // finished another number ...
       
    sum += result;
       
  }
       
 
       
  // solved another problem !
       
  std::cout << sum << std::endl;
       
  return 0;
       
}

This solution contains 18 empty lines, 25 comments and 18 preprocessor commands.