OFFSET
1,4
COMMENTS
The number of these divisors is given by A323308.
LINKS
Amiram Eldar, Table of n, a(n) for n = 1..10000
Vaclav Kotesovec, Plot of Sum_{k=1..n} a(k) / n^(3/2) for n = 1..10^8.
FORMULA
Multiplicative with a(p) = 1 and a(p^e) = p^e + 1 for e > 1.
Dirichlet g.f.: zeta(s)*zeta(s-1)*Product_{p prime} (1 - p^(1-s) + p^(2-2*s) - p^(2-3*s)).
From Vaclav Kotesovec, Feb 18 2023: (Start)
Dirichlet g.f.: zeta(s) * zeta(2*s-2) * Product_{primes p} (1 - p^(3-4*s) - p^(2-3*s) + p^(3-3*s)).
Sum_{k=1..n} a(k) ~ c * zeta(3/2) * n^(3/2) / 3, where c = Product_{primes p} (1 + 1/p^(3/2) - 1/p^(5/2) - 1/p^3) = 1.48039182258752809541724060173644... (End)
a(n) = A034448(A057521(n)) (the sum of unitary divisors of the powerful part of n). - Amiram Eldar, Dec 12 2023
MATHEMATICA
f[p_, e_] := If[e == 1, 1, p^e + 1]; a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[n]; Array[a, 100]
PROG
(PARI) a(n) = {my(f = factor(n)); prod(i = 1, #f~, if(f[i, 2] == 1, 1, f[i, 1]^f[i, 2] + 1)); }
(PARI) for(n=1, 100, print1(direuler(p=2, n, (1 - p^3*X^4 - p^2*X^3 + p^3*X^3) / ((1 - X) * (1 - p^2*X^2)))[n], ", ")) \\ Vaclav Kotesovec, Feb 18 2023
CROSSREFS
KEYWORD
nonn,mult,easy
AUTHOR
Amiram Eldar, Feb 18 2023
STATUS
approved