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 <deque>
       
#include <algorithm>
       
 
       
// ---------- powmod from my toolbox ----------

       
 
       
// return (base^exponent) % modulo for 32-bit values, no need for mulmod

       
unsigned int powmod(unsigned int base, unsigned int exponent, unsigned int modulo)
       
{
       
  unsigned int result = 1;
       
  while (exponent > 0)
       
  {
       
    // fast exponentation:
       
    // odd exponent ? a^b = a*a^(b-1)
       
    if (exponent & 1)
       
      result = (result * (unsigned long long)base) % modulo;
       
 
       
    // even exponent ? a^b = (a*a)^(b/2)
       
    base = (base * (unsigned long long)base) % modulo;
       
    exponent >>= 1;
       
  }
       
  return result;
       
}
       
 
       
// ---------- problem-specific code ----------

       
 
       
// a 2D point

       
struct Location
       
{
       
  // coordinates
       
  unsigned int x;
       
  unsigned int y;
       
 
       
  // create new point
       
  Location(unsigned int x_ = 0, unsigned y_ = 0)
       
  : x(x_), y(y_)
       
  {}
       
 
       
  // return point at (2^index % modulo, 3^index % modulo)
       
  static Location generate(unsigned int index, unsigned int modulo)
       
  {
       
    auto x = powmod(2, index, modulo);
       
    auto y = powmod(3, index, modulo);
       
    return Location(x, y);
       
  }
       
 
       
  // for sorting
       
  bool operator<(const Location& other) const
       
  {
       
    if (x != other.x)
       
      return x < other.x;
       
    else
       
      return y < other.y;
       
  }
       
  // for sorting
       
  bool operator==(const Location& other) const
       
  {
       
    return x == other.x && y == other.y;
       
  }
       
};
       
 
       
// return length of longest increasing subsequence

       
template <typename T>
       
unsigned int getLIS(T first, T last) // begin/end of a container
       
{
       
  // no data ?
       
  if (first == last)
       
    return 0;
       
 
       
  // seqBack[i] contains the smallest last element of a subsequence with length i + 1 (plus one because buffer is zero-indexed)
       
  std::vector<typename T::value_type> seqBack;
       
  seqBack.reserve(1 << 16);
       
 
       
  // initialize subsequence with length = 1 with the first element
       
  seqBack.push_back(*first++);
       
  while (first != last)
       
  {
       
    // locate shortest sequence where the last element >= current element
       
    auto current = *first++;
       
    auto pos = std::upper_bound(seqBack.begin(), seqBack.end(), current);
       
    // current element is the largest element observed so far ? => it extends the longest sequence
       
    if (pos == seqBack.end())
       
      seqBack.push_back(current);
       
    else // current element is smaller for that sequence length => "optimize"
       
      *pos = current; // could skip this step if *pos == current
       
  }
       
 
       
  return seqBack.size();
       
}
       
 
       
// compute S(n)

       
unsigned int minPathLength(unsigned int n)
       
{
       
  // generate all stations
       
  std::deque<Location> stations;
       
  //if (n > 10000000)
       
    //stations.reserve(n * 1.05);
       
 
       
  size_t lastUnsortedSize = 0;
       
  size_t lastSortedSize   = 0;
       
 
       
  for (unsigned int i = 0; i <= 2 * n; i++)
       
  {
       
    auto current = Location::generate(i, n);
       
    // the route always starts at (0,0) so discard any stations placed there
       
    if (current.x == 0 && current.y == 0)
       
      continue;
       
 
       
    stations.push_back(current);
       
 
       
    // 0x100000 is the sweet spot:
       
    // - if I remove duplicates more often, then I spent too much time sorting
       
    // - if I remove duplicates less often, then memory consumption grows (and sort/unique slow down as well)
       
    if (i % 0x100000 == 0 && i > 0)
       
    {
       
      // prune duplicates
       
      std::sort(stations.begin(), stations.end());
       
      auto garbage = std::unique(stations.begin(), stations.end());
       
      size_t unsortedSize = stations.size();
       
      size_t sortedSize   = std::distance(stations.begin(), garbage);
       
      stations.erase(garbage, stations.end());
       
 
       
      // same sizes as in last pass ? => assume that we entered a cycle/loop
       
      if (unsortedSize == lastUnsortedSize &&
       
          sortedSize   == lastSortedSize)
       
        break;
       
 
       
      lastUnsortedSize = unsortedSize;
       
      lastSortedSize   = sortedSize;
       
    }
       
  }
       
 
       
  // prune duplicates
       
  std::sort(stations.begin(), stations.end());
       
  auto garbage = std::unique(stations.begin(), stations.end());
       
  stations.erase(garbage, stations.end());
       
 
       
  // reduce locations to their y-components
       
  std::vector<unsigned int> onlyY;
       
  // simple but memory consuming
       
  //for (auto station : stations)
       
  //  onlyY.push_back(station.y);
       
  // a bit slower but keeps me below 256 MByte
       
  while (!stations.empty())
       
  {
       
    onlyY.push_back(stations.begin()->y);
       
    stations.erase(stations.begin());
       
  }
       
 
       
  // find shortest route
       
  return getLIS(onlyY.begin(), onlyY.end());
       
}
       
 
       
int main()
       
{
       
  unsigned int limit = 30;
       
  std::cin >> limit;
       
 
       
  unsigned int sum = 0;
       
//#define PARALLEL

       
#ifdef  PARALLEL
       
#pragma omp parallel for reduction(+:sum) schedule(dynamic)
       
#endif
       
  for (int k = limit; k >= 1; k--) // run loop in reverse to optimize load for OpenMP
       
    // evaluate k^5
       
    sum += minPathLength(k*k*k*k*k);
       
 
       
  std::cout << sum << std::endl;
       
  return 0;
       
}

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