%I #41 Mar 03 2023 10:41:38
%S 0,1,1,4,1,13,1,16,9,29,1,52,1,53,34,64,1,117,1,116,58,125,1,208,25,
%T 173,81,212,1,361,1,256,130,293,74,468,1,365,178,464,1,673,1,500,306,
%U 533,1,832,49,725,298,692,1,1053,146,848,370,845,1,1444,1,965,522
%N a(n) = n^2 * Sum_{p|n} 1/p^2, where p are primes dividing n.
%C Generalized formula is f(n,m) = n^m * Sum_{p|n} 1/p^m, where f(n,0) = A001221(n) and f(n,1) = A069359(n).
%H Daniel Suteu, <a href="/A322078/b322078.txt">Table of n, a(n) for n = 1..10000</a>
%F Sum_{k=1..n} a(k) ~ A085541 * A000330(n).
%F G.f.: Sum_{k>=1} x^prime(k) * (1 + x^prime(k)) / (1 - x^prime(k))^3. - _Ilya Gutkovskiy_, Oct 10 2019
%F a(n) = 1 <=> n is prime. _ _Alois P. Heinz_, Oct 11 2019
%F Dirichlet g.f.: zeta(s-2)*primezeta(s). This follows because Sum_{n>=1} a(n)/n^s = Sum_{n>=1} (n^2/n^s) Sum_{p|n} 1/p^2. Since n = p*j, rewrite the sum as Sum_{p} Sum_{j>=1} 1/(p^2*(p*j)^(s-2)) = Sum_{p} 1/p^s Sum_{j>=1} 1/j^(s-2) = zeta(s-2)*primezeta(s). The result generalizes to higher powers of p. - _Michael Shamos_, Mar 02 2023
%e a(40) = 464 because the prime factors of 40 are 2 and 5, so we have 40^2 * (1/2^2 + 1/5^2) = 464.
%p a:= n-> n^2*add(1/i[1]^2, i=ifactors(n)[2]):
%p seq(a(n), n=1..70); # _Alois P. Heinz_, Oct 11 2019
%t f[p_, e_] := 1/p^2; a[n_] := If[n==1, 0, n^2*Plus@@f@@@FactorInteger[n]]; Array[a, 60] (* _Amiram Eldar_, Nov 26 2018 *)
%o (PARI) a(n) = my(f=factor(n)[,1]~); sum(k=1, #f, n^2\f[k]^2);
%o (Magma) [0] cat [n^2*&+[1/p^2:p in PrimeDivisors(n)]:n in [2..70]]; // _Marius A. Burtea_, Oct 10 2019
%Y Cf. A000040, A085541.
%Y Sequences of the form n^k * Sum_{p|n, p prime} 1/p^k for k = 0..10: A001221 (k=0), A069359 (k=1), this sequence (k=2), A351242 (k=3), A351244 (k=4), A351245 (k=5), A351246 (k=6), A351247 (k=7), A351248 (k=8), A351249 (k=9), A351262 (k=10).
%K nonn
%O 1,4
%A _Daniel Suteu_, Nov 25 2018