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

Additionally, a similar identity:

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

`\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: