%I #33 Jun 23 2023 16:59:53
%S 0,0,0,0,1,0,0,4,0,0,6,0,4,4,0,12,1,0,12,4,6,4,12,12,0,12,6,12,12,0,
%T 24,16,0,12,18,12,13,16,12,28,6,0,36,16,12,24,24,24,4,16,30,24,18,12,
%U 36,36,0,28,42,12,36,16,24,52,1,24,48,28,18,24,60,36,12
%N Number of solutions to w^2 + x^2 + y^2 + z^2 = n in positive integers.
%H T. D. Noe, <a href="/A063730/b063730.txt">Table of n, a(n) for n = 0..10000</a>
%H <a href="/index/Su#ssq">Index entries for sequences related to sums of squares</a>
%F G.f.: (Sum_{m>=1} x^(m^2))^4.
%F a(n) = ( A000118(n) - 4*A005875(n) + 6*A004018(n) - 4*A000122(n) + A000007(n) )/16. - _Max Alekseyev_, Sep 29 2012
%F G.f.: ((theta_3(q) - 1)/2)^4, where theta_3() is the Jacobi theta function. - _Ilya Gutkovskiy_, Aug 08 2018
%t r[n_] := Reduce[ w > 0 && x > 0 && y > 0 && z > 0 && w^2 + x^2 + y^2 + z^2 == n, {w, x, y, z}, Integers]; a[n_] := Which[rn = r[n]; rn === False, 0, Head[rn] === Or, Length[rn], True, 1]; Table[a[n], {n, 0, 72}] (* _Jean-François Alcover_, Jul 22 2013 *)
%t a[n_ ] := Length[FindInstance[{n == w^2 + x^2 + y^2 + z^2, w > 0, x > 0, y > 0, z > 0}, {w, x, y, z}, Integers, 10^18]]; (* _Michael Somos_, Jun 23 2023 *)
%o (PARI) seq(n)=Vec((sum(k=1, sqrtint(n), x^(k^2)) + O(x*x^n))^4 + O(x*x^n), -(n+1)) \\ _Andrew Howroyd_, Aug 08 2018
%Y Cf. A063691, A063725.
%Y Column k=4 of A337165.
%K nonn
%O 0,8
%A _N. J. A. Sloane_, Aug 23 2001