%I #37 Feb 16 2025 08:32:53

%S 1,8,21,64,145,168,301,512,621,1160,1221,1344,2353,2408,3045,4096,

%T 5185,4968,6517,9280,6321,9768,11661,10752,18625,18824,16281,19264,

%U 25201,24360,28861,32768,25641,41480,43645,39744,51985,52136,49413,74240

%N Number of Pythagorean quadruples mod n; i.e., number of solutions to w^2 + x^2 + y^2 = z^2 mod n.

%H Amiram Eldar, <a href="/A096018/b096018.txt">Table of n, a(n) for n = 1..10000</a> (terms 1..127 from R. J. Mathar, terms 128..2500 from Andrew Howroyd)

%H A. H. Hakami, <a href="https://www.researchgate.net/publication/29868875_Small_zeros_of_quadratic_congruences_to_a_prime_power_modulus">Small zeros of quadratic congruences to a prime power modulus</a>, PhD Thesis (2009), Lemma 4.4.

%H A. H. Hakami, <a href="http://dx.doi.org/10.1007/s11139-014-9614-3">Small primitive zeros of quadratic forms mod p^m</a>, Raman. J. 38 (2015) 189-198, Lemma 2.1 for n=4, det Q=-1, omega_j(y')= p^(m-j)-p^(m-j-1).

%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.

%H Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/PythagoreanQuadruple.html">Pythagorean Quadruple</a>.

%H <a href="/index/Su#sums_of_squares">Index to sequences related to sums of squares</a>.

%F a(n) is multiplicative. For the powers of primes p, there are several cases. For p=2, we have a(2^e) = 2^(3e). For odd primes p with p==1 (mod 4), we have a(p^e) = p^(2*e-1)*(p^(e+1)+p^e-1). For odd primes p with p==3 (mod 4) and even e we have a(p^e) = p^(3*e) +(p-1)*p^(2*e-1)*(1-p^e)/(1+p). For odd primes p == 3 (mod 4) and odd e we have a(p^e) = p^(3*e) -(p-1)*p^(2*e-1)*(1+p^e)/(1+p). [Corrected Jun 24 2018, _R. J. Mathar_]

%F Sum_{k=1..n} a(k) ~ c * n^4 / 4, where c = A334425 * A334426 /(A088539 * A243381) = 0.94532146880744347512... . - _Amiram Eldar_, Nov 21 2023

%p A096018 := proc(n)

%p a := 1;

%p for pe in ifactors(n)[2] do

%p p := op(1,pe) ;

%p e := op(2,pe) ;

%p if p = 2 then

%p a := a*p^(3*e) ;

%p elif modp(p,4) = 1 then

%p a := a* p^(2*e-1)*(p^(e+1)+p^e-1) ;

%p else

%p if type(e,'even') then

%p a := a* (p^(3*e)+(p-1)*p^(2*e-1)*(1-p^e)/(1+p)) ;

%p else

%p a := a* (p^(3*e)-(p-1)*p^(2*e-1)*(1+p^e)/(1+p)) ;

%p end if;

%p end if;

%p end do:

%p a ;

%p end proc:

%p seq(A096018(n),n=1..50) ; # _R. J. Mathar_, Jun 24 2018

%t Table[cnt=0; Do[If[Mod[w^2+x^2+y^2-z^2, n]==0, cnt++ ], {w, 0, n-1}, {x, 0, n-1}, {y, 0, n-1}, {z, 0, n-1}]; cnt, {n, 50}]

%t f[2, e_] := 2^(3*e); f[p_, e_] := If[Mod[p, 4] == 1, p^(2*e - 1)*(p^(e + 1) + p^e - 1), If[EvenQ[e], p^(3*e) + (p - 1)*p^(2*e - 1)*(1 - p^e)/(1 + p), p^(3*e) - (p - 1)*p^(2*e - 1)*(1 + p^e)/(1 + p)]]; a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[n]; Array[a, 100] (* _Amiram Eldar_, Sep 11 2020 *)

%o (PARI)

%o M(n,f)={sum(i=0, n-1, Mod(x^(f(i)%n), x^n-1))}

%o a(n)={polcoeff(lift(M(n, i->i^2)^3 * M(n, i->-(i^2))), 0)} \\ _Andrew Howroyd_, Jun 23 2018

%o (PARI) a(n) = {my(f = factor(n)); prod(i = 1, #f~, p = f[i, 1]; e = f[i, 2]; if(p == 2, 2^(3*e), if(p%4 == 1, p^(2*e-1)*(p^(e+1) + p^e - 1), if(e%2, p^(3*e) - (p - 1)*p^(2*e - 1)*(1 + p^e)/(1 + p), p^(3*e) + (p - 1)*p^(2*e - 1)*(1 - p^e)/(1 + p)))));} \\ _Amiram Eldar_, Nov 21 2023

%Y Cf. A062775 (number of solutions to x^2 + y^2 = z^2 mod n), A240547.

%Y Cf. A088539, A243381, A334425, A334426.

%K mult,nonn,easy

%O 1,2

%A _T. D. Noe_, Jun 15 2004