%I #19 Jul 01 2024 10:41:48
%S 0,4,9,16,25,13,49,64,81,29,121,25,169,53,34,256,289,85,361,41,58,125,
%T 529,73,625,173,729,65,841,38,961,1024,130,293,74,97,1369,365,178,89,
%U 1681,62,1849,137,106,533,2209,265,2401,629,298,185,2809,733,146,113
%N a(n) = A008475(n^2).
%H Alois P. Heinz, <a href="/A081403/b081403.txt">Table of n, a(n) for n = 1..10000</a>
%F Additive with a(p^e) = p^(2e).
%e a(1) = 0 since 1 has no prime factor; n = p^2: a(p^2) = p^2; n = 6: a(6) = 4+9 = 13; a(u*w) = a(u)+a(w) if gcd(u,w) = 1; a(21) = a(7)+a(3) = 49+9 = 58; additive with respect of unitary prime divisor decompositions.
%p a:= n-> add(i[1]^i[2], i=ifactors(n^2)[2]):
%p seq(a(n), n=1..60); # _Alois P. Heinz_, Sep 03 2019
%t ffi[x_] := Flatten[FactorInteger[x]] lf[x_] := Length[FactorInteger[x]] ba[x_] := Table[Part[ffi[x], 2*w-1], {w, 1, lf[x]}] ep[x_] := Table[Part[ffi[x], 2*w], {w, 1, lf[x]}] supo[x_] := Apply[Plus, ba[x]^ep[x]] Table[supo[w], {w, 1, 25}]
%o (PARI) f(n) = { my(f=factor(n)); vecsum(vector(#f~, i, f[i, 1]^f[i, 2])); }; \\ A008475
%o a(n) = f(n^2); \\ _Michel Marcus_, Sep 03 2019
%o (Python)
%o from sympy import factorint
%o def A081403(n): return sum(p**(e<<1) for p,e in factorint(n).items()) # _Chai Wah Wu_, Jul 01 2024
%Y Cf. A008475, A000142, A081403, A081404.
%K nonn
%O 1,2
%A _Labos Elemer_, Mar 31 2003