### Partial sums of arithmetical functions

In this post I would like to present some interesting generalized formulas for computing the partial sums of some arithmetical functions.

We provide sub-linear algorithms for computing the partial sums of the following functions:

### # Definitions

Throughout this post, F_n(x) represents the Faulhaber polynomials, which are used to efficiently sum the m-th powers of the first n positive integers:

\sum_{k=1}^n k^m = (B_(m+1)(n+1) - B_(m+1)(1)) / (m+1)

where B_n(x) are the Bernoulli polynomials, for m >= 0.

The Faulhaber polynomials F_n(x) are defined as:

F_n(x) = (B_(n+1)(x+1) - B_(n+1)(1)) / (n+1)

where B_n(x) are the Bernoulli polynomials, and B_n(1) are the Bernoulli numbers B_n with B_1 = 1/2.

A few special cases of F_n(x) include:

F_0(x) = x

F_1(x) = 1/2 x (x+1)

F_2(x) = 1/6 x (x + 1) (2 x + 1)

An asymptotic approximation for F_n(x) is given by the following formula:

F_n(x)  ~  x^(n+1)/(n+1)

### # General formula for partial sums

For any arithmetical function f(k), we have the following identity:

\sum_{k=1}^n \sum_{d|k} f(d) * (k/d)^m = \sum_{k=1}^n \sum_{d|k} f(k/d) * d^m = \sum_{k=1}^n f(k) * F_m(floor(n/k))

where floor(x) is the floor function and F_n(x) are Faulhaber's polynomials.

We are interested in computing the following partial sum in sub-linear time complexity:

a(n,m) = \sum_{k=1}^n f(k) * F_m(floor(n/k))

If the partial sums of f(k) can be computed efficiently, then we can create a sub-linear algorithm for computing a(n,m), for any integer m >= 0.

Let:

R(n) = \sum_{k=1}^n f(k)

which gives us the following sub-linear formula, assuming that R(n) can be computed efficiently:

a(n,m) = \sum_{k=1}^(floor(sqrt(n))) (R(floor(n/k)) - R(floor(n / (k+1)))) * F_m(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * F_m(floor(n/k))

More generally, for any two functions f(k) and g(k), we have:

\sum_{k=1}^n \sum_{d|k} f(d) * g(k/d) = \sum_{k=1}^n f(k) * \sum_{j=1}^(floor(n/k)) g(j)

If partial sums of f(k) and g(k) can be computed efficiently, then we can also create a sub-linear formula.

Let:

R(n) = \sum_{k=1}^n f(k)

S(n) = \sum_{k=1}^n g(k)

which gives us:

\sum_{k=1}^n \sum_{d|k} f(d) * g(k/d) = \sum_{k=1}^(floor(sqrt(n))) (R(floor(n/k)) - R(floor(n / (k+1)))) * S(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * S(floor(n/k))

Alternatively, using the Dirichlet hyperbola method, we can rewrite the formula as:

\sum_{k=1}^n \sum_{d|k} f(d) * g(k/d) = \sum_{k=1}^floor(sqrt(n)) (f(k) * S(floor(n/k)) + g(k) * R(floor(n/k))) - R(floor(sqrt(n)))*S(floor(sqrt(n)))

### # General recursive formula for partial sums

For a given function f(n), let R(n) and D(n) be defined as:

R(n) = \sum_{k=1}^n f(k)

D(n) = \sum_{k=1}^n \sum_{d|k} f(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) f(j)

By noticing that we can represent D(n) in terms of R(n):

D(n) = \sum_{k=1}^n R(floor(n/k))

gives us the following recursive formula for R(n), with R(1) = f(1):

R(n) = D(n) - \sum_{k=2}^n R(floor(n/k))

which can be computed in sub-linear time if we can compute D(n) efficiently:

R(n) = D(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) R(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) R(k) * (floor(n/k) - floor(n/(k+1)))

To achieve the best sub-linear time, it's required to precompute the values of R(k) for all small values of k, up to about k = O(n^(2/3)), and performing a simple O(1) lookup when k is bellow this bound.

### # Partial sums of the Möbius function μ(k)

The Mertens function M(n) is defined as the partial sums of the Möbius function μ(k):

M(n) = \sum_{k=1}^n μ(k)

We're interested in computing this function in sub-linear time.

Based on the identity:

sum_{k=1}^n \sum_{d|k} μ(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) μ(j) = \sum_{k=1}^n M(floor(n/k)) = 1

we find the following recursive formula for M(n):

M(n) = 1 - \sum_{k=2}^n M(floor(n/k))

which can be computed in sub-linear time as:

M(n) = 1 - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) M(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) M(k) * (floor(n/k) - floor(n/(k+1)))

This requires precomputing the values of M(k) up to about k = O(n^(2/3)) (preferably using a sieve for computing the values of the Möbius function bellow this bound), allowing us to perform a simple O(1) lookup for small k. The use of memoization is also recommended for reducing the time even further. A decent implementation of this algorithm should be able to compute M(10^10) = -33722 in about 10 seconds (or less).

### # Partial sums of the Liouville function λ(k)

The Liouville function λ(n) is defined as:

λ(n) = (-1)^(\Omega(n))

where \Omega(n) is the prime big omega function (the number of prime factors of n counted with multiplicity).

We're interested in computing the following partial sum (A002819):

L(n) = \sum_{k=1}^n λ(k)

which can also be defined in terms of the Mertens function M(x):

L(n) = \sum_{k=1}^floor(sqrt(n)) M(floor(n / k^2))

However, we can do a little better than this, using the following identity:

sum_{k=1}^n \sum_{d|k} λ(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) λ(j) = \sum_{k=1}^n L(floor(n/k)) = floor(sqrt(n))

which allows us to create the following sub-linear recursive formula:

L(n) = floor(sqrt(n)) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) L(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) L(k) * (floor(n/k) - floor(n/(k+1)))

