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>
       
 
       
const unsigned int SliceSize = 1000*1000;
       
 
       
// bitmasks of for a^2 + k * b^2

       
const unsigned char One   = 1 << 0; // k = 1
       
const unsigned char Two   = 1 << 1; // k = 2
       
const unsigned char Three = 1 << 2; // k = 3
       
const unsigned char Seven = 1 << 3; // k = 7
       
// every matching number must have all four bits set

       
const unsigned char All   = One | Two | Three | Seven;
       
 
       
int main()
       
{
       
  unsigned int limit = 2*1000*1000*1000;
       
  std::cin >> limit;
       
 
       
  // how many numbers match the criteria
       
  unsigned int count = 0;
       
 
       
  // contains bitmasks for the current slice
       
  std::vector<unsigned char> used(SliceSize, 0);
       
 
       
  // store start values of b for each factor k=1,2,3,7
       
  unsigned int maxA = sqrt(limit);
       
  std::vector<unsigned int> b1(maxA + 1, 1);
       
  std::vector<unsigned int> b2(maxA + 1, 1);
       
  std::vector<unsigned int> b3(maxA + 1, 1);
       
  std::vector<unsigned int> b7(maxA + 1, 1);
       
 
       
  // split into slices
       
  unsigned int from = 0;
       
  while (from < limit)
       
  {
       
    // size of current slice
       
    unsigned int to = from + SliceSize;
       
    if (to > limit)
       
      to = limit;
       
 
       
    // process all a and b where a^2+kb^2 is within the current slice
       
    for (unsigned int a = 1; a*a + b1[a]*b1[a] < to; a++)
       
    {
       
      // a^2 + b^2
       
      unsigned int b = b1[a];
       
      for (; a*a +   b*b < to; b++)
       
        used[a*a +   b*b - from] |= One;
       
      b1[a] = b;
       
 
       
      // a^2 + 2b^2
       
      b = b2[a];
       
      for (; a*a + 2*b*b < to; b++)
       
        used[a*a + 2*b*b - from] |= Two;
       
      b2[a] = b;
       
 
       
      // a^2 + 3b^2
       
      b = b3[a];
       
      for (; a*a + 3*b*b < to; b++)
       
        used[a*a + 3*b*b - from] |= Three;
       
      b3[a] = b;
       
 
       
      // a^2 + 7b^2
       
      b = b7[a];
       
      for (; a*a + 7*b*b < to; b++)
       
        used[a*a + 7*b*b - from] |= Seven;
       
      b7[a] = b;
       
    }
       
 
       
    // count all matching numbers
       
    for (auto& x : used)
       
    {
       
      if (x == All)
       
        count++;
       
 
       
      // reset for next iteration
       
      x = 0;
       
    }
       
 
       
    from = to;
       
  }
       
 
       
  std::cout << count << std::endl;
       
  return 0;
       
}

This solution contains 15 empty lines, 14 comments and 3 preprocessor commands.