cp's OEIS Frontend

A316148 Number of non-congruent solutions of x^2+y^2 == z^2+w^2 (mod n).

%I A316148 #14 Dec 18 2023 01:43:18
%S A316148 1,8,33,96,145,264,385,896,945,1160,1441,3168,2353,3080,4785,7680,
%T A316148 5185,7560,7201,13920,12705,11528,12673,29568,18625,18824,26001,36960,
%U A316148 25201,38280,30721,63488,47553,41480,55825,90720,51985,57608,77649,129920,70561,101640,81313
%N A316148 Number of non-congruent solutions of x^2+y^2 == z^2+w^2 (mod n).
%H A316148 Amiram Eldar, <a href="/A316148/b316148.txt">Table of n, a(n) for n = 1..10000</a>
%H A316148 <a href="/index/Su#sums_of_squares">Index to sequences related to sums of squares</a>.
%F A316148 Multiplicative with a(2^e) = 2^(2e+1)*(2^e-1), a(p^e) = p^(3e)+p^(3e-1)-p^(2e-1) for odd primes p.
%F A316148 Sum_{k=1..n} a(k) ~ c * n^4 / 4, where c = zeta(2)/zeta(3) = 1.368432... (A306633). - _Amiram Eldar_, Dec 18 2023
%p A316148 A316148 := proc(n)
%p A316148     a := 1;
%p A316148     for pe in ifactors(n)[2] do
%p A316148         p := op(1,pe) ;
%p A316148         e := op(2,pe) ;
%p A316148         if p = 2 then
%p A316148             a := a*p^(2*e+1)*(p^e-1) ;
%p A316148         else
%p A316148             a := a*p^(2*e-1)*(p^(e+1)+p^e-1) ;
%p A316148         end if;
%p A316148     end do:
%p A316148     a ;
%p A316148 end proc:
%p A316148 seq(A316148(n),n=1..100) ;
%t A316148 f[2, e_] :=  2^(2*e+1)*(2^e-1); f[p_, e_] := p^(3*e)+p^(3*e-1)-p^(2*e-1); a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[n]; Array[a, 100] (* _Amiram Eldar_, Sep 11 2020 *)
%o A316148 (PARI) a(n) = {my(f = factor(n), p, e); prod(i = 1, #f~, p = f[i,1]; e = f[i,2]; if(p == 2, 2^(2*e+1)*(2^e-1), p^(3*e)+p^(3*e-1)-p^(2*e-1)));} \\ _Amiram Eldar_, Dec 18 2023
%Y A316148 Cf. A096018, A240547, A306633.
%K A316148 nonn,easy,mult
%O A316148 1,2
%A A316148 _R. J. Mathar_, Jun 25 2018