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

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

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

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

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

`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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.

`ω(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.

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.

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

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

`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

`\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:- Constant e to n decimal places
- Largest factorial less than 10^(10^n)
- Conditional Euler's totient function
- Generalized prime omega function
- Generalized prime big omega function
- Partial sums of the Dedekind psi function
- Partial sums of the Euler totient function
- Partial sums of the Jordan totient function
- Partial sums of the prime omega function
- Partial sums of the sigma function times k^m

### # Solutions in Perl

Additionally, some efficient implementations in the Perl programming language:

- Partial sums of the sigma function
- Partial sums of the sigma function times k
- Partial sums of the sigma function times k^m
- Partial sums of the Euler totient function
- Partial sums of the Jordan totient function
- Partial sums of the Dedekind psi function
- Partial sums of the prime omega_function
- Partial sums of the gcd-sum function

## Comments

## Post a Comment