A180114 a(n) = sigma(A001694(n)), sum of divisors of the powerful number A001694(n).

1, 7, 15, 13, 31, 31, 40, 63, 91, 57, 127, 195, 121, 217, 280, 133, 156, 255, 403, 183, 399, 465, 600, 403, 364, 511, 819, 307, 847, 400, 381, 855, 961, 1240, 741, 931, 1092, 1023, 553, 1651, 781, 1815, 1240, 1281, 1093, 1767, 1953, 871, 2520, 2821, 993, 1995

Offset: 1

Author

Walter Kehowski, Aug 10 2010

Keywords

• easy
• nonn

Examples

			Sigma(2^2) = 7, sigma(2^3) = 15, sigma(3^2) = 13.
		

Links

• Amiram Eldar, Table of n, a(n) for n = 1..10000
• Rafael Jakimczuk and Matilde Lalín, Asymptotics of sums of divisor functions over sequences with restricted factorization structure, Notes on Number Theory and Discrete Mathematics, Vol. 28, No. 4 (2022), pp. 617-634, eq. (7).

Crossrefs

Cf. A001694, A000203, A090699, A180090, A180097, A180098, A362984.

Programs

Maple

emin := proc(n::posint) local L; if n>1 then L:=ifactors(n)[2]; L:=map(z->z[2],L); min(L) else 0 fi end: L:=[]: for w to 1 do for n from 1 to 144 do sn:=sigma(n); if emin(n)>1 then L:=[op(L),sn]; print(n,ifactor(n),sn,ifactor(sn)) fi; od; od;

Mathematica

pwfQ[n_] := n == 1 || Min[Last /@ FactorInteger[n]] > 1; DivisorSigma[1, Select[ Range@ 1000, pwfQ]] (* Giovanni Resta, Feb 06 2018 *)

PARI

lista(kmax) = for(k = 1, kmax, if(k==1 || vecmin(factor(k)[, 2]) > 1, print1(sigma(k), ", "))); \\ Amiram Eldar, May 12 2023

Python

from itertools import count, islice
from math import prod
from sympy import factorint
def A180114_gen(): # generator of terms
    for n in count(1):
        f = factorint(n)
        if all(e>1 for e in f.values()):
            yield prod((p**(e+1)-1)//(p-1) for p,e in f.items())
A180114_list = list(islice(A180114_gen(),20)) # Chai Wah Wu, May 21 2023

Python

from math import isqrt
from sympy import mobius, integer_nthroot, divisor_sigma
def A180114(n):
    def squarefreepi(n):
        return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
    def bisection(f, kmin=0, kmax=1):
        while f(kmax) > kmax: kmax <<= 1
        while kmax-kmin > 1:
            kmid = kmax+kmin>>1
            if f(kmid) <= kmid:
                kmax = kmid
            else:
                kmin = kmid
        return kmax
    def f(x):
        c, l = n+x, 0
        j = isqrt(x)
        while j>1:
            k2 = integer_nthroot(x//j**2, 3)[0]+1
            w = squarefreepi(k2-1)
            c -= j*(w-l)
            l, j = w, isqrt(x//k2**3)
        c -= squarefreepi(integer_nthroot(x, 3)[0])-l
        return c
    return divisor_sigma(bisection(f, n, n)) # Chai Wah Wu, Sep 10 2024

Formula

From Amiram Eldar, May 12 2023: (Start)

Sum_{A001694(k) < x} a(k) = c * x^(3/2) + O(x^(23/18 + eps)), where c = A362984 * A090699 / 3 = 1.5572721108... (Jakimczuk and Lalín, 2022). [corrected Sep 21 2024]

Sum_{k=1..n} a(k) ~ c * n^3, where c = A362984 / (3 * A090699^2) = 0.151716514097... . (End)

Extensions

a(1)=1 prepended by and more terms from Giovanni Resta, Feb 06 2018