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