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 <algorithm>
       
#include <cmath>
       
 
       
const double PI = 3.1415926535897932;
       
 
       
int main()
       
{
       
  // example was n = 200, problem asks for n = 10000
       
  unsigned int order = 10000;
       
  std::cin >> order;
       
 
       
  // find maximum value such that e^(maximum/n) - 1 <= PI
       
  // ... and while at it: store all values in f
       
  std::vector<double> f;
       
  unsigned int maximum = 0;
       
  while (true)
       
  {
       
    auto current = exp(maximum / double(order)) - 1;
       
    // too large ?
       
    if (current > PI)
       
      break;
       
 
       
    f.push_back(current);
       
    maximum++;
       
  }
       
 
       
  // generate all pairs f[i] + f[j] <= PI (with i <= j)
       
  std::vector<double> pairs;
       
  pairs.reserve(maximum * maximum / 10); // some heuristic to avoid too frequent allocations
       
  for (size_t i = 0; i < maximum; i++)
       
    for (size_t j = i; j < maximum; j++)
       
    {
       
      // abort, sum exceeds 3.1415
       
      if (f[i] + f[j] > PI)
       
        break;
       
 
       
      pairs.push_back(f[i] + f[j]);
       
    }
       
 
       
  // sort all values
       
  //std::sort(pairs.begin(), pairs.end());
       
  // a bit faster: compare double as 64 bit integers => works on CPU, not FPU
       
  // only possible because all values are positive and there are no NaN or other special values
       
  std::sort(pairs.begin(), pairs.end(), [] (double a, double b)
       
                                      {
       
                                        // cast from 64 bit double to 64 bit integer
       
                                        auto aa = (const long long*) &a;
       
                                        auto bb = (const long long*) &b;
       
                                        return *aa < *bb;
       
                                      });
       
 
       
  // optimal value for f_a + f_b
       
  size_t left  = 0;
       
  // optimal value for f_c + f_d
       
  size_t right = 0;
       
  // minimize error |f_a + f_b - PI|
       
  double minError  = PI;
       
  for (size_t i = 0; i < pairs.size(); i++)
       
  {
       
    // binary search for the best match
       
    auto current = pairs[i];
       
    auto need    = PI - current;
       
    if (need < current) // same as need < PI/2
       
      break;
       
 
       
    // this value and its predecessor are candidates
       
    auto match = std::upper_bound(pairs.begin(), pairs.end(), need); // actually "sums.begin() + i" works, too
       
 
       
    // "match" is a value slightly too large
       
    auto error = fabs(need - *match);
       
    if (error < minError)
       
    {
       
      minError = error;
       
      left  = i;
       
      right = std::distance(pairs.begin(), match);
       
    }
       
 
       
    // now "match" is a value slightly too small (or an exact match => spoiler alert: there is no perfect match)
       
    match--;
       
    error = fabs(need - *match);
       
    if (error < minError)
       
    {
       
      minError = error;
       
      left  = i;
       
      right = std::distance(pairs.begin(), match);
       
    }
       
  }
       
 
       
  // resolve the left pair f_a and f_b
       
  unsigned int result = 0;
       
  bool resolvedLeft = false;
       
  for (size_t a = 0; a < maximum && !resolvedLeft; a++)
       
    for (size_t b = a; b < maximum; b++)
       
      if (pairs[left]  == f[a] + f[b])
       
      {
       
        result += a*a + b*b;
       
        resolvedLeft  = true;
       
        break;
       
      }
       
 
       
  // resolve the left pair f_c and f_d
       
  bool resolvedRight = false;
       
  for (size_t c = 0; c < maximum && !resolvedRight; c++)
       
    for (size_t d = c; d < maximum; d++)
       
      if (pairs[right] == f[c] + f[d])
       
      {
       
        result += c*c + d*d;
       
        resolvedRight = true;
       
        break;
       
      }
       
 
       
  std::cout << result << std::endl;
       
  return 0;
       
}

This solution contains 14 empty lines, 20 comments and 4 preprocessor commands.