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 <cmath>
       
 
       
// true if n is a perfect square

       
template <typename T>
       
bool isSquare(T n)
       
{
       
  T root = sqrt(n);
       
  return root * root == n;
       
}
       
 
       
// iterate over all remainders and quotients, it's okay for limit = 100000 but too slow for 10^12

       
void bruteForce(unsigned int limit)
       
{
       
  unsigned int sum = 0;
       
 
       
  for (unsigned int remainder = 1; remainder < limit; remainder++)
       
    for (unsigned int quotient = remainder + 2; quotient < limit; quotient++)
       
    {
       
      // d / r = q / d, therefore d^2 = r * q
       
      auto divisorSquared = remainder * quotient;
       
      unsigned int divisor = sqrt(divisorSquared);
       
 
       
      // find n
       
      auto n = divisor * quotient + remainder;
       
      // too large ?
       
      if (n > limit)
       
        break;
       
      // not a perfect square ?
       
      if (!isSquare(n))
       
        continue;
       
 
       
      // most d^2 are actually NOT a perfect square, therefore all computations above yield garbage
       
      if (!isSquare(divisorSquared))
       
        continue;
       
 
       
      // yes, found another combination
       
      //std::cout << n << " = " << divisor << " * " << quotient << " + " << remainder << std::endl;
       
      sum += n;
       
    }
       
 
       
  std::cout << sum << std::endl;
       
}
       
 
       
int main()
       
{
       
  unsigned long long limit = 1000000000000ULL;
       
  std::cin >> limit;
       
 
       
  //bruteForce(limit); return 0;
       
 
       
  unsigned long long sum = 0;
       
 
       
  // remainder < divisor <= quotient
       
  for (unsigned long long divisor = 1; divisor * divisor <= limit; divisor++)
       
  {
       
    // it seems there is at most one solution per divisor
       
    bool found = false;
       
 
       
    // all possible perfect squares ...
       
    for (unsigned long long y = sqrt(divisor); y >= 1 && !found; y--)
       
      // ... with all possible factors ...
       
      for (unsigned long long factor = 1; factor * y * y <= divisor; factor++)
       
      {
       
        // ... generate all possible remainders
       
        auto remainder = factor * y * y;
       
 
       
        // n = d^3 / r + r
       
        auto n = divisor * divisor * divisor / remainder + remainder;
       
 
       
        // heuristic (MIGHT FAIL !) to skip very large n
       
        if (n > 2*limit)
       
          break;
       
        // ignore if beyond limit
       
        if (n > limit)
       
          continue;
       
 
       
        // must be a perfect square (I moved this check in front of the next check because that turned out to be faster)
       
        if (!isSquare(n))
       
          continue;
       
 
       
        // there's a high likelihood that d^3 isn't a multiple of r and therefore d^3/r isn't an integer
       
        auto quotient = n / divisor;
       
        // cross-check whether all variables produce a valid result
       
        if (divisor * quotient + remainder != n)
       
          continue;
       
 
       
        sum += n;
       
 
       
        // it seems there is at most one solution per divisor
       
        found = true;
       
        break;
       
        // note: I kept track of various remainders (to avoid repetitions for the same divisor)
       
        //       aborting after a valid solution eliminates the need
       
      }
       
  }
       
 
       
  // display result
       
  std::cout << sum << std::endl;
       
  return 0;
       
}

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