Same as in computing the Mertens function, in order to achieve the best sub-linear time, it's recommended to create a lookup table of L(k) for all small values of k, up to about k = O(n^(2/3)), and perform a simple O(1) operation when k exists in the lookup table.

### # Partial sums of σ_m(k)

The sigma function σ_m(n) is defined as:

σ_m(n) = \sum_{d|n} d^m

where d runs over all the positive divisors of n.

Some formulas for computing the partial sums of the sigma function, include:

\sum_{k=1}^n σ_0(k) = \sum_{k=1}^n floor(n/k)

\sum_{k=1}^n σ_1(k) = 1/2 \sum_{k=1}^n floor(n/k) * floor(1 + n/k)

Can we have a general formula for any partial sum of σ_m(k), with 1 <= k <= n and m >= 0?

Based on the following identity for any arithmetic function f(x):

\sum_{k=1}^n \sum_{d|k} f(d) = \sum_{k=1}^n f(k) * floor(n/k) = \sum_{k=1}^n \sum_{j=1}^(floor(n/k)) f(j)

we have:

\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n \sum_{j=1}^(floor(n/k)) j^m

allowing us to use Faulhaber's formula, which gives us the following general formula, for any integer m>=0:

\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n F_m(floor(n/k))

We can simplify this formula by taking advantage of the following symmetric identity for any fixed integers j>=0 and m>=0:

\sum_{k=1}^n σ_j(k) * F_m(floor(n/k)) = \sum_{k=1}^n σ_m(k) * F_j(floor(n/k))

implying:

\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n k^m * \sum_{j=1}^(floor(n/k)) 1

giving us the simpler formula, for any m ∈ ℂ :

\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n k^m * floor(n/k)

One key point here is noticing the fact that the number of summation terms can be reduced in half as:

\sum_{k=1}^n σ_m(k) = F_m(n) - F_m(floor(n/2)) + \sum_{k=1}^(floor(n/2)) k^m * floor(n/k)

Continuing this process up to sqrt(n), the sum gets split into two sums, where the number of terms in each sum is O(sqrt(n)):

\sum_{k=1}^n σ_m(k) = \sum_{k=1}^(floor(sqrt(n))) k * (F_m(floor(n/k)) - F_m(floor(n/(k+1)))) + \sum_{k=1}^(floor(n/(1+floor(sqrt(n))))) k^m * floor(n/k)

where F_n(x) are Faulhaber's polynomials, for any integer m >= 0.

### # Partial sums of k^m * σ_j(k)

Given the following function:

a(n,m,j) = \sum_{k=1}^n k^m * σ_j(k)

is there an efficient algorithm to compute a(n,m,j) for relatively large integers n?

In the simple case, with m=1 and j=1 (A143128), we have:

\sum_{k=1}^n k * σ(k) = 1/2 \sum_{k=1}^n k^2 * floor(n/k) * floor(1 + n/k) = \sum_{k=1}^n k^2 * F_1(floor(n/k))

In general, for any j ∈ ℂ, we have the following identity:

\sum_{k=1}^n k * σ_j(k) = 1/2 \sum_{k=1}^n k^(j+1) * floor(n/k) * floor(1 + n/k) = \sum_{k=1}^n k^(j+1) * F_1(floor(n/k))

A nice generalization for partial sums of k^m * σ(k) (see A319086), with m >= 0, is:

\sum_{k=1}^n k^m * σ(k) = \sum_{k=1}^n k^(m+1) * F_m(floor(n/k))

Putting this two formulas together, gives us the following full generalization for k^m * σ_j(k), with m >= 0 and j ∈ ℂ:

\sum_{k=1}^n k^m * σ_j(k) = \sum_{k=1}^n k^(m+j) * F_m(floor(n/k))

This formula can be computed in O(sqrt(n)) steps, as:

