OFFSET
1,4
COMMENTS
LINKS
Daniel Suteu, Table of n, a(n) for n = 1..10000
FORMULA
G.f.: Sum_{k>=1} x^prime(k) * (1 + x^prime(k)) / (1 - x^prime(k))^3. - Ilya Gutkovskiy, Oct 10 2019
a(n) = 1 <=> n is prime. _ Alois P. Heinz, Oct 11 2019
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
EXAMPLE
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.
MAPLE
a:= n-> n^2*add(1/i[1]^2, i=ifactors(n)[2]):
seq(a(n), n=1..70); # Alois P. Heinz, Oct 11 2019
MATHEMATICA
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 *)
PROG
(PARI) a(n) = my(f=factor(n)[, 1]~); sum(k=1, #f, n^2\f[k]^2);
(Magma) [0] cat [n^2*&+[1/p^2:p in PrimeDivisors(n)]:n in [2..70]]; // Marius A. Burtea, Oct 10 2019
CROSSREFS
KEYWORD
nonn
AUTHOR
Daniel Suteu, Nov 25 2018
STATUS
approved