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:

$φ(n)$ is the Euler's totient function
$J_k(n)$ is the Jordan's totient function
$ψ_k(n)$ is the Dedekind psi function
$σ_k(n)$ is the Divisor (sigma) function
$c_q(n)$ is the Ramanujan's sum function
$F_n(x)$ is the Faulhaber formula
$μ(n)$ is the Möbius function
$ω(n)$ is the Prime omega function
$\pi(n)$ is the Prime-counting function

# 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$ .

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:

Mathematics and computer science