A096018 Number of Pythagorean quadruples mod n; i.e., number of solutions to w^2 + x^2 + y^2 = z^2 mod n.
1, 8, 21, 64, 145, 168, 301, 512, 621, 1160, 1221, 1344, 2353, 2408, 3045, 4096, 5185, 4968, 6517, 9280, 6321, 9768, 11661, 10752, 18625, 18824, 16281, 19264, 25201, 24360, 28861, 32768, 25641, 41480, 43645, 39744, 51985, 52136, 49413, 74240
Offset: 1
Links
- Amiram Eldar, Table of n, a(n) for n = 1..10000 (terms 1..127 from R. J. Mathar, terms 128..2500 from Andrew Howroyd)
- A. H. Hakami, Small zeros of quadratic congruences to a prime power modulus, PhD Thesis (2009), Lemma 4.4.
- A. H. Hakami, Small primitive zeros of quadratic forms mod p^m, 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).
- László Tóth, Counting Solutions of Quadratic Congruences in Several Variables Revisited, J. Int. Seq. 17 (2014), Article 14.11.6.
- Eric Weisstein's World of Mathematics, Pythagorean Quadruple.
- Index to sequences related to sums of squares.
Crossrefs
Programs
-
Maple
A096018 := proc(n) a := 1; for pe in ifactors(n)[2] do p := op(1,pe) ; e := op(2,pe) ; if p = 2 then a := a*p^(3*e) ; elif modp(p,4) = 1 then a := a* p^(2*e-1)*(p^(e+1)+p^e-1) ; else if type(e,'even') then a := a* (p^(3*e)+(p-1)*p^(2*e-1)*(1-p^e)/(1+p)) ; else a := a* (p^(3*e)-(p-1)*p^(2*e-1)*(1+p^e)/(1+p)) ; end if; end if; end do: a ; end proc: seq(A096018(n),n=1..50) ; # R. J. Mathar, Jun 24 2018
-
Mathematica
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}] 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 *) -
PARI
M(n,f)={sum(i=0, n-1, Mod(x^(f(i)%n), x^n-1))} a(n)={polcoeff(lift(M(n, i->i^2)^3 * M(n, i->-(i^2))), 0)} \\ Andrew Howroyd, Jun 23 2018 -
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
Formula
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]
Sum_{k=1..n} a(k) ~ c * n^4 / 4, where c = A334425 * A334426 /(A088539 * A243381) = 0.94532146880744347512... . - Amiram Eldar, Nov 21 2023