%I #26 Oct 18 2022 07:27:09
%S 1,4,6,24,30,24,42,96,54,120,110,144,182,168,180,384,306,216,342,720,
%T 252,440,506,576,750,728,486,1008,870,720,930,1536,660,1224,1260,1296,
%U 1406,1368,1092,2880,1722,1008,1806,2640,1620,2024,2162,2304,2058,3000
%N Number of solutions to x^2 + y^2 + z^2 = 1 mod n.
%H Michael De Vlieger, <a href="/A087784/b087784.txt">Table of n, a(n) for n = 1..10000</a>
%H László Tóth, <a href="http://arxiv.org/abs/1404.4214">Counting solutions of quadratic congruences in several variables revisited</a>, arXiv preprint arXiv:1404.4214 [math.NT], 2014.
%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>, Journal of Integer Sequences, 17 (2014), Article 14.11.6.
%F a(n) = n^2 * (3/2 if 4|n) * Product_{primes p == 1 mod 4 dividing n} (1+1/p) * Product_{primes p == 3 mod 4 dividing n} (1-1/p). - Bjorn Poonen, Dec 09 2003
%F Sum_{k=1..n} a(k) ~ c * n^3 + O(n^2 * log(n)), where c = 36*G/Pi^4 = 0.338518..., where G is Catalan's constant (A006752) (Tóth, 2014). - _Amiram Eldar_, Oct 18 2022
%t Table[With[{f = FactorInteger[n][[All, 1]]}, Apply[Times, Map[1 + 1/# &, Select[f, Mod[#, 4] == 1 &]]] Apply[Times, Map[1 - 1/# &, Select[f, Mod[#, 4] == 3 &]]] (1 + Boole[Divisible[n, 4]]/2) n^2] - Boole[n == 1], {n, 50}] (* _Michael De Vlieger_, Feb 15 2018 *)
%o (PARI) a(n) = {my(f=factor(n)); if ((n % 4), 1, 3/2)*n^2*prod(k=1, #f~, p = f[k,1]; m = p % 4; if (m==1, 1+1/p, if (m==3, 1-1/p, 1)));} \\ _Michel Marcus_, Feb 14 2018
%Y Cf. A006752, A060968, A087687, A208895, A264083.
%K nonn,mult
%O 1,2
%A Yuval Dekel (dekelyuval(AT)hotmail.com), Oct 06 2003
%E More terms from _David Wasserman_, Jun 17 2005