OFFSET
1,4
COMMENTS
The terms of this sequence can be calculated efficiently using the 1st formula. The value of the function f(n, k) is equal to the number of integers i from 1 to n such that the greatest divisor of k that is a unitary divisor of n is 1 if k = 1, or k-free if k >= 2 (k-free numbers are numbers that are not divisible by a k-th power other than 1). E.g., f(n, 1) = A047994(n), f(n, 2) = A384048(n), and f(n, 3) = A384049(n).
LINKS
Amiram Eldar, Table of n, a(n) for n = 1..10000
FORMULA
a(n) = Sum_{k=1..A051903(n)} (n - f(n, k)) = A051903(n) * n - Sum_{k=1..A051903(n)} f(n, k), where f(n, k) is multiplicative for a given k, with f(p^e, k) = p^e - 1 if e >= k and f(p^e, k) = p^e if e < k.
a(n) = 1 if and only if n is prime.
a(n) >= 2 if and only if n is composite.
a(n) >= 2*n - A047994(n) - A384048(n) with equality if and only if n is cubefree that is not squarefree (i.e., n in A067259, or equivalently, A051903(n) = 2).
a(n) < 5*n/4 and lim sun_{n->oo} a(n)/n = 5/4.
Sum_{k=1..n} a(k) ~ c * n^2 / 2, where c = Sum{k>=1} (1 - Product_{p prime} (1 - 1/(p^(2*k-1)*(p+1)))) = 0.36292303251495264373... .
EXAMPLE
MATHEMATICA
f[p_, e_, k_] := p^e - If[e < k, 0, 1]; a[n_] := Module[{fct = FactorInteger[n], emax, s}, emax = Max[fct[[;; , 2]]]; s = emax * n; Do[s -= Times @@ (f[#1, #2, k] & @@@ fct), {k, 1, emax}]; s]; a[1] = 0; Array[a, 100]
PROG
(PARI) a(n) = if(n == 1, 0, my(f = factor(n), p = f[, 1], e = f[, 2], emax = vecmax(e), s = emax*n); for(k = 1, emax, s -= prod(i = 1, #p, p[i]^e[i] - if(e[i] < k, 0, 1))); s);
CROSSREFS
KEYWORD
nonn
AUTHOR
Amiram Eldar, Jun 06 2025
STATUS
approved
