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>
       
 
       
// ---------- taken from my toolbox ----------

       
 
       
// greatest common divisor

       
template <typename T>
       
T gcd(T a, T b)
       
{
       
  while (a != 0)
       
  {
       
    T c = a;
       
    a = b % a;
       
    b = c;
       
  }
       
  return b;
       
}
       
 
       
// ---------- brute force solution (not used anymore) ----------

       
 
       
// try all combinations of a,b,c

       
unsigned int bruteForce(unsigned int limit)
       
{
       
  // accept solution if left side of the equation of Descartes theorem
       
  // doesn't differ by more than epsilon from the right side
       
  const double Epsilon = 0.0000000000001;
       
 
       
  // sum of all valid a,b,c
       
  unsigned int result = 0;
       
  for (unsigned int a = 1; a <= limit; a++)
       
    for (unsigned int b = 1; b <= a; b++)
       
      for (unsigned int c = 1; c <= b; c++)
       
      {
       
        // curvatures
       
        double k1 = 1.0 / a;
       
        double k2 = 1.0 / b;
       
        double k3 = 1.0 / c;
       
        // both sides of the equation
       
        double leftSide  = (k1 + k2 + k3) * (k1 + k2 + k3);
       
        double rightSide = 2 * (k1*k1 + k2*k2 + k3*k3);
       
        // acceptable error ?
       
        if (fabs(leftSide - rightSide) <  Epsilon)
       
        {
       
          auto sum = a + b + c;
       
          //std::cout << a << "+" << b << "+" << c << "=" << sum << std::endl;
       
          result += sum;
       
        }
       
      }
       
 
       
  return result;
       
}
       
 
       
// ---------- relevant code ----------

       
 
       
// compute triangular number

       
unsigned long long triangle(unsigned int n)
       
{
       
  // sum of 1+2+3+...+n
       
  // see https://en.wikipedia.org/wiki/Triangular_number
       
  unsigned long long nn = n;
       
  return nn * (nn + 1) / 2;
       
}
       
 
       
typedef long long Number;
       
unsigned long long evaluate(unsigned int limit)
       
{
       
  // all multiples of the basic 4+4+1 solution
       
  unsigned long long result = 0;
       
 
       
  // all multiples of a^2, b^2, c^2
       
  // where 1/c2 = (1/a2 + 1/b2) + 2sqrt(1/a2 * 1/b2)
       
  // => sqrt(1/c2) = sqrt(1/a2) + sqrt(1/b2)
       
  // => 1/c = 1/a + 1/b
       
  // => c = a*b / (a+b)
       
  for (unsigned long long a = 1; a*a < limit; a++)
       
    for (unsigned long long b = 1; b <= a; b++)
       
    {
       
      // numerator = a*a * b*b
       
      auto a2 = a*a;
       
      auto b2 = b*b;
       
      auto numerator = a2 * b2;
       
 
       
      // denominator = (a+b)^2
       
      auto denominator = (a+b) * (a+b);
       
 
       
      // c^2 must be an integer
       
      auto c2 = numerator / denominator;
       
      // no remainder left by the previous divison ?
       
      if (c2 * denominator != numerator)
       
        continue;
       
 
       
      // all values must be coprime
       
      if (gcd(gcd(a2, b2), c2) == 1)
       
      {
       
        auto multiples = limit / (unsigned int)a2;
       
        auto sum = (a2 + b2 + c2) * triangle(multiples);
       
        result += sum;
       
      }
       
    }
       
 
       
  return result;
       
}
       
 
       
int main()
       
{
       
  unsigned int limit = 1000000000;
       
  std::cin >> limit;
       
 
       
  //std::cout << bruteForce(limit) << std::endl;
       
  std::cout << evaluate  (limit) << std::endl;
       
  return 0;
       
}

This solution contains 16 empty lines, 27 comments and 2 preprocessor commands.