%I #34 Oct 18 2022 07:27:18
%S 1,4,12,8,30,48,56,0,108,120,132,96,182,224,360,0,306,432,380,240,672,
%T 528,552,0,750,728,972,448,870,1440,992,0,1584,1224,1680,864,1406,
%U 1520,2184,0,1722,2688,1892,1056,3240,2208,2256,0,2744,3000,3672,1456
%N Number of solutions of x^2 + y^2 + z^2 == -1 (mod n) with x, y, and z in 0..n-1.
%H Andrew Howroyd and David A. Corneth, <a href="/A229179/b229179.txt">Table of n, a(n) for n = 1..10000</a> (first 2500 terms from Andrew Howroyd).
%H László Tóth, <a href="https://cs.uwaterloo.ca/journals/JIS/VOL17/Toth/toth12.html">Counting Solutions of Quadratic Congruences in Several Variables Revisited</a>, J. Int. Seq. 17 (2014), Article 14.11.6.
%F a(8 * n) = 0; for odd prime p, a(p^k) = p^(2 * k - 1) * (p + 1); a(2) = 4, a(4) = 8. - _David A. Corneth_, Jun 24 2018
%F Sum_{k=1..n} a(k) ~ c * n^3, where c = 13/(4*Pi^2) = 0.329293... . - _Amiram Eldar_, Oct 18 2022
%e As 60 = 4 * 3 * 5, a(60) = a(4) * a(3) * a(5) = 8 * (3 * (3 + 1)) * (5 * (5 + 1)) = 8 * 12 * 30 = 2880. - _David A. Corneth_, Jun 24 2018
%t Table[Sum[ If[Mod[a^2 + b^2 + c^2 + 1, n] == 0, 1, 0], {c, 0, n - 1}, {b, 0, n - 1}, {a, 0, n - 1}], {n, 14}]
%t f[p_, e_] := If[p == 2, Which[e == 1, 4, e == 2, 8, e > 2, 0], (p + 1)*p^(2*e - 1)]; a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[n]; Array[a, 50] (* _Amiram Eldar_, Oct 18 2022 *)
%o (PARI) a(n)={my(p=Mod(sum(i=0, n-1, x^(i^2 % n)), x^n-1)); polcoeff(lift(p^3), n-1)} \\ _Andrew Howroyd_, Jun 24 2018
%o (PARI) first(n) = {my(res = vector(n)); forstep(i = 1, n, 2, f = factor(i); res[i] = 1; for(j = 1, #f~, res[i] *= f[j, 1] * (f[j, 1] + 1) * f[j, 1] ^ ((f[j, 2] - 1) << 1)); res); for(k = 1, 2, forstep(i = 1, n >> k, 2, res[i << k] = res[i] << (k+1))); res} \\ _David A. Corneth_, Jun 24 2018
%Y Cf. A086932, A229180.
%K nonn,easy,mult
%O 1,2
%A _José María Grau Ribas_, Sep 29 2013