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,2

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).

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)

Example

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

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 *)

Prog

(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

Crossrefs

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

Sequence in context: A290629 A234861 A076401 * A204740 A244536 A237054

Adjacent sequences: A180111 A180112 A180113 * A180115 A180116 A180117

Keyword

easy,nonn

Author

Walter Kehowski, Aug 10 2010

Extensions

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

Status

approved