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>
       
 
       
// cache gcd(a,b)

       
std::vector<std::vector<unsigned int>> pointsOnEdge;
       
 
       
// 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;
       
}
       
 
       
// count lattice points with Pick's theorem
       
// A(a, 0), B(b, 0), C(-c, 0), D(-d, 0)

       
unsigned int countLatticePoints(unsigned int a, unsigned int b, unsigned int c, unsigned int d)
       
{
       
  // Pick's theorem: A = i + b/2 - 1
       
  // where i is the area inside the quadrilateral
       
  // and b is the number of lattice points on the boundary
       
  // https://en.wikipedia.org/wiki/Pick%27s_theorem
       
 
       
  // total area of the bounding box
       
  auto rectangle = (a + c) * (b + d);
       
  // only half of the area belongs to the quadrilateral
       
  auto inside = rectangle / 2;
       
 
       
  // points on the edges of the quadrilateral
       
  auto boundary = pointsOnEdge[a][b] + pointsOnEdge[b][c] + pointsOnEdge[c][d] + pointsOnEdge[a][d];
       
  // the four corners are not inside the bounding box (they are on its edges)
       
  boundary -= 4;
       
 
       
  // Pick's formula
       
  return inside - (boundary / 2) - 1;
       
}
       
 
       
int main()
       
{
       
  unsigned int limit = 100;
       
  std::cin >> limit;
       
 
       
  // precompute number of lattice points lying exactly on the sides/edges
       
  pointsOnEdge.resize(limit + 1);
       
  for (auto& x : pointsOnEdge)
       
    x.resize(limit + 1);
       
  for (unsigned int a = 1; a <= limit; a++)
       
    for (unsigned int b = 1; b <= limit; b++)
       
      // note: only one endpoint is included for each edge
       
      //       e.g. a=4 b=2 gcd(4,2)=2
       
      //       => has actually 3 lattice points:
       
      //          both endpoints plus point (3,1)
       
      //       but each endpoint is part of two edges therefore it would be counted twice
       
      //       that means that keeping only one endpoint is okay and much simpler
       
      pointsOnEdge[a][b] = pointsOnEdge[b][a] = gcd(a, b);
       
 
       
  // precompute square numbers
       
  // note: std::vector<bool> would be sufficient but std::vector<char> is about 20% faster
       
  std::vector<unsigned char> isSquare(countLatticePoints(limit, limit, limit, limit) + 1, false);
       
  for (size_t i = 1; i*i < isSquare.size(); i++)
       
    isSquare[i*i] = true;
       
 
       
  // plain enumeration of all combinations, 100^4 = 10^8 iterations
       
  auto count = 0;
       
  for (unsigned int a = 1; a <= limit; a++)
       
    for (unsigned int b = 1; b <= limit; b++)
       
      for (unsigned int c = 1; c <= limit; c++)
       
        for (unsigned int d = 1; d <= limit; d++)
       
          // count only perfect squares
       
          if (isSquare[countLatticePoints(a, b, c, d)])
       
            count++;
       
 
       
  // display result
       
  std::cout << count << std::endl;
       
  return 0;
       
}

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