|
%I #22 Nov 21 2023 05:23:06
%S 1,1,3,4,1,3,13,4,9,1,1,12,25,13,3,16,1,9,37,4,39,1,1,12,25,25,27,52,
%T 1,3,61,16,3,1,13,36,73,37,75,4,1,39,85,4,9,1,1,48,133,25,3,100,1,27,
%U 1,52,111,1,1,12,121,61,117,64,25,3,133,4,3,13
%N Number of solutions to x^2 + xy + y^2 == 0 (mod n).
%H Andrew Howroyd, <a href="/A087694/b087694.txt">Table of n, a(n) for n = 1..10000</a>
%F Multiplicative with a(3^e) = 3^e, a(p^e) = ((p-1)*e+p)*p^(e-1) if p mod 3 = 1, a(p^e) = p^(2*floor(e/2)) if p mod 3 = 2. - _Vladeta Jovovic_, Sep 27 2003
%F Sum_{k=1..n} a(k) ~ c * n^2 / 2, where c = A073010/A086724 = 0.77383581325017004332... . - _Amiram Eldar_, Nov 21 2023
%p A087694 := proc(n) option remember; local pf,p,f,e ; if n = 1 then 1; else pf := ifactors(n)[2] ; if nops(pf) = 1 then f := op(1,pf) ; p := op(1,f) ; e := op(2,f) ; if p = 3 then n ; elif p mod 3 =1 then ((p-1)*e+p)*p^(e-1) ; else p^(2*floor(e/2)) ; end if; else mul(procname(op(1,p)^op(2,p)),p=pf) ; end if; end if; end proc:
%p seq(A087694(n),n=1..70) ; # _R. J. Mathar_, Jan 07 2011
%t a[n_] := If[n==1, 1, Product[{p, e} = pe; Which[p==3, 3^e, Mod[p, 3] == 2, (p^2)^Quotient[e, 2], True, ((p-1) e + p) p^(e-1)], {pe, FactorInteger[n] }]];
%t a /@ Range[1, 100] (* _Jean-François Alcover_, Sep 20 2019, from PARI *)
%o (PARI) a(n)={my(f=factor(n)); prod(i=1, #f~, my(p=f[i,1], e=f[i,2]); if(p==3, 3^e, if(p%3==2, (p^2)^(e\2), ((p-1)*e+p)*p^(e-1))))} \\ _Andrew Howroyd_, Jul 09 2018
%Y Cf. A000086, A073010, A086724.
%K mult,nonn,easy
%O 1,3
%A Yuval Dekel (dekelyuval(AT)hotmail.com), Sep 27 2003
|