\sum_{k=1}^n k^m * σ_j(k) = \sum_{k=1}^(floor(sqrt(n))) F_m(k) * (F_(m+j)(floor(n/k)) - F_(m+j)(floor(n/(k+1)))) + \sum_{k=1}^(floor(n / (floor(sqrt(n))+1))) k^(m+j) * F_m(floor(n/k))

which is also equivalent with:

\sum_{k=1}^n k^m * σ_j(k) = \sum_{k=1}^(floor(sqrt(n))) F_(m+j)(k) * (F_m(floor(n/k)) - F_m(floor(n/(k+1)))) + \sum_{k=1}^(floor(n / (floor(sqrt(n))+1))) k^m * F_(m+j)(floor(n/k))

for any fixed integers m >= 0 and j >= 0, where F_n(x) are Faulhaber's polynomials.

Proof:
\sum_{k=1}^n k^m * σ_j(k) = \sum_{k=1}^n k^m * \sum_{d|k} d^j = \sum_{k=1}^n k^m * \sum_{b=1}^(floor(n/k)) b^(m+j) = \sum_{k=1}^n k^(m+j) * \sum_{b=1}^(floor(n/k)) b^m

Since we have an inner sum of m-th powers, we can use Faulhaber's formula to compute \sum_{b=1}^(floor(n/k)) b^m, concluding the proof:

\sum_{k=1}^n k^m * σ_j(k) = \sum_{k=1}^n k^(m+j) * F_m(floor(n/k))

By rewriting the sum using the Dirichlet hyperbola method, we have:

\sum_{k=1}^n k^m * σ_j(k) = \sum_{k=1}^floor(sqrt(n)) (k^m * F_(m+j)(floor(n/k)) + k^(m+j) * F_m(floor(n/k))) - F_(m+j)(floor(sqrt(n))) * F_m(floor(sqrt(n)))

An asymptotic approximation to this partial sum is given by the following formula:

\sum_{k=1}^n k^m * σ_j(k)  ~  (zeta(j+1) * n^(j+m+1)) / (j+m+1)

where zeta(s) is the Riemann zeta function, for m >= 0 and j >= 1.

### # Sum of the number of divisors of gcd(x,y) with x y <= n

Given the following interesting sum (A268732):

a(n) = \sum_{x y <= n} σ_0(gcd(x,y)) = \sum_{k=1}^n \sum_{d|k} σ_0(gcd(d, k/d))

where σ_0(k) is the number of divisors of k, we're interested in computing a(n) in sub-linear time.

This sum was considered by Adrian W. Dudek in his paper (see remark 1, page 3), in which he stated that "it would be interesting to see an asymptotic formula for this.".

Bellow we present a fast sub-linear formula for computing a(n) up to about n=10^15:

a(n) = \sum_{k=1}^floor(sqrt(n)) (2 * \sum_{j=1}^floor(sqrt(n/k^2)) floor(n/(j k^2)) - floor(sqrt(n/k^2))^2)

Evaluating a(10^k) for k=1..15, we get:

a(10^1)  = 31
a(10^2)  = 629
a(10^3)  = 9823
a(10^4)  = 135568
a(10^5)  = 1732437
a(10^6)  = 21107131
a(10^7)  = 248928748
a(10^8)  = 2867996696
a(10^9)  = 32467409097
a(10^10) = 362549612240
a(10^11) = 4004254692640
a(10^12) = 43830142939380
a(10^13) = 476177421208658
a(10^14) = 5140534231877816
a(10^15) = 55192942833495679

The n-th term can be described asymptotically as:

a(n) = n * zeta(2) * (log(n) + 2 γ - 1 + 2 \frac{zeta'(2)}{zeta(2)}) + O(sqrt(n) log(n))

where γ is the Euler-Mascheroni constant.

Alternate form:

a(n) = n * zeta(2) * (-24 log(A) + 2 γ + 2 log(2 π) + log(n) + 2 γ - 1) + O(sqrt(n) log(n))

where γ is the Euler-Mascheroni constant and A is the Glaisher-Kinkelin constant.

### # Partial sums of the number of unitary divisors

A closely related sum to the previous one, is the partial sum of the number of unitary divisors of k with 1 <= k <= n (A064608):

a(n) = \sum_{k=1}^n 2^(ω(k)) = \sum_{k=1}^n μ(k)^2 * floor(n/k)

where ω(k) is the prime omega function and μ(k) is the Möbius function.

The n-th term can be computed in sub-linear time, using the following formula:

a(n) = \sum_{k=1}^floor(sqrt(n)) μ(k) * (2 * \sum_{j=1}^floor(sqrt(n/k^2)) floor(n/(j k^2)) - floor(sqrt(n/k^2))^2)

where μ(k) is the Möbius function.

An asymptotic formula for a(n) is given as:

a(n) = \frac{n}{zeta(2)} * (log(n) + 2 γ - 1 - 2 \frac{zeta'(2)}{zeta(2)}) + O(sqrt(n) log(n))

where γ is the Euler-Mascheroni constant.

### # Partial sums of the gcd-sum function

The Pillai's arithmetical function (A018804), also known as the gcd-sum function, is defined as:

P(n) = \sum_{k=1}^n gcd(n, k)

A simple and efficient formula for computing this function, using the divisors of n and Euler's totient function, is given as:

P(n) = \sum_{d|n} d * φ(n/d)

Now, let's consider the partial sums of this function (A272718):

G(n) = \sum_{k=1}^n P(k) = \sum_{k=1}^n \sum_{j=1}^k gcd(k, j)

Can we compute G(n) efficiently?

Given that we can compute P(n) efficiently, we can use the following formula to compute G(n):

G(n) = \sum_{k=1}^n \sum_{d|k} d * φ(k/d)

By flattening the double-summation, we get:

G(n) = 1/2 \sum_{k=1}^n φ(k) * floor(n/k) * floor(1+n/k)

It's possible to approach this problem in the same way as we did the for the partial-sums of the sigma function, based on the formula:

\sum_{k=1}^n φ(k) = 1/2 \sum_{k=1}^n μ(k) * floor(n/k) * floor(1 + n/k)

where μ(k) is the Möbius function, but this requires a fast algorithm for computing the Mertens function (A002321).

However, we can compute partial sums of Euler's totient function more efficiently, without using the Mertens function.

Let:

Φ(n) = \sum_{k=1}^n φ(k)

Based on the identity:

\sum_{d|n} φ(d) = n

we have:

\sum_{k=1}^n \sum_{d|k} φ(k) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) φ(j) = \sum_{k=1}^n Φ(floor(n/k)) = \sum_{k=1}^n k = F_1(n)

implying the following recurrence:

Φ(n) = F_1(n) - \sum_{k=2}^n Φ(floor(n/k))

which can be computed in sub-linear time as:

Φ(n) = F_1(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) Φ(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) Φ(k) * (floor(n/k) - floor(n/(k+1)))

This requires precomputing Φ(k) up to about k = O(n^(2/3)) and doing a simple O(1) lookup when k is less than the precomputed upper-bound.

In combination with the Dirichlet hyperbola method, we can compute G(n) in sub-linear time as:

G(n) = \sum_{k=1}^floor(sqrt(n)) (k * Φ(floor(n/k)) + φ(k) * F_1(floor(n/k))) - Φ(floor(sqrt(n)))*F_1(floor(sqrt(n)))

where F_n(x) are the Faulhaber polynomials.

### # Partial sums of Jordan's totient function

Jordan's totient function is a generalization of Euler's totient function and is defined as:

J_m(n) = \sum_{d|n} μ(d) (n/d)^m = n^m * \prod_{p|n} (1 - 1/p^m)

where the product is taken over the distinct prime factors of n, and is multiplicative with:

J_m(p^k) = p^(k m) - p^((k - 1) m)

with the beautiful property:

\sum_{d|n} J_m(d) = n^m

Jordan's totient function J_k(n) counts the number of k-tuples (x_1, x_2, ..., x_k) for which gcd(n, x1, x2, ..., x_k) = 1, with 1 <= x_k <= n.

For m=1, the Jordan totient function is equivalent with Euler's totient function:

J_1(n) = φ(n)

For m >= 1, we want to compute:

\sum_{k=1}^n J_m(k)

where J_m(k) is the Jordan totient function.

There is an elegant formula for computing this partial sum, defined in terms of the Möbius function and Faulhaber's formula:

\sum_{k=1}^n J_m(k) = \sum_{k=1}^n \sum_{d|k} μ(d) (k/d)^m = \sum_{k=1}^n μ(k) * F_m(floor(n/k))

where μ(k) is the Möbius function and F_n(x) are Faulhaber's polynomials.

For m=1, we have A002088.
For m=2, we have A321879.

This partial sum can be computed somewhat efficiently, assuming that we have an efficient method for computing the Mertens function, which is defined as:

M(n) = \sum_{k=1}^n μ(k)

There exists a sub-linear time algorithm for computing the Mertens function M(n), due to Marc Deléglise and Joël Rivat, with a time complexity of O(n^(2/3) * (ln ln n)^(1/3)), described in their paper here.

Assuming that we use this algorithm (or a similar one which has roughly the same complexity), then we can also compute this partial sum in sub-linear time:

\sum_{k=1}^n J_m(k) = \sum_{k=1}^(floor(sqrt(n))) (M(floor(n/k)) - M(floor(n / (k+1)))) * F_m(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) μ(k) * F_m(floor(n/k))

However, a more practical and efficient approach uses the following recursive formula:

Let:

Φ_m(n) = \sum_{k=1}^n J_m(k)

which can be expressed recursively as:

Φ_m(n) = F_m(n) - \sum_{k=2}^n Φ_m(floor(n/k))

and computed in sub-linear time as:

Φ_m(n) = F_m(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) Φ_m(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) Φ_m(k) * (floor(n/k) - floor(n/(k+1)))

where F_n(x) are Faulhaber's polynomials.

This approach requires precomputing Φ_m(k) up to about k = 2 * floor(n^(2/3)) and doing a simple lookup when k is less than the precomputed upper-bound.

An asymptotic approximation for partial sums of the Jordan's totient function, is given by the following formula:

\sum_{k=1}^n J_m(k)  ~  n^(m+1) / ((m+1) * zeta(m+1))

or, alternatively:

\sum_{k=1}^n J_m(k)  ~  \frac{F_m(n)}{zeta(m+1)}

where zeta(s) is the Riemann zeta function, for m >= 1.

Additionally, the following symmetric identity is satisfied for all integers m>=0 and u>=0:

\sum_{k=1}^n J_m(k) * F_u(floor(n/k)) = \sum_{k=1}^n J_u(k) * F_m(floor(n/k))

where F_n(x) are Faulhaber's polynomials.

More generally, we also consider the partial sums of J_m(k) * k^p, for m >= 1 and p >= 0:

R(n) = \sum_{k=1}^n J_m(k) * k^p

for which we have the following beautiful recursive formula:

R(n) = F_(m+p)(n) - \sum_{k=2}^n k^p * R(floor(n/k))

which can be computed in sub-linear time as:

R(n) = F_(m+p)(n) - \sum_{k=2}^floor(n / (1+floor(sqrt(n)))) k^p * R(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) R(k) * (F_p(floor(n/k)) - F_p(floor(n/(k+1))))

with the asymptotic formula:

\sum_{k=1}^n J_m(k) * k^p  ~  \frac{F_(m+p)(n)}{\zeta(m+1)}

where J_m(k) is the Jordan totient function, zeta(s) is the Riemann zeta function and F_n(x) are Faulhaber's polynomials, for m >= 1 and p >= 0.

### # Partial sums of the Dedekind psi function

The generalized Dedekind psi function is defined as:

ψ_m(n) = \sum_{d|n} |μ(d)| (n/d)^m = (J_(2 m)(n)) / (J_m(n))

where μ(k) is the Möbius function and J_m(n) is the Jordan totient function, and is multiplicative with:

ψ_m(p^k) = p^(k m) + p^((k - 1) m)

The Möbius transform of the Dedekind psi function can be expressed in terms of the Jordan totient function J_m(k):

\sum_{d|n} μ(d) ψ_m(n/d) = \sum_{d|n} |μ(d)| J_m(n/d)

The inverse Möbius transform of the Dedekind psi function satisfies the following identity:

\sum_{d|n} ψ_m(d) = \sum_{d|n} 2^(ω(d)) (n/d)^m

which gives us the identity:

ψ_m(n) = \sum_{d|n} 2^(ω(d)) J_m(n/d)

where J_m(k) is the Jordan totient function and ω(n) is the prime omega function.

For m >= 1, we're interested in computing the sum (A173290):

\sum_{k=1}^n ψ_m(k)

We can express this sum in terms of the Möbius function:

\sum_{k=1}^n  ψ_m(k) = \sum_{k=1}^n \sum_{d|k} |μ(d)| (k/d)^m = \sum_{k=1}^n |μ(k)| * F_m(floor(n/k))

with F_n(x) is Faulhaber's formula and |μ(k)| is the absolute value of the Möbius function, defined as:

|μ(k)| = 1  iff k is a square-free integer
|μ(k)| = 0  otherwise

Next we define S(n) as the sum over the absolute value of the Möbius function |μ(k)| for 1 <= k <= n, which represents the number of square-free integers between 1 and n. Importantly, S(n) can be computed in O(sqrt(n)) steps:

S(n) = \sum_{k=1}^n |μ(k)| = \sum_{k=1}^(floor(sqrt(n))) μ(k) * floor(n / k^2)

This allows us to create a sub-linear formula for computing the partial sums of the Dedekind psi function:

\sum_{k=1}^n ψ_m(k) = \sum_{k=1}^(floor(sqrt(n))) (S(floor(n/k)) - S(floor(n / (k+1)))) * F_m(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) |μ(k)| * F_m(floor(n/k))

By using the Dirichlet hyperbola method, we can rewrite the formula as:

\sum_{k=1}^n ψ_m(k) = \sum_{k=1}^floor(sqrt(n)) (|μ(k)| * F_m(floor(n/k)) + k^m * S(floor(n/k))) - F_m(floor(sqrt(n)))*S(floor(sqrt(n)))

where F_n(x) are the Faulhaber polynomials.

Additionally, if we define E_m(n) and D_m(n) as:

E_m(n) = \sum_{k=1}^n ψ_m(k)

D_m(n) = \sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k)  ψ_m(j) = \sum_{k=1}^n E_m(floor(n/k))

then we can create the following recursive formula to express partial sums of the Dedekind psi function ψ_m(n):

E_m(n) = D_m(n) - \sum_{k=2}^n E_m(floor(n/k))

which can be computed in sublinear time if we can compute D_m(n) efficiently:

E_m(n) = D_m(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) E_m(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) E_m(k) * (floor(n/k) - floor(n/(k+1)))

An asymptotic approximation to this partial sum is given by:

\sum_{k=1}^n ψ_m(k)  ~  (n^(m+1) * zeta(m+1)) / ((m+1) * zeta(2*(m+1)))

where zeta(s) is the Riemann zeta function, for m >= 1.

Additionally, the following symmetric identity is satisfied for all integers m>=1 and u>=1:

\sum_{k=1}^n ψ_m(k) * F_u(floor(n/k)) = \sum_{k=1}^n ψ_u(k) * F_m(floor(n/k))

### # Partial sums of the inverse Möbius transform of the Dedekind psi function

The inverse Möbius transform of the generalized Dedekind psi function ψ_m(n), is defined as:

\sum_{d|n} ψ_m(d) = \sum_{d|n} 2^(ω(d)) * (n/d)^m

where ω(n) is the prime omega function.

For any fixed integer m >= 0, we're interested in computing the following partial sum (A306775):

\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n \sum_{d|k} 2^(ω(d)) * (k/d)^m

which can be expressed in terms of the Faulhaber polynomials F_n(x), as:

\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n 2^(ω(k)) * F_m(floor(n/k))

which is also equivalent with:

\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n k^m * \sum_{j=1}^floor(n/k) 2^(ω(j))

By defining R(n) as the partial sums of 2^(ω(k)) (see: A064608):

R(n) = \sum_{k=1}^n 2^(ω(k))

we can create the following sublinear formula, assuming that we can compute R(n) efficiently:

\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^floor(sqrt(n)) (2^(ω(k)) * F_m(floor(n/k)) + k^m * R(floor(n/k))) - F_m(floor(sqrt(n))) * R(floor(sqrt(n)))

An asymptotic approximation to this partial sum is given by the following formula:

\sum_{k=1}^n \sum_{d|k} ψ_m(d)  ~  (n^(m+1) * zeta(m+1)^2) / ((m+1) * zeta(2*(m+1)))

where zeta(s) is the Riemann zeta function, for m >= 1.

### # Partial sums of the prime omega function

The prime omega function ω(n) counts the number of distinct prime factors of n.

By definition, the following identity is implied:

ω(n) = ω(rad(n))

where rad(n) is the radical function.

We're interested in efficiently computing the following partial sum (A013939):

a(n) = \sum_{k=1}^n ω(k)

which is equivalent with:

a(n) = \sum_{p <= n} floor(n/p)

where p runs over the prime numbers.

Using the same trick as in the previous sections, we can compute this partial sum in O(sqrt(n)) steps, as:

a(n) = \sum_{k=1}^(floor(sqrt(n))) k*(\pi(floor(n/k)) - \pi(floor(n/(k+1)))) + \sum_{p}^(floor(n/(1+floor(sqrt(n))))) floor(n/p)

where \pi(x) is the prime-counting function and p runs over the prime numbers.

This formula can be generalized as:

A_m(n) = \sum_{p}^n F_m(floor(n/p))

and computed efficiently as:

A_m(n) = \sum_{k=1}^(floor(sqrt(n))) (\pi(floor(n/k)) - \pi(floor(n / (k+1)))) * F_m(k) + \sum_{p}^(floor(n / (floor(sqrt(n)) + 1))) F_m(floor(n/p))

where F_n(x) are Faulhaber's polynomials.

The partial sums of the prime omega function is a special case of A_m(n) with m=0.

This also allows us to generalize the prime omega function as:

ω_m(n) = n^m * \sum_{p|n} 1/p^m

where p are the distinct prime factors of n.

The partial sums of ω_1(n) (A069359, A322068), is given by:

A_1(n) = \sum_{k=1}^n ω_1(k) = 1/2 \sum_{p <= n} floor(n/p) * floor(1 + n/p)

where p runs over the prime numbers.

Evaluating this function at A_1(10^n), we observe a clear asymptotic behavior:

A_1(10^1)  = 25
A_1(10^2)  = 2298
A_1(10^3)  = 226342
A_1(10^4)  = 22616110
A_1(10^5)  = 2261266482
A_1(10^6)  = 226124236118
A_1(10^7)  = 22612374197143
A_1(10^8)  = 2261237139656553
A_1(10^9)  = 226123710243814636
A_1(10^10) = 22612371006991736766
A_1(10^11) = 2261237100241987653515
A_1(10^12) = 226123710021083492369813
A_1(10^13) = 22612371002056432695022703
A_1(10^14) = 2261237100205367824451036203

which can be described as:

A_1(n)  ~  (n (n+1))/2 * 0.4522474200410654985065... (see: A085548)

For m=2 we have:

A_2(n) = \sum_{k=1}^n ω_2(k) = \sum_{p <= n} F_2(floor(n/p))

Evaluating A_2(10^n) for n=1..13, we get:

A_2(10^1)  = 75
A_2(10^2)  = 59962
A_2(10^3)  = 58403906
A_2(10^4)  = 58270913442
A_2(10^5)  = 58255785988898
A_2(10^6)  = 58254390385024132
A_2(10^7)  = 58254229074894448703
A_2(10^8)  = 58254214780225801032503
A_2(10^9)  = 58254213248247357411667320
A_2(10^10) = 58254213116747777047390609694
A_2(10^11) = 58254213101385832019517484266265
A_2(10^12) = 58254213099991292350208499967189227
A_2(10^13) = 58254213099830361065330294973944269431

with the asymptotic behavior:

A_2(n)  ~  (n (n + 1) (2 n + 1))/6 * 0.1747626392994435364231... (see: A085541)

In general, for m >= 1, the asymptotic behavior of A_m(n) can be described in terms of the prime zeta function, as:

A_m(n)  ~  F_m(n) * P(m+1)

where F_m(n) are Faulhaber's polynomials and P(s) is defined as:

P(s) = \sum_{p} 1/p^s

with p running over the prime numbers.

### # Partial sums of the prime big omega function

The prime big omega function, denoted as Ω(n), is defined to be the number of prime factors of n counted with multiplicity.

We're interested in computing the following partial sum:

\sum_{k=1}^n Ω(k)

which is equivalent with:

\sum_{k=1}^n f(k) * floor(n/k)

where f(k) is defined as:

f(k) = 1 if k is a prime power with exponent >= 1
f(k) = 0 otherwise

This formula can be computed in sub-linear time, because we can use the prime-counting function \pi(n) to efficiently compute partial sums of f(n), which represent the number of prime powers <= n:

V(n) = \sum_{k=1}^n f(k) = \sum_{k=1}^(floor(log_2(n))) \pi(floor(n^(1/k)))

allowing us to create the following sub-linear formula:

\sum_{k=1}^n f(k) * floor(n/k) = \sum_{k=1}^(floor(sqrt(n))) (V(floor(n/k)) - V(floor(n / (k+1)))) * k + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * floor(n/k)

with the full generalization:

G_m(n) = \sum_{k=1}^n Ω_m(k) = \sum_{k=1}^n f(k) * F_m(floor(n/k))

computed efficiently as:

G_m(n) = \sum_{k=1}^(floor(sqrt(n))) (V(floor(n/k)) - V(floor(n / (k+1)))) * F_m(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * F_m(floor(n/k))

where F_n(x) are Faulhaber's polynomials.

Evaluating G_1(10^n), we observe an asymptotic behavior:

G_1(10^1)  = 30
G_1(10^2)  = 2815
G_1(10^3)  = 276337
G_1(10^4)  = 27591490
G_1(10^5)  = 2758525172
G_1(10^6)  = 275847515154
G_1(10^7)  = 27584671195911
G_1(10^8)  = 2758466558498626
G_1(10^9)  = 275846649393437566
G_1(10^10) = 27584664891073330599
G_1(10^11) = 2758466488352698209587

which can be described as:

G_1(n)  ~  (n (n+1))/2 * 0.55169329765699918... (see: A154945)

For m=2 we have:

G_2(10^1)  = 82
G_2(10^2)  = 66799
G_2(10^3)  = 64901405
G_2(10^4)  = 64727468210
G_2(10^5)  = 64708096890744
G_2(10^6)  = 64706281936598588
G_2(10^7)  = 64706077322294843451
G_2(10^8)  = 64706058761567362618628
G_2(10^9)  = 64706056807390376400359474
G_2(10^10) = 64706056632561375736945155965
G_2(10^11) = 64706056612919470606889256184409

with the asymptotic behavior:

G_2(n)  ~  (n (n + 1) (2 n + 1))/6 * 0.19411816983263379... (see: A286229)

In general, for m >= 1, the asymptotic behavior of G_m(n) can be described in terms of the prime zeta function, as:

G_m(n)  ~  F_m(n) * Q(m+1)

where F_m(n) are Faulhaber's polynomials and Q(m) is defined as:

Q(m) = \sum_{k=1}^\infty \sum_{p} 1/p^(m*k) = \sum_{p} 1/(p^m - 1)

where p runs over the prime numbers.

### # Partial sums of Dirichlet convolutions of  (-1)^(ω(n)) and (-1)^(Ω(n))

In this section, we consider the following partial sums:

\sum_{k=1}^n \sum_{d|k} (-1)^(ω(d)) * (k/d)^m = \sum_{k=1}^n (-1)^(ω(k)) * F_m(floor(n/k))

and

\sum_{k=1}^n \sum_{d|k} (-1)^(Ω(d)) * (k/d)^m = \sum_{k=1}^n (-1)^(Ω(k)) * F_m(floor(n/k))

where ω(n) and Ω(n) are the prime omega functions, and F_m(n) are the Faulhaber polynomials.

Several well-known identities include:

(-1)^(ω(n)) = μ(rad(n))

(-1)^(Ω(n)) = λ(n)

where rad(n) is the radical function, μ(n) is the Möbius function and λ(n) is the Liouville function.

For a square-free integer n > 1, we have:

\sum_{d|n} (-1)^(ω(d)) = 0

More generally, if r^k is a perfect power, where r is square-free and k >= 2, then:

\sum_{d|r^k} (-1)^(ω(d)) = μ(r) * (k-1)^(ω(r))

For m=0, we have the easy-to-compute identity:

\sum_{k=1}^n (-1)^(Ω(k)) * floor(n/k) = floor(sqrt(n))

For m >= 1, the partial sums can be described asymptotically as:

 \sum_{k=1}^n (-1)^(ω(k)) * F_m(floor(n/k))  ~  F_m(n) * Q(m+1)

\sum_{k=1}^n (-1)^(Ω(k)) * F_m(floor(n/k))  ~  F_m(n) * \frac{zeta(2 (m+1))}{zeta(m+1)}

where zeta(s) is the Riemann zeta function,  and Q(s) is defined as:

Q(s) = \prod_{p} (1 - 1/(p^s - 1))

where p runs over the prime numbers.

### # Generalization of the prime omega functions

The prime big omega function Ω(n) is defined as the number of prime factors of n counted with multiplicity, while ω(n) is defined as the number of distinct prime factors of n:

ω(n) = \sum_{p^k | n} 1

Ω(n) = \sum_{p^k | n} k

where p^k are the prime power factors of n, where n = \prod_{e=1} p_e^(k_e).

This functions can be generalized as:

ω_m(n) = \sum_{p^k | n} n^m / p^m

Ω_m(n) =  \sum_{p^k|n} \sum_{j=1}^k n^m / p^(j * m)

Where for ω_m(n) we have the following OEIS sequences:

ω_0(n) = A001221(n)
ω_1(n) = A069359(n)
ω_2(n) = A322078(n)

While for Ω_m(n) we have:

Ω_0(n) = A001222(n)
Ω_1(n) = A095112(n)
Ω_2(n) = A322664(n)

The generalized prime omega function, ω_m(n) also satisfies the symmetric identity for m>=0 and u>=0:

\sum_{k=1}^n ω_m(k) * F_u(floor(n/k)) = \sum_{k=1}^n ω_u(k) * F_m(floor(n/k))

where F_n(x) are Faulhaber's polynomials.

For taking partial sums of the generalized prime omega function ω_m(k), we have:

\sum_{k=1}^n ω_m(k) = \sum_{p <= n} F_m(floor(n/p))

where p runs over the prime numbers <= n.

Very similar, for Ω_m(k) we have:

\sum_{k=1}^n Ω_m(k) = \sum_{p^k <= n} F_m(floor(n/p^k))

where p^k runs over the prime powers <= n, with k >= 1.

Additionally, we can also consider the following formula, which has some interesting properties:

a_m(n) = \sum_{p^k | n} (k * n^m) / p^m

where p^k are the prime power factors of n, where n = \prod_{e=1} p_e^(k_e).

Several sequences include:

a_0(n) = A001222(n)
a_1(n) = A003415(n)
a_2(n) = A322079(n)

For m >= 1, partial sums of this function can be described asymptotically as:

\sum_{k=1}^n a_m(k)  ~  F_m(n) * \sum_{p} 1/(p^m * (p-1))  ; (see: A136141A152441)

where F_n(x) are Faulhaber's polynomials and p runs over the prime numbers.

### # OEIS sequences

Let's consider the following general partial sum, where f(k) is an arithmetic function:

R_m(n) = \sum_{k=1}^n f(k) * F_m(floor(n/k))

where F_n(x) are Faulhaber's polynomials.

Bellow we have a table of several interesting OEIS sequences for f(k), which produce R_0(n) and R_1(n):

+------------+------------+------------+
|    f(k)    |   R_0(n)    |   R_1(n)    |
+------------+------------+------------+
| A000012(k) | A006218(n) | A024916(n) |
| A000027(k) | A024916(n) | A143127(n) |
| A010051(k) | A013939(n) | A322068(n) |
| A008966(k) | A064608(n) | A173290(n) |
| A008683(k) | A000012(n) | A002088(n) |
| A000290(k) | A064602(n) | A143128(n) |
| A000010(k) | A000217(n) | A272718(n) |
| A034444(k) | A061503(n) | A306775(n) |
| A007434(k) | A000330(n) |            |
| A069513(k) | A022559(n) |            |
+------------+------------+------------+

### # Implementations in Sidef

This section includes several Sidef scripts, implementing the formulas described in this post:

### # Implementations in Perl

Additionally, some efficient implementations in the Perl programming language: