### Curiosities in number theory

In this post I would like to present some interesting exercises in number theory, along with some curious formulas and identities for some number-theoretic functions.

The following list includes the notations used in this post:

### # Generalized Pillai arithmetical function

Given an arithmetical function f(x), we want to compute the following sum:

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

For example, if f(x) = 1 if x=1, and 0 otherwise, then the problem reduces to the Euler's totient function: a(n) = φ(n).

This is interesting, because if the prime factorization of n can be computed, then we can efficiently compute φ(n) using the following formula:

φ(n) = n * \prod_{p|n} (1 - 1/p)

where the product is taken over the distinct prime factors of n.

Now consider that f(x) returns 1 only when x is a prime number (A117494), or only when x is a perfect cube (A078429), etc... Can we have an efficient algorithm to compute a(n) for arbitrary large values of n, assuming that we can factorize n?

#### ## Solution

There exists an elegant elementary formula which uses the φ(n) function and the divisors of n to efficiently compute a(n):

a(n) = \sum_{d|n} f(d) * φ(n/d),

where the sum is taken over the divisors d of n, which can also be stated as:

\sum_{k=1}^n f(gcd(n, k)) * gcd(n, k)^m = \sum_{d|n} f(d) * d^m * φ(n/d)

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

### # Sum of denominators of k^n / n

Given the following function (A279912):

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

can it be computed efficiently, assuming that we know the prime factorization of n?

#### ## Solution

By observing that we can express a(n) as:

a(n) = \sum_{k=1}^n gcd(n, k) * φ(n / gcd(n, k))

where gcd(n, k) is greater than 1 only when k shares a divisor with n, gives us the following simple formula for a prime number p:

a(p) = (p-1)^2 + p

using the fact that φ(p) = p-1 for prime p.

Also, it turns out that a(n) is multiplicative, so we can take the product over the primes of any square-free integer n, using our little formula:

a(n) = \prod_{p|n} ((p-1)^2 + p)

More generally, by using the identity of Euler's totient function for prime powers:

φ(p^k) = p^(k-1) * (p-1)

we can express a(n) as a product over the prime powers of n, for any n >= 1:

a(n) = \prod_{p^k | n} p^(k-1) * ((p-1) * p^k + 1)

Alternatively, we can have a simpler general formula (at the cost of being slightly less efficient to compute), as a sum over the divisors of n:

a(n) = \sum_{d|n} d * (φ(n/d))^2

where φ(k) is the Euler's totient function.

### # Sum of denominators of n^k / k

In the same spirit as the above problem, let (A279911):

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

How fast can we compute a(n)? What about when n is a prime power?

#### ## Solution

This is a slightly harder problem, and unlike the previous function, this function is not multiplicative.

However, we can still say something interesting about it:

a(p^k) = 1/2 (2 + \sum_{j=1}^(2k) (-1)^j * p^j) = (p^(2k+1) + p + 2) / (2(p+1))

for prime powers p^k, where p is prime and k >= 1.

In general, for any n, we have the following formula:

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

where f(n,k) is the largest divisor d of k for which gcd(d, n) = 1, with the optimization f(n,k) = k if gcd(n, k) = 1.

By rewriting the formula in a slightly different way, we have:

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

where m = denom(\frac{n^n}{n!}) =  A095996(n).

Yet another formula, more efficient than the previous ones, is given as:

a(n) = \sum_{k=1}^n gcd({k / p^(v_p(k))  })

where v_p(k) is the p-adic valuation function, which gives the maximal exponent e of p such that p^e divides k, where p are the unique primes dividing n.

### # Number of partitions of n into 2 parts which are not relatively prime

Given a positive integer n, find the number of ways it can be represented as a sum of two integers a + b = n, satisfying the conditions: 2 <= a <= b and gcd(a, n) > 1 and gcd(b, n) > 1. (A082023)

For example, if we let n=24, there are 8 different ways to represent it:

2 + 22 = 24
3 + 21 = 24
4 + 20 = 24
6 + 18 = 24
8 + 16 = 24
9 + 15 = 24
10 + 14 = 24
12 + 12 = 24

Therefore a(24) = 8.

Can we compute a(n) for arbitrary large n?

#### ## Solution

There exists an efficient formula to count the number of such partitions, if we know the prime factorization of n, based on the cototient function:

a(n) = floor((n - φ(n))/2)

### # Factorial less than 10^(10^n)

Find the unique integer k such that k! < 10^(10^n) <= (k+1)!, given a non-negative integer n. (A119906)

#### ## Solution

By taking the logarithm of 10^(10^n) and of Stirling's approximation for n!, we end up with:

k*(log(k)-1) + 1/2 log(2 \pi k) < log(10)*10^n <= (k+1)*(log(k+1)-1) + 1/2 log(2 \pi (k+1))

Now we can use the binary search algorithm to efficiently find a solution in integers to a(n) = k, where the first few solutions are:

a(0) = 3
a(1) = 13
a(2) = 69
a(3) = 449
a(4) = 3248
a(5) = 25205
a(6) = 205022
a(7) = 1723507
a(8) = 14842906
a(9) = 130202808

Notice that a(n) is also the upper limit in computing the e constant to 10^n decimal places, using the following formula:

\sum_{j=0}^(\infty) 1/(j!) = e

For example, if we want to compute the first 10^4 digits of e, we need to sum the terms up to j = a(4) = 3248:

\sum_{j=0}^(a(4)) 1/(j!) = e correct to the first 10^4 digits

### # Summation modulo a Carmichael number

Let's consider the following sum:

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

Can we compute S(n) efficiently for Carmichael numbers, assuming that we can factorize n?

#### ## Solution

The answer is yes, and it follows from the properties of a Carmichael number: a^(n-1) ≡ 1 mod n for every a coprime to n.

By using the Euler totient function φ(n), we can efficiently count the number of values 1 <= a <= n such that gcd(a, n) = 1.

This allows us to create the following efficient formula:

S(n) = (\sum_{d|n} d^(n-1) φ(n/d)) mod n

where n is a Carmichael number.

More generally:

\sum_{k=1}^n k^(φ(n)) ≡ \sum_{k=1}^n gcd(n,k)^(φ(n)) ≡ \sum_{d|n} φ(n/d) d^(φ(n)) mod n

which is based on the identity:

k^(φ(n)) ≡ gcd(n,k)^(φ(n)) mod n

and Euler's theorem, which states that:

a^(φ(n)) ≡ 1 mod n, for gcd(a,n) = 1.

This can be generalized to any arithmetical function f(x), as:

\sum_{k=1}^n f(gcd(n,k)) = \sum_{d|n} f(d) φ(n/d)

### # Curious identities for the sum of divisors of n

Known to G.H. Hardy and E.M. Wright, the following infinite series uses the Ramanujan's sum function c_q(n) to compute the sum of divisors of n, denoted as σ(n):

σ(n) = pi^2 n/6 \sum_{q=1}^(\infty) (c_q(n))/q^2

A closely related formula, also using Ramanujan's sum function, computes σ(n) in a finite number of steps:

σ(n) = \sum_{q=1}^n c_q(n) * floor(n/q)

where floor(x) is the floor function.

A few more interesting identities include:

σ(n) = \sum_{k=1}^n \sum_{j=1}^k cos((2 \pi n j)/k)

σ(n) = \sum_{k=1}^n gcd(n, k) / (φ(n / gcd(n, k)))

Furthermore, the last formula can be generalized for any e ∈ ℂ:

σ_e(n) = \sum_{k=1}^n gcd(n, k)^e / (φ(n / gcd(n, k)))

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

Some useful 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 summation of σ_m(k), with 1 <= k <= n and m >= 0?

#### ## Solution

First, let's take a look at Faulhaber's formula, which is used to 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.

By noticing that:

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

let F_n(x) be the Faulhaber's formula:

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

where B_n(x) are the Bernoulli polynomials. This gives us the following general formula, for any fixed positive integer m:

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

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

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

We can reduce the number of summation terms 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)

If we continue this process up to sqrt(n), we split the sum into two sums, as illustrated in the formula bellow, which takes only O(sqrt(n)) steps to compute:

\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) is Faulhaber's formula defined above, for any integer m >= 0.

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

Given the following function (A143128):

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

is there an efficient algorithm to compute a(n) for large integers n?

#### ## Solution

In the simple case, with m=1 and j=1, the following simple formula can be computed in O(n) arithmetical operations:

a(n) = 1/2 \sum_{k=1}^n k^2 * floor(n/k) * floor(1 + 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^(e+1) * floor(n/k) * floor(1 + 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) * (B_(m+1)(floor(1 + n/k)) - B_(m+1)(1)) / (m+1)

where B_n(x) are the Bernoulli polynomials.

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) * (B_(m+1)(floor(1 + n/k)) - B_(m+1)(1)) / (m+1)

where B_n(x) are the Bernoulli polynomials.

This formula can be computed in O(sqrt(n)) steps.

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

where B_n(x) are the Bernoulli polynomials, then we have:

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

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

An asymptotic approximation for 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.

### # 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 the 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):

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

Can we compute a(n) efficiently?

#### ## Solution

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

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

By flattening the double-summation, we get:

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

There exists an algorithm  to compute a(n) in O(sqrt(n)) steps, but it's a bit too messy, so I will not include it here.

Also, it's very well 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, unfortunately, this would result in a sub-optimal algorithm, since it requires an efficient algorithm for computing the partial-sums of the Euler totient function (A002088), which in turn requires a very fast algorithm for computing the Mertens function.

The complexity of this approach, however, is better than linear time, but worse than O(sqrt(n)) (see: Perl program).

### # Curious identities for the Jordan's totient function J_2(n)

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

J_k(n) = n^k * \prod_{p|n} (1 - 1/p^k)

where the product is taken over the distinct prime factors of n, with the beautiful property:

\sum_{d|n} J_k(d) = n^k

The Jordan's totient function 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 k=1, the function is equivalent with Euler's totient function:

J_1(n) = φ(n)

For k=2, we have the following interesting identities (A007434):

J_2(n) = \sum_{d|n} d^2 * μ(n/d)

J_2(n) = \sum_{k=1}^n gcd(n, k) * φ(gcd(n, k))

J_2(n) = \sum_{k=1}^n gcd(n, k)^2 * cos((2 \pi k)/n)

The first formula uses the Möbius inversion formula, which is defined as:

g(n) = \sum_{d|n} f(d)

f(n) = \sum_{d|n} g(d) * μ(n/d)

In our case, since sum_{d|n} J_2(d) = n^2, we have g(d) = d^2.

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

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 a beautiful formula for computing this partial sum, defined in terms of the Möbius function and Bernoulli polynomials:

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

where μ(k) is the Möbius function and B_n(x) are the Bernoulli 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)

Indeed, 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.

Next, we define the Faulhaber's formula in terms of Bernoulli polynomials, as:

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

then we have:

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

Alternatively, if we define:

M(a, b) = \sum_{k=a}^b μ(k)

then we can rewrite the formula as:

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

An asymptotic approximation is given by the following formula:

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

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

Additionally, the following symmetric identity is satisfied for any 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))

### # Partial sums of the Dedekind psi function

The generalized Dedekind psi function is defined as:

ψ_m(n) = (J_(2 m)(n)) / (J_m(n))

where J_m(n) is the Jordan's totient 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 |μ(k)| * F_m(floor(n/k))

where |μ(k)| is the absolute value of the Möbius function, defined as:

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

with F_n(x) being the Faulhaber's formula, defined in terms of Bernoulli polynomials as:

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

Next, we define S(n) as the sum over the absolute value of the Möbius function |μ(k)| for 1 <= k <= n, which is equivalent with the number of square-free integers between 1 and n. Luckily, 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 more efficient 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))

An asymptotic approximation for 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 any 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 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) is Faulhaber's formula, defined in terms of Bernoulli polynomials as:

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

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) is Faulhaber's formula and P(s) is defined as:

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

with p running 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 following:

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

where for m=0, we have:

ω_0(n) = ω(n)
Ω_0(n) = Ω(n)

While for m=1, we have:

ω_1(n) = A069359(n)
Ω_1(n) = A003415(n)  (the arithmetic derivative of n)

The 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) is Faulhaber's formula.

### # Solutions in Sidef

Bellow I included several Sidef scripts with solutions to the presented exercises:

### # Solutions in Perl

Additionally, some efficient implementations in the Perl programming language: