Find the number of $n$ such that $n!$ is a sum of three squares

Brute-force, but compact:

DeleteCases[Table[{k, PowersRepresentations[k!, 3, 2]}, {k, 10}], {___, 0, ___}, {3}]

{{1, {}}, {2, {}}, {3, {{1, 1, 2}}}, {4, {{2, 2, 4}}}, {5, {{2, 4, 10}}},
 {6, {{8, 16, 20}}}, {7, {{4, 20, 68}, {12, 36, 60}, {20, 44, 52}}},
 {8, {{8, 16, 200}, {8, 80, 184}, {40, 88, 176}, {72, 120, 144}, {80, 104, 152}}},
 {9, {{8, 304, 520}, {16, 200, 568}, {24, 48, 600}, {24, 240, 552}, {40, 104, 592},
      {40, 272, 536}, {80, 184, 568}, {80, 344, 488}, {120, 264, 528}, {136, 272, 520},
      {152, 400, 424}, {176, 248, 520}, {184, 368, 440}, {200, 328, 464},
      {216, 360, 432}, {240, 312, 456}, {248, 376, 400}}}, {10, {}}}

For larger factorials, it might take quite longer; e.g. Length[PowersRepresentations[11!, 3, 2]] == 126.

You can Try

Table[set = {x, y, z} /. 
          NSolve[{n! == x^2 + y^2 + z^2, x > 0, y > 0, z > 0},{x, y, z},Integers];
set = Union[Sort /@ set];
Join[{n}, set],
 {n, 3,8}]

This will give you how you can express a factorial as a sum of three squares. Now you can use further conditions (like $x\neq y \neq z$) with Select or Cases to filter them.

If you only care about counting and enumerating, why not use Legendre's three-square theorem?

SetAttributes[SumOf3SquaresQ, Listable];

SumOf3SquaresQ[n_] := Mod[n/4^IntegerExponent[n, 4], 8] != 7

FactorialSumOf3SquaresPi[x_] := Total[Boole[SumOf3SquaresQ[Range[x]!]]]

FactorialSumOf3SquaresPi[10000]

8746

Or test your conjecture by subtracting $7x/8$ and dividing by $x^{2/3}$:

With[{x = 1000},
  ListLinePlot[
    (Accumulate[Boole[SumOf3SquaresQ[Range[x]!]]] - 7 Range[x]/8)/Range[x]^(2/3)
  ]
]

Or list such numbers:

Pick[Range[50], SumOf3SquaresQ[Range[50]!]]

{1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 
  26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 
  45, 46, 47, 50}

Edit

If we keep a running tally (instead of computing things separately), we can get results much quicker:

FactorialSumOf3SquaresPiFast = Compile[{{x, _Integer}},
  Module[{count = 0, m, offset = 1, prod = 1},
    Do[
      m = n;
      While[BitAnd[m, 3] == 0, m = Floor[m/4]];
      If[EvenQ[m], m = Floor[m/2]; offset++];
      prod = BitAnd[prod * m, 7];
      If[prod != 7 - Boole[EvenQ[offset]], count++],
      {n, 1, x}
    ];
    count
  ],
  CompilationTarget -> "C",
  RuntimeOptions -> "Speed"
];

Compare timings of old vs new:

FactorialSumOf3SquaresPi[10000] // AbsoluteTiming

{1.57493, 8746}

FactorialSumOf3SquaresPiFast[10000] // AbsoluteTiming

{0.000495, 8746}

Here's timings on large $x$:

FactorialSumOf3SquaresPiFast[10^7] // AbsoluteTiming

{0.488944, 8751280}

FactorialSumOf3SquaresPiFast[10^8] // AbsoluteTiming

{4.8907, 87502143}

FactorialSumOf3SquaresPiFast[10^9] // AbsoluteTiming

{49.8406, 875001736}