Partial sums of arithmetical functions
In this post we're going to look at some very interesting generalized formulas for computing the partial sums of some arithmetical functions.
`\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`.
`F_n(x) = (B_(n+1)(x+1) - B_(n+1)(1)) / (n+1)`
where `B_n(x)` are the Bernoulli polynomials, and `B_n(1)` are the Bernoulli numbers `B_n` with `B_1 = 1/2`.
Therefore, for `m >= 0`, we have:
`\sum_{k=1}^n k^m = F_m(n)`
Furthermore, the Faulhaber polynomials can also be used in efficiently computing second partial sums of `m`-th powers (sums of power sums), for `m >= 0`:
`\sum_{k=1}^n \sum_{j=1}^k j^m = (n+1) * F_m(n) - F_(m+1)(n)`
A few special cases of `F_n(x)` include:
`F_0(x) = x`
`F_1(x) = 1/2 x (x+1)`
`F_2(x) = 1/6 x (x + 1) (2 x + 1)`
An asymptotic approximation for `F_n(x)` is given by the following formula:
`F_n(x) ~ x^(n+1)/(n+1)`
`\sum_{k=1}^n \sum_{d|k} f(d) * (k/d)^m = \sum_{k=1}^n \sum_{d|k} f(k/d) * d^m = \sum_{k=1}^n f(k) * F_m(floor(n/k))`
where `floor(x)` is the floor function and `F_n(x)` are Faulhaber's polynomials.
We are interested in computing the following partial sum in sublinear time complexity:
`a(n,m) = \sum_{k=1}^n f(k) * F_m(floor(n/k))`
If the partial sums of `f(k)` can be computed efficiently, then we can create a sub-linear algorithm for computing `a(n,m)`, for any integer `m >= 0`.
Let:
`R(n) = \sum_{k=1}^n f(k)`
which gives us the following sub-linear formula, assuming that `R(n)` can be computed efficiently:
`a(n,m) = \sum_{k=1}^(floor(sqrt(n))) (R(floor(n/k)) - R(floor(n / (k+1)))) * F_m(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * F_m(floor(n/k))`
More generally, for any two functions `f(k)` and `g(k)`, we have:
`\sum_{k=1}^n \sum_{d|k} f(d) * g(k/d) = \sum_{k=1}^n f(k) * \sum_{j=1}^(floor(n/k)) g(j)`
If partial sums of `f(k)` and `g(k)` can be computed efficiently, then we can also create a sub-linear formula.
Let:
`R(n) = \sum_{k=1}^n f(k)`
`S(n) = \sum_{k=1}^n g(k)`
which gives us:
`\sum_{k=1}^n \sum_{d|k} f(d) * g(k/d) = \sum_{k=1}^(floor(sqrt(n))) (R(floor(n/k)) - R(floor(n / (k+1)))) * S(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * S(floor(n/k))`
Alternatively, by using the Dirichlet hyperbola method, we can rewrite the formula as:
`\sum_{k=1}^n \sum_{d|k} f(d) * g(k/d) = \sum_{k=1}^floor(sqrt(n)) (f(k) * S(floor(n/k)) + g(k) * R(floor(n/k))) - R(floor(sqrt(n)))*S(floor(sqrt(n)))`
We're interested in computing the following partial sum (A002088) in sub-linear time:
`R(n) = \sum_{k=1}^n \phi(k)`
`\sum_{d|n} \phi(d) = n`
To achieve the best sub-linear time, it's recommended to create a lookup table of `R(k)` for all small values of `k`, up to about `k = floor(n^(2/3))`, and perform a simple `O(1)` operation when `k` exists in the lookup table.
`σ_m(n) = \sum_{d|n} d^m`
where `d` runs over all the positive divisors of `n`.
Some 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 sum of `σ_m(k)`, with `1 <= k <= n` and `m >= 0`?
Based on the following identity for any arithmetic function `f(x)`:
`\sum_{k=1}^n \sum_{d|k} f(d) = \sum_{k=1}^n f(k) * floor(n/k) = \sum_{k=1}^n \sum_{j=1}^(floor(n/k)) f(j)`
we have:
`\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n \sum_{j=1}^(floor(n/k)) j^m`
allowing us to use Faulhaber's formula, which gives us the following general formula, for any integer `m>=0`:
`\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))`
implying:
`\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n k^m * \sum_{j=1}^(floor(n/k)) 1`
giving us the simpler formula, for any `m ∈ ℂ` :
`\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n k^m * floor(n/k)`
One key point here is noticing the fact that the number of summation terms can be reduced 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)`
Continuing this process up to `sqrt(n)`, the sum gets split into two sums, where the number of terms in each sum is `O(sqrt(n))`:
`\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)` are Faulhaber's polynomials, for any integer `m >= 0`.
`a(n,m,j) = \sum_{k=1}^n k^m * σ_j(k)`
is there an efficient algorithm to compute `a(n,m,j)` for relatively large integers `n`?
In the simple case, with `m=1` and `j=1` (A143128), we have:
`\sum_{k=1}^n k * σ(k) = 1/2 \sum_{k=1}^n k^2 * floor(n/k) * floor(1 + n/k) = \sum_{k=1}^n k^2 * F_1(floor(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^(j+1) * floor(n/k) * floor(1 + n/k) = \sum_{k=1}^n k^(j+1) * F_1(floor(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) * F_m(floor(n/k))`
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) * F_m(floor(n/k))`
This formula can be computed in `O(sqrt(n))` steps, as:
`\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`, where `F_n(x)` are Faulhaber's polynomials.
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))`
By rewriting the sum using the Dirichlet hyperbola method, we have:
`\sum_{k=1}^n k^m * σ_j(k) = \sum_{k=1}^floor(sqrt(n)) (k^m * F_(m+j)(floor(n/k)) + k^(m+j) * F_m(floor(n/k))) - F_(m+j)(floor(sqrt(n))) * F_m(floor(sqrt(n)))`
An asymptotic approximation to 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`.
`a(n) = \sum_{x y <= n} σ_0(gcd(x,y)) = \sum_{k=1}^n \sum_{d|k} σ_0(gcd(d, k/d))`
where `σ_0(k)` is the number of divisors of `k`, we're interested in computing `a(n)` in sub-linear time.
This sum was considered by Adrian W. Dudek in his paper (see remark 1, page 3), in which he stated that "it would be interesting to see an asymptotic formula for this.".
Below we present a fast sub-linear formula for computing `a(n)` up to about `n=10^15`:
`a(n) = \sum_{k=1}^floor(sqrt(n)) (2 * \sum_{j=1}^floor(sqrt(n/k^2)) floor(n/(j k^2)) - floor(sqrt(n/k^2))^2)`
Evaluating `a(10^k)` for `k=1..15`, we get:
`a(10^1) = 31`
`a(10^2) = 629`
`a(10^3) = 9823`
`a(10^4) = 135568`
`a(10^5) = 1732437`
`a(10^6) = 21107131`
`a(10^7) = 248928748`
`a(10^8) = 2867996696`
`a(10^9) = 32467409097`
`a(10^10) = 362549612240`
`a(10^11) = 4004254692640`
`a(10^12) = 43830142939380`
`a(10^13) = 476177421208658`
`a(10^14) = 5140534231877816`
`a(10^15) = 55192942833495679`
The `n`-th term can be described asymptotically as:
`a(n) = n * zeta(2) * (log(n) + 2 γ - 1 + 2 \frac{zeta'(2)}{zeta(2)}) + O(sqrt(n) log(n))`
where `γ` is the Euler-Mascheroni constant.
Alternate form:
`a(n) = n * zeta(2) * (-24 log(A) + 2 γ + 2 log(2 π) + log(n) + 2 γ - 1) + O(sqrt(n) log(n))`
where `γ` is the Euler-Mascheroni constant and `A` is the Glaisher-Kinkelin constant.
`a(n) = \sum_{k=1}^n 2^(ω(k)) = \sum_{k=1}^n μ(k)^2 * floor(n/k)`
where `ω(k)` is the prime omega function and `μ(k)` is the Möbius function.
The `n`-th term can be computed in sub-linear time, using the following formula:
`a(n) = \sum_{k=1}^floor(sqrt(n)) μ(k) * (2 * \sum_{j=1}^floor(sqrt(n/k^2)) floor(n/(j k^2)) - floor(sqrt(n/k^2))^2)`
where `μ(k)` is the Möbius function.
An asymptotic formula for `a(n)` is given as:
`a(n) = \frac{n}{zeta(2)} * (log(n) + 2 γ - 1 - 2 \frac{zeta'(2)}{zeta(2)}) + O(sqrt(n) log(n))`
where `γ` is the Euler-Mascheroni constant.
`P(n) = \sum_{k=1}^n gcd(n, k)`
A simple and efficient formula for computing this function, using the divisors of `n` and 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):
`G(n) = \sum_{k=1}^n P(k) = \sum_{k=1}^n \sum_{j=1}^k gcd(k, j)`
Can we compute `G(n)` efficiently?
Given that we can compute `P(n)` efficiently, we can use the following formula to compute `G(n)`:
`G(n) = \sum_{k=1}^n \sum_{d|k} d * φ(k/d)`
By flattening the double-summation, we get:
`G(n) = 1/2 \sum_{k=1}^n φ(k) * floor(n/k) * floor(1+n/k)`
It's 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 this requires a fast algorithm for computing the Mertens function (A002321).
However, we can compute partial sums of Euler's totient function more efficiently, without using the Mertens function.
Let:
`Φ(n) = \sum_{k=1}^n φ(k)`
Based on the identity:
`\sum_{d|n} φ(d) = n`
we have:
`\sum_{k=1}^n \sum_{d|k} φ(k) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) φ(j) = \sum_{k=1}^n Φ(floor(n/k)) = \sum_{k=1}^n k = F_1(n)`
implying the following recurrence:
`Φ(n) = F_1(n) - \sum_{k=2}^n Φ(floor(n/k))`
which can be computed in sub-linear time as:
`Φ(n) = F_1(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) Φ(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) Φ(k) * (floor(n/k) - floor(n/(k+1)))`
This requires precomputing `Φ(k)` up to about `k = floor(n^(2/3))` and performing a simple `O(1)` lookup when `k` is less than the precomputed bound.
In combination with the Dirichlet hyperbola method, we can compute `G(n)` in sub-linear time as:
`G(n) = \sum_{k=1}^floor(sqrt(n)) (k * Φ(floor(n/k)) + φ(k) * F_1(floor(n/k))) - Φ(floor(sqrt(n)))*F_1(floor(sqrt(n)))`
where `F_n(x)` are the Faulhaber polynomials.
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 an elegant formula for computing this partial sum, defined in terms of the Möbius function and Faulhaber's formula:
`\sum_{k=1}^n J_m(k) = \sum_{k=1}^n \sum_{d|k} μ(d) (k/d)^m = \sum_{k=1}^n μ(k) * F_m(floor(n/k))`
where `μ(k)` is the Möbius function and `F_n(x)` are Faulhaber's 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)`
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.
Assuming that we use this algorithm (or a similar one which has roughly the same complexity), then we can also compute this partial sum in sub-linear time:
`\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))`
However, a more practical and efficient approach uses the following recursive formula:
Let:
`Φ_m(n) = \sum_{k=1}^n J_m(k)`
which can be expressed recursively as:
`Φ_m(n) = F_m(n) - \sum_{k=2}^n Φ_m(floor(n/k))`
and computed in sub-linear time as:
`Φ_m(n) = F_m(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) Φ_m(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) Φ_m(k) * (floor(n/k) - floor(n/(k+1)))`
where `F_n(x)` are Faulhaber's polynomials.
This approach requires precomputing `Φ_m(k)` up to about `k = floor(n^(2/3))` and doing a simple lookup when `k` is less than the precomputed bound.
An asymptotic approximation for partial sums of the Jordan's totient function, is given by the following formula:
`\sum_{k=1}^n J_m(k) ~ n^(m+1) / ((m+1) * zeta(m+1))`
or, alternatively:
`\sum_{k=1}^n J_m(k) ~ \frac{F_m(n)}{zeta(m+1)}`
where `zeta(s)` is the Riemann zeta function, for `m >= 1`.
Additionally, the following symmetric identity is satisfied for all 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))`
where `F_n(x)` are Faulhaber's polynomials.
More generally, we also consider the partial sums of `J_m(k) * k^p`, for `m >= 1` and `p >= 0`:
`R(n) = \sum_{k=1}^n J_m(k) * k^p`
for which we have the following beautiful recursive formula:
`R(n) = F_(m+p)(n) - \sum_{k=2}^n k^p * R(floor(n/k))`
which can be computed in sub-linear time as:
`R(n) = F_(m+p)(n) - \sum_{k=2}^floor(n / (1+floor(sqrt(n)))) k^p * R(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) R(k) * (F_p(floor(n/k)) - F_p(floor(n/(k+1))))`
with the asymptotic formula:
`\sum_{k=1}^n J_m(k) * k^p ~ \frac{F_(m+p)(n)}{\zeta(m+1)}`
where `J_m(k)` is the Jordan totient function, `zeta(s)` is the Riemann zeta function and `F_n(x)` are Faulhaber's polynomials, for `m >= 1` and `p >= 0`.
`ψ_m(n) = \sum_{d|n} |μ(d)| (n/d)^m = (J_(2 m)(n)) / (J_m(n))`
where `μ(k)` is the Möbius function and `J_m(n)` is the Jordan totient function, and is multiplicative with:
`ψ_m(p^k) = p^(k m) + p^((k - 1) m)`
The Möbius transform of the Dedekind psi function can be expressed in terms of the Jordan totient function `J_m(k)`:
`\sum_{d|n} μ(d) ψ_m(n/d) = \sum_{d|n} |μ(d)| J_m(n/d)`
The inverse Möbius transform of the Dedekind psi function satisfies the following identity:
`\sum_{d|n} ψ_m(d) = \sum_{d|n} 2^(ω(d)) (n/d)^m`
which gives us the identity:
`ψ_m(n) = \sum_{d|n} 2^(ω(d)) J_m(n/d)`
where `J_m(k)` is the Jordan totient function and `ω(n)` is the prime omega 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 \sum_{d|k} |μ(d)| (k/d)^m = \sum_{k=1}^n |μ(k)| * F_m(floor(n/k))`
with `F_n(x)` is Faulhaber's formula and `|μ(k)|` is the absolute value of the Möbius function, defined as:
`|μ(k)| = 1` iff `k` is a square-free integer
`|μ(k)| = 0` otherwise
Next we define `S(n)` as the sum over the absolute value of the Möbius function `|μ(k)|` for `1 <= k <= n`, which represents the number of square-free integers between `1` and `n`. Importantly, `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 sub-linear 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))`
By using the Dirichlet hyperbola method, we can rewrite the formula as:
`\sum_{k=1}^n ψ_m(k) = \sum_{k=1}^floor(sqrt(n)) (|μ(k)| * F_m(floor(n/k)) + k^m * S(floor(n/k))) - F_m(floor(sqrt(n)))*S(floor(sqrt(n)))`
where `F_n(x)` are the Faulhaber polynomials.
Additionally, if we define `E_m(n)` and `D_m(n)` as:
`E_m(n) = \sum_{k=1}^n ψ_m(k)`
`D_m(n) = \sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) ψ_m(j) = \sum_{k=1}^n E_m(floor(n/k))`
then we can create the following recursive formula to express partial sums of the Dedekind psi function `ψ_m(n)`:
`E_m(n) = D_m(n) - \sum_{k=2}^n E_m(floor(n/k))`
which can be computed in sublinear time if we can compute `D_m(n)` efficiently:
`E_m(n) = D_m(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) E_m(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) E_m(k) * (floor(n/k) - floor(n/(k+1)))`
An asymptotic approximation to 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 all 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))`
`\sum_{d|n} ψ_m(d) = \sum_{d|n} 2^(ω(d)) * (n/d)^m`
where `ω(n)` is the prime omega function.
For any fixed integer `m >= 0`, we're interested in computing the following partial sum (A306775):
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n \sum_{d|k} 2^(ω(d)) * (k/d)^m`
which can be expressed in terms of the Faulhaber polynomials `F_n(x)`, as:
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n ψ_m(k) * floor(n/k) = \sum_{k=1}^n 2^(ω(k)) * F_m(floor(n/k))`
which is also equivalent with:
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n k^m * \sum_{j=1}^floor(n/k) 2^(ω(j))`
By defining `R(n)` as the partial sums of `2^(ω(k))` (see: A064608):
`R(n) = \sum_{k=1}^n 2^(ω(k))`
we can create the following sublinear formula, assuming that we can compute `R(n)` efficiently:
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^floor(sqrt(n)) (2^(ω(k)) * F_m(floor(n/k)) + k^m * R(floor(n/k))) - F_m(floor(sqrt(n))) * R(floor(sqrt(n)))`
An asymptotic approximation to this partial sum is given by the following formula:
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) ~ (n^(m+1) * zeta(m+1)^2) / ((m+1) * zeta(2*(m+1)))`
where `zeta(s)` is the Riemann zeta function, for `m >= 1`.
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))`
or, alternatively:
`A_m(n) = \sum_{k=1}^floor(sqrt(n)) k^m * \pi(floor(n/k)) + \sum_{p}^floor(sqrt(n)) F_m(floor(n/p)) - F_m(floor(sqrt(n))) * \pi(floor(sqrt(n)))`
where `F_n(x)` are Faulhaber's polynomials, `\pi(x)` is the prime-counting function and `p` runs over the prime numbers.
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)` are Faulhaber's polynomials and `P(s)` is defined as:
`P(s) = \sum_{p} 1/p^s`
with `p` running over the prime numbers.
We're interested in computing the following partial sum:
`\sum_{k=1}^n Ω(k)`
which is equivalent with:
`\sum_{k=1}^n f(k) * floor(n/k)`
where `f(k)` is defined as:
`f(k) = 1` if `k` is a prime power with exponent `>= 1`
`f(k) = 0` otherwise
Let:
`V(n) = \sum_{k=1}^n f(k)`
This formula can be computed in sub-linear time, because we can use the prime-counting function `\pi(n)` to efficiently compute partial sums of `f(n)`, which represent the number of prime powers `<= n`:
`V(n) = \sum_{k=1}^(floor(log_2(n))) \pi(floor(n^(1/k)))`
allowing us to create the following sub-linear formula:
`\sum_{k=1}^n f(k) * floor(n/k) = \sum_{k=1}^(floor(sqrt(n))) (V(floor(n/k)) - V(floor(n / (k+1)))) * k + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * floor(n/k)`
with the full generalization:
`G_m(n) = \sum_{k=1}^n Ω_m(k) = \sum_{k=1}^n f(k) * F_m(floor(n/k))`
computed efficiently as:
`G_m(n) = \sum_{k=1}^floor(sqrt(n)) (V(floor(n/k)) - V(floor(n / (k+1)))) * F_m(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * F_m(floor(n/k))`
or, alternatively:
`G_m(n) = \sum_{k=1}^floor(sqrt(n)) (k^m * V(floor(n/k)) + f(k) * F_m(floor(n/k))) - F_m(floor(sqrt(n))) * V(floor(sqrt(n)))`
where `F_n(x)` are Faulhaber's polynomials.
Evaluating `G_1(10^n)`, we observe an asymptotic behavior:
`G_1(10^1) = 30`
`G_1(10^2) = 2815`
`G_1(10^3) = 276337`
`G_1(10^4) = 27591490`
`G_1(10^5) = 2758525172`
`G_1(10^6) = 275847515154`
`G_1(10^7) = 27584671195911`
`G_1(10^8) = 2758466558498626`
`G_1(10^9) = 275846649393437566`
`G_1(10^10) = 27584664891073330599`
`G_1(10^11) = 2758466488352698209587`
which can be described as:
`G_1(n) ~ (n (n+1))/2 * 0.55169329765699918...` (see: A154945)
For `m=2` we have:
`G_2(10^1) = 82`
`G_2(10^2) = 66799`
`G_2(10^3) = 64901405`
`G_2(10^4) = 64727468210`
`G_2(10^5) = 64708096890744`
`G_2(10^6) = 64706281936598588`
`G_2(10^7) = 64706077322294843451`
`G_2(10^8) = 64706058761567362618628`
`G_2(10^9) = 64706056807390376400359474`
`G_2(10^10) = 64706056632561375736945155965`
`G_2(10^11) = 64706056612919470606889256184409`
with the asymptotic behavior:
`G_2(n) ~ (n (n + 1) (2 n + 1))/6 * 0.19411816983263379...` (see: A286229)
In general, for `m >= 1`, the asymptotic behavior of `G_m(n)` can be described in terms of the prime zeta function, as:
`G_m(n) ~ F_m(n) * Q(m+1)`
where `F_m(n)` are Faulhaber's polynomials and `Q(m)` is defined as:
`Q(m) = \sum_{k=1}^\infty \sum_{p} 1/p^(m*k) = \sum_{p} 1/(p^m - 1)`
where `p` runs over the prime numbers.
`\sum_{k=1}^n \sum_{d|k} (-1)^(ω(d)) * (k/d)^m = \sum_{k=1}^n (-1)^(ω(k)) * F_m(floor(n/k))`
and
`\sum_{k=1}^n \sum_{d|k} (-1)^(Ω(d)) * (k/d)^m = \sum_{k=1}^n (-1)^(Ω(k)) * F_m(floor(n/k))`
where `ω(n)` and `Ω(n)` are the prime omega functions, and `F_m(n)` are the Faulhaber polynomials.
Several well-known identities include:
`(-1)^(ω(n)) = μ(rad(n))`
`(-1)^(Ω(n)) = λ(n)`
where `rad(n)` is the radical function, `μ(n)` is the Möbius function and `λ(n)` is the Liouville function.
For a square-free integer `n > 1`, we have:
`\sum_{d|n} (-1)^(ω(d)) = 0`
More generally, if `r^k` is a perfect power, where `r` is square-free and `k >= 2`, then:
`\sum_{d|r^k} (-1)^(ω(d)) = μ(r) * (k-1)^(ω(r))`
For `m=0`, we have the easy-to-compute identity:
`\sum_{k=1}^n (-1)^(Ω(k)) * floor(n/k) = floor(sqrt(n))`
For `m >= 1`, the partial sums can be described asymptotically as:
` \sum_{k=1}^n (-1)^(ω(k)) * F_m(floor(n/k)) ~ F_m(n) * Q(m+1)`
`\sum_{k=1}^n (-1)^(Ω(k)) * F_m(floor(n/k)) ~ F_m(n) * \frac{zeta(2 (m+1))}{zeta(m+1)}`
where `zeta(s)` is the Riemann zeta function, and `Q(s)` is defined as:
`Q(s) = \prod_{p} (1 - 1/(p^s - 1))`
where `p` runs 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:
`ω_m(n) = \sum_{p^k | n} n^m / p^m`
`Ω_m(n) = \sum_{p^k|n} \sum_{j=1}^k n^m / p^(j * m)`
Where for `ω_m(n)` we have the following OEIS sequences:
`ω_0(n)` = A001221(n)
`ω_1(n)` = A069359(n)
`ω_2(n)` = A322078(n)
While for `Ω_m(n)` we have:
`Ω_0(n)` = A001222(n)
`Ω_1(n)` = A095112(n)
`Ω_2(n)` = A322664(n)
The generalized 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)` are Faulhaber's polynomials.
For taking partial sums of the generalized prime omega function `ω_m(k)`, we have:
`\sum_{k=1}^n ω_m(k) = \sum_{p <= n} F_m(floor(n/p))`
where `p` runs over the prime numbers `<= n`.
Very similar, for `Ω_m(k)` we have:
`\sum_{k=1}^n Ω_m(k) = \sum_{p^k <= n} F_m(floor(n/p^k))`
where `p^k` runs over the prime powers `<= n`, with `k >= 1`.
Additionally, we can also consider the following formula, which has some interesting properties:
`a_m(n) = \sum_{p^k | n} (k * n^m) / p^m`
where `p^k` are the prime power factors of `n`, where `n = \prod_{e=1} p_e^(k_e)`.
Several sequences include:
`a_0(n)` = A001222(n)
`a_1(n)` = A003415(n)
`a_2(n)` = A322079(n)
For `m >= 1`, partial sums of this function can be described asymptotically as:
`\sum_{k=1}^n a_m(k) ~ F_m(n) * \sum_{p} 1/(p^m * (p-1))` ; (see: A136141, A152441)
where `F_n(x)` are Faulhaber's polynomials and `p` runs over the prime numbers.
`R_m(n) = \sum_{k=1}^n f(k) * F_m(floor(n/k))`
where `F_n(x)` are Faulhaber's polynomials.
Below we have a table of several interesting OEIS sequences for `f(k)`, which produce `R_0(n)` and `R_1(n)`:
+------------+------------+------------+
| `f(k)` | `R_0(n)` | `R_1(n)` |
+------------+------------+------------+
| A000012(k) | A006218(n) | A024916(n) |
| A000027(k) | A024916(n) | A143127(n) |
| A010051(k) | A013939(n) | A322068(n) |
| A008966(k) | A064608(n) | A173290(n) |
| A008683(k) | A000012(n) | A002088(n) |
| A000290(k) | A064602(n) | A143128(n) |
| A000010(k) | A000217(n) | A272718(n) |
| A034444(k) | A061503(n) | A306775(n) |
We provide sub-linear algorithms for computing the partial sums of the following functions:
- Möbius function `μ(n)`
- Euler's totient function `φ(n)`
- Jordan's totient function `J_k(n)`
- Dedekind psi function `ψ_k(n)`
- Divisor (sigma) function `σ_k(n)`
- Prime omega function `ω(n)`
- Liouville function `λ(n)`
- Pillai's arithmetical function
# Definitions
Throughout this post, `F_n(x)` represents the Faulhaber polynomials, which are used to efficiently 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`.
The Faulhaber polynomials `F_n(x)` are defined as:
`F_n(x) = (B_(n+1)(x+1) - B_(n+1)(1)) / (n+1)`
where `B_n(x)` are the Bernoulli polynomials, and `B_n(1)` are the Bernoulli numbers `B_n` with `B_1 = 1/2`.
Therefore, for `m >= 0`, we have:
`\sum_{k=1}^n k^m = F_m(n)`
Furthermore, the Faulhaber polynomials can also be used in efficiently computing second partial sums of `m`-th powers (sums of power sums), for `m >= 0`:
`\sum_{k=1}^n \sum_{j=1}^k j^m = (n+1) * F_m(n) - F_(m+1)(n)`
A few special cases of `F_n(x)` include:
`F_0(x) = x`
`F_1(x) = 1/2 x (x+1)`
`F_2(x) = 1/6 x (x + 1) (2 x + 1)`
An asymptotic approximation for `F_n(x)` is given by the following formula:
`F_n(x) ~ x^(n+1)/(n+1)`
# General formula for partial sums
For any arithmetical function `f(k)`, we have the following identity:`\sum_{k=1}^n \sum_{d|k} f(d) * (k/d)^m = \sum_{k=1}^n \sum_{d|k} f(k/d) * d^m = \sum_{k=1}^n f(k) * F_m(floor(n/k))`
where `floor(x)` is the floor function and `F_n(x)` are Faulhaber's polynomials.
We are interested in computing the following partial sum in sublinear time complexity:
`a(n,m) = \sum_{k=1}^n f(k) * F_m(floor(n/k))`
If the partial sums of `f(k)` can be computed efficiently, then we can create a sub-linear algorithm for computing `a(n,m)`, for any integer `m >= 0`.
Let:
`R(n) = \sum_{k=1}^n f(k)`
which gives us the following sub-linear formula, assuming that `R(n)` can be computed efficiently:
`a(n,m) = \sum_{k=1}^(floor(sqrt(n))) (R(floor(n/k)) - R(floor(n / (k+1)))) * F_m(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * F_m(floor(n/k))`
More generally, for any two functions `f(k)` and `g(k)`, we have:
`\sum_{k=1}^n \sum_{d|k} f(d) * g(k/d) = \sum_{k=1}^n f(k) * \sum_{j=1}^(floor(n/k)) g(j)`
If partial sums of `f(k)` and `g(k)` can be computed efficiently, then we can also create a sub-linear formula.
Let:
`R(n) = \sum_{k=1}^n f(k)`
`S(n) = \sum_{k=1}^n g(k)`
which gives us:
`\sum_{k=1}^n \sum_{d|k} f(d) * g(k/d) = \sum_{k=1}^(floor(sqrt(n))) (R(floor(n/k)) - R(floor(n / (k+1)))) * S(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * S(floor(n/k))`
Alternatively, by using the Dirichlet hyperbola method, we can rewrite the formula as:
`\sum_{k=1}^n \sum_{d|k} f(d) * g(k/d) = \sum_{k=1}^floor(sqrt(n)) (f(k) * S(floor(n/k)) + g(k) * R(floor(n/k))) - R(floor(sqrt(n)))*S(floor(sqrt(n)))`
# Möbius inversion formula
When we cannot compute partial sums of `f(k)` efficiently, we can use the following trick to find the inverse Möbius transform of `f` and hope that we can compute partial sums of `g(k)` efficiently:
`g(n) = \sum_{d|n} \mu(d) * f(n/d)`
which gives us the following summation formula:
`\sum_{k=1}^n f(k) = \sum_{k=1}^n g(k) * \floor(n/k)`
`S(n) = \sum_{k=1}^n \sum_{j=1}^k f(j) = \sum_{k=1}^n f(k) * (n-k+1)`
can be represented as:
`S(n) = (n+1) * A(n) - B(n)`
where:
`A(n) = \sum_{k=1}^n f(k)`
`B(n) = \sum_{k=1}^n k * f(k)`
By computing `A(n)` and `B(n)` in sub-linear time, allows us to compute `S(n)` in sub-linear time as well.
`R(n) = \sum_{k=1}^n f(k)`
`D(n) = \sum_{k=1}^n \sum_{d|k} f(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) f(j)`
By noticing that we can represent `D(n)` in terms of `R(n)`:
`D(n) = \sum_{k=1}^n R(floor(n/k))`
gives us the following recursive formula for `R(n)`, with `R(1) = f(1)`:
`R(n) = D(n) - \sum_{k=2}^n R(floor(n/k))`
which can be computed in sub-linear time if we can compute `D(n)` efficiently:
`R(n) = D(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) R(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) R(k) * (floor(n/k) - floor(n/(k+1)))`
To achieve the best sub-linear time, it's required to precompute a lookup table of `R(k)` for all small values of `k`, up to about `k = floor(n^(2/3))`, and perform a fast `O(1)` lookup when `k` is below this bound.
`M(n) = \sum_{k=1}^n μ(k)`
We're interested in computing this function in sub-linear time.
Based on the identity:
`sum_{k=1}^n \sum_{d|k} μ(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) μ(j) = \sum_{k=1}^n M(floor(n/k)) = 1`
we find the following recursive formula for `M(n)`:
`M(n) = 1 - \sum_{k=2}^n M(floor(n/k))`
which can be computed in sub-linear time as:
`M(n) = 1 - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) M(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) M(k) * (floor(n/k) - floor(n/(k+1)))`
This requires precomputing the values of `M(k)` up to about `k = floor(n^(2/3))` (preferably using a sieve for computing the values of the Möbius function below this bound), allowing us to perform a simple `O(1)` lookup for small `k`. The use of memoization is also recommended for reducing the time even further. A decent implementation of this algorithm should be able to compute `M(10^10) = -33722` in about `10` seconds (or less).
`λ(n) = (-1)^(\Omega(n))`
where `\Omega(n)` is the prime big omega function (the number of prime factors of `n` counted with multiplicity).
We're interested in computing the following partial sum (A002819):
`L(n) = \sum_{k=1}^n λ(k)`
which can also be defined in terms of the Mertens function `M(x)`:
`L(n) = \sum_{k=1}^floor(sqrt(n)) M(floor(n / k^2))`
However, we can do a little better than this, using the following identity:
`sum_{k=1}^n \sum_{d|k} λ(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) λ(j) = \sum_{k=1}^n L(floor(n/k)) = floor(sqrt(n))`
which allows us to create the following sub-linear recursive formula:
`L(n) = floor(sqrt(n)) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) L(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) L(k) * (floor(n/k) - floor(n/(k+1)))`
Same as in computing the Mertens function, in order to achieve the best sub-linear time, it's recommended to create a lookup table of `L(k)` for all small values of `k`, up to about `k = floor(n^(2/3))`, and perform a simple `O(1)` operation when `k` exists in the lookup table.
# Double-summation
Any double-summation of the form:`S(n) = \sum_{k=1}^n \sum_{j=1}^k f(j) = \sum_{k=1}^n f(k) * (n-k+1)`
can be represented as:
`S(n) = (n+1) * A(n) - B(n)`
where:
`A(n) = \sum_{k=1}^n f(k)`
`B(n) = \sum_{k=1}^n k * f(k)`
By computing `A(n)` and `B(n)` in sub-linear time, allows us to compute `S(n)` in sub-linear time as well.
# General recursive formula for partial sums
For a given function `f(n)`, let `R(n)` and `D(n)` be defined as:`R(n) = \sum_{k=1}^n f(k)`
`D(n) = \sum_{k=1}^n \sum_{d|k} f(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) f(j)`
By noticing that we can represent `D(n)` in terms of `R(n)`:
`D(n) = \sum_{k=1}^n R(floor(n/k))`
gives us the following recursive formula for `R(n)`, with `R(1) = f(1)`:
`R(n) = D(n) - \sum_{k=2}^n R(floor(n/k))`
which can be computed in sub-linear time if we can compute `D(n)` efficiently:
`R(n) = D(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) R(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) R(k) * (floor(n/k) - floor(n/(k+1)))`
To achieve the best sub-linear time, it's required to precompute a lookup table of `R(k)` for all small values of `k`, up to about `k = floor(n^(2/3))`, and perform a fast `O(1)` lookup when `k` is below this bound.
# Partial sums of the Möbius function `μ(k)`
The Mertens function `M(n)` is defined as the partial sums of the Möbius function `μ(k)`:`M(n) = \sum_{k=1}^n μ(k)`
We're interested in computing this function in sub-linear time.
Based on the identity:
`sum_{k=1}^n \sum_{d|k} μ(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) μ(j) = \sum_{k=1}^n M(floor(n/k)) = 1`
we find the following recursive formula for `M(n)`:
`M(n) = 1 - \sum_{k=2}^n M(floor(n/k))`
which can be computed in sub-linear time as:
`M(n) = 1 - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) M(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) M(k) * (floor(n/k) - floor(n/(k+1)))`
This requires precomputing the values of `M(k)` up to about `k = floor(n^(2/3))` (preferably using a sieve for computing the values of the Möbius function below this bound), allowing us to perform a simple `O(1)` lookup for small `k`. The use of memoization is also recommended for reducing the time even further. A decent implementation of this algorithm should be able to compute `M(10^10) = -33722` in about `10` seconds (or less).
# Partial sums of the Liouville function `λ(k)`
The Liouville function `λ(n)` is defined as:`λ(n) = (-1)^(\Omega(n))`
where `\Omega(n)` is the prime big omega function (the number of prime factors of `n` counted with multiplicity).
We're interested in computing the following partial sum (A002819):
`L(n) = \sum_{k=1}^n λ(k)`
which can also be defined in terms of the Mertens function `M(x)`:
`L(n) = \sum_{k=1}^floor(sqrt(n)) M(floor(n / k^2))`
However, we can do a little better than this, using the following identity:
`sum_{k=1}^n \sum_{d|k} λ(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) λ(j) = \sum_{k=1}^n L(floor(n/k)) = floor(sqrt(n))`
which allows us to create the following sub-linear recursive formula:
`L(n) = floor(sqrt(n)) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) L(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) L(k) * (floor(n/k) - floor(n/(k+1)))`
Same as in computing the Mertens function, in order to achieve the best sub-linear time, it's recommended to create a lookup table of `L(k)` for all small values of `k`, up to about `k = floor(n^(2/3))`, and perform a simple `O(1)` operation when `k` exists in the lookup table.
# Partial sums of the Euler totient function `\phi(k)`
Euler's totient function `\phi(n)` is defined as the number of positive integers `<= n` that are coprime to `n`.
We're interested in computing the following partial sum (A002088) in sub-linear time:
`R(n) = \sum_{k=1}^n \phi(k)`
Based on the following identity:
`\sum_{d|n} \phi(d) = n`
we have:
`\sum_{k=1}^n \sum_{d|k} \phi(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) \phi(j) = \sum_{k=1}^n R(floor(n/k)) = F_1(n)`
from which we can create the following sub-linear recursive formula:
`R(n) = F_1(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) R(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) R(k) * (floor(n/k) - floor(n/(k+1)))`
from which we can create the following sub-linear recursive formula:
`R(n) = F_1(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) R(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) R(k) * (floor(n/k) - floor(n/(k+1)))`
To achieve the best sub-linear time, it's recommended to create a lookup table of `R(k)` for all small values of `k`, up to about `k = floor(n^(2/3))`, and perform a simple `O(1)` operation when `k` exists in the lookup table.
An alternative sub-linear formula, using the Mertens function `M(n)`, is given as:
`R(n) = \sum_{k=1}^floor(sqrt(n)) (k*M(floor(n/k)) + \mu(k)*F_1(floor(n/k))) - M(floor(sqrt(n)))*F_1(floor(sqrt(n)))`
# Partial sums of alternating sum of divisors function
Let `f(p^k) = p^k - p^(k-1) + p^(k-2) - ... +- 1`, and then extend by multiplicativity. This is sequence A206369, for which we have the following simpler identities:
`f(n) = \sum_{d|n} (-1)^{\Omega(n/d)} * d`
`f(n) = \sum_{d^2|n} \phi(n/d^2)`
Now let's consider partial sums of `f(n)`. This is the sequence A370905:
`a(n) = \sum_{k=1}^n f(k) = 1/2 \sum_{k=1}^n (-1)^{\Omega(k)} * floor(n/k) * floor(1 + n/k)`
By using the second identity of `f(n)`, we can compute `a(n)` in sublinear-time by applying the Dirichlet hyperbola method:
`a(n) = \sum_{k=1}^{floor(n^(1/4))} S(floor(n/{k^2})) + \sum_{k=1}^{floor(\sqrt(n))} \phi(k) * floor(\sqrt(floor(n/k))) - S(floor(\sqrt(n))) * floor(n^(1/4))`
where `S(n)` are the partial sums of the Euler totient function:
`S(n) = \sum_{k=1}^n \phi(k)`
...which can be computed in sublinear time, as shown in the previous section.
# Partial sums of `σ_m(k)`
The sigma function `σ_m(n)` is defined as:`σ_m(n) = \sum_{d|n} d^m`
where `d` runs over all the positive divisors of `n`.
Some 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 sum of `σ_m(k)`, with `1 <= k <= n` and `m >= 0`?
Based on the following identity for any arithmetic function `f(x)`:
`\sum_{k=1}^n \sum_{d|k} f(d) = \sum_{k=1}^n f(k) * floor(n/k) = \sum_{k=1}^n \sum_{j=1}^(floor(n/k)) f(j)`
we have:
`\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n \sum_{j=1}^(floor(n/k)) j^m`
allowing us to use Faulhaber's formula, which gives us the following general formula, for any integer `m>=0`:
`\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))`
implying:
`\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n k^m * \sum_{j=1}^(floor(n/k)) 1`
giving us the simpler formula, for any `m ∈ ℂ` :
`\sum_{k=1}^n σ_m(k) = \sum_{k=1}^n k^m * floor(n/k)`
One key point here is noticing the fact that the number of summation terms can be reduced 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)`
Continuing this process up to `sqrt(n)`, the sum gets split into two sums, where the number of terms in each sum is `O(sqrt(n))`:
`\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)` are Faulhaber's polynomials, for any integer `m >= 0`.
# Partial sums of `k^m * σ_j(k)`
Given the following function:`a(n,m,j) = \sum_{k=1}^n k^m * σ_j(k)`
is there an efficient algorithm to compute `a(n,m,j)` for relatively large integers `n`?
In the simple case, with `m=1` and `j=1` (A143128), we have:
`\sum_{k=1}^n k * σ(k) = 1/2 \sum_{k=1}^n k^2 * floor(n/k) * floor(1 + n/k) = \sum_{k=1}^n k^2 * F_1(floor(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^(j+1) * floor(n/k) * floor(1 + n/k) = \sum_{k=1}^n k^(j+1) * F_1(floor(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) * F_m(floor(n/k))`
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) * F_m(floor(n/k))`
This formula can be computed in `O(sqrt(n))` steps, as:
`\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`, where `F_n(x)` are Faulhaber's polynomials.
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))`
By rewriting the sum using the Dirichlet hyperbola method, we have:
`\sum_{k=1}^n k^m * σ_j(k) = \sum_{k=1}^floor(sqrt(n)) (k^m * F_(m+j)(floor(n/k)) + k^(m+j) * F_m(floor(n/k))) - F_(m+j)(floor(sqrt(n))) * F_m(floor(sqrt(n)))`
An asymptotic approximation to 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`.
# Sum of the number of divisors of `gcd(x,y)` with `x y <= n`
Given the following interesting sum (A268732):`a(n) = \sum_{x y <= n} σ_0(gcd(x,y)) = \sum_{k=1}^n \sum_{d|k} σ_0(gcd(d, k/d))`
where `σ_0(k)` is the number of divisors of `k`, we're interested in computing `a(n)` in sub-linear time.
This sum was considered by Adrian W. Dudek in his paper (see remark 1, page 3), in which he stated that "it would be interesting to see an asymptotic formula for this.".
Below we present a fast sub-linear formula for computing `a(n)` up to about `n=10^15`:
`a(n) = \sum_{k=1}^floor(sqrt(n)) (2 * \sum_{j=1}^floor(sqrt(n/k^2)) floor(n/(j k^2)) - floor(sqrt(n/k^2))^2)`
Evaluating `a(10^k)` for `k=1..15`, we get:
`a(10^1) = 31`
`a(10^2) = 629`
`a(10^3) = 9823`
`a(10^4) = 135568`
`a(10^5) = 1732437`
`a(10^6) = 21107131`
`a(10^7) = 248928748`
`a(10^8) = 2867996696`
`a(10^9) = 32467409097`
`a(10^10) = 362549612240`
`a(10^11) = 4004254692640`
`a(10^12) = 43830142939380`
`a(10^13) = 476177421208658`
`a(10^14) = 5140534231877816`
`a(10^15) = 55192942833495679`
The `n`-th term can be described asymptotically as:
`a(n) = n * zeta(2) * (log(n) + 2 γ - 1 + 2 \frac{zeta'(2)}{zeta(2)}) + O(sqrt(n) log(n))`
where `γ` is the Euler-Mascheroni constant.
Alternate form:
`a(n) = n * zeta(2) * (-24 log(A) + 2 γ + 2 log(2 π) + log(n) + 2 γ - 1) + O(sqrt(n) log(n))`
where `γ` is the Euler-Mascheroni constant and `A` is the Glaisher-Kinkelin constant.
# Partial sums of the number of unitary divisors
A closely related sum to the previous one, is the partial sum of the number of unitary divisors of `k` with `1 <= k <= n` (A064608):`a(n) = \sum_{k=1}^n 2^(ω(k)) = \sum_{k=1}^n μ(k)^2 * floor(n/k)`
where `ω(k)` is the prime omega function and `μ(k)` is the Möbius function.
The `n`-th term can be computed in sub-linear time, using the following formula:
`a(n) = \sum_{k=1}^floor(sqrt(n)) μ(k) * (2 * \sum_{j=1}^floor(sqrt(n/k^2)) floor(n/(j k^2)) - floor(sqrt(n/k^2))^2)`
where `μ(k)` is the Möbius function.
An asymptotic formula for `a(n)` is given as:
`a(n) = \frac{n}{zeta(2)} * (log(n) + 2 γ - 1 - 2 \frac{zeta'(2)}{zeta(2)}) + O(sqrt(n) log(n))`
where `γ` is the Euler-Mascheroni constant.
# 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 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):
`G(n) = \sum_{k=1}^n P(k) = \sum_{k=1}^n \sum_{j=1}^k gcd(k, j)`
Can we compute `G(n)` efficiently?
Given that we can compute `P(n)` efficiently, we can use the following formula to compute `G(n)`:
`G(n) = \sum_{k=1}^n \sum_{d|k} d * φ(k/d)`
By flattening the double-summation, we get:
`G(n) = 1/2 \sum_{k=1}^n φ(k) * floor(n/k) * floor(1+n/k)`
It's 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 this requires a fast algorithm for computing the Mertens function (A002321).
However, we can compute partial sums of Euler's totient function more efficiently, without using the Mertens function.
Let:
`Φ(n) = \sum_{k=1}^n φ(k)`
Based on the identity:
`\sum_{d|n} φ(d) = n`
we have:
`\sum_{k=1}^n \sum_{d|k} φ(k) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) φ(j) = \sum_{k=1}^n Φ(floor(n/k)) = \sum_{k=1}^n k = F_1(n)`
implying the following recurrence:
`Φ(n) = F_1(n) - \sum_{k=2}^n Φ(floor(n/k))`
which can be computed in sub-linear time as:
`Φ(n) = F_1(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) Φ(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) Φ(k) * (floor(n/k) - floor(n/(k+1)))`
This requires precomputing `Φ(k)` up to about `k = floor(n^(2/3))` and performing a simple `O(1)` lookup when `k` is less than the precomputed bound.
In combination with the Dirichlet hyperbola method, we can compute `G(n)` in sub-linear time as:
`G(n) = \sum_{k=1}^floor(sqrt(n)) (k * Φ(floor(n/k)) + φ(k) * F_1(floor(n/k))) - Φ(floor(sqrt(n)))*F_1(floor(sqrt(n)))`
where `F_n(x)` are the Faulhaber polynomials.
# Partial sums of Jordan's totient function
Jordan's totient function is a generalization of Euler's totient function and is defined as:
`J_m(n) = \sum_{d|n} μ(d) (n/d)^m = n^m * \prod_{p|n} (1 - 1/p^m)`
where the product is taken over the distinct prime factors of `n`, and is multiplicative with:
`J_m(p^k) = p^(k m) - p^((k - 1) m)`
with the beautiful property:
`\sum_{d|n} J_m(d) = n^m`
Jordan's totient function `J_k(n)` 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 `m=1`, the Jordan totient function is equivalent with Euler's totient function:
`J_1(n) = φ(n)`
`J_m(n) = \sum_{d|n} μ(d) (n/d)^m = n^m * \prod_{p|n} (1 - 1/p^m)`
where the product is taken over the distinct prime factors of `n`, and is multiplicative with:
`J_m(p^k) = p^(k m) - p^((k - 1) m)`
with the beautiful property:
`\sum_{d|n} J_m(d) = n^m`
Jordan's totient function `J_k(n)` 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 `m=1`, the Jordan totient function is equivalent with Euler's totient function:
`J_1(n) = φ(n)`
`\sum_{k=1}^n J_m(k)`
where `J_m(k)` is the Jordan totient function.
There is an elegant formula for computing this partial sum, defined in terms of the Möbius function and Faulhaber's formula:
`\sum_{k=1}^n J_m(k) = \sum_{k=1}^n \sum_{d|k} μ(d) (k/d)^m = \sum_{k=1}^n μ(k) * F_m(floor(n/k))`
where `μ(k)` is the Möbius function and `F_n(x)` are Faulhaber's 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)`
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.
Assuming that we use this algorithm (or a similar one which has roughly the same complexity), then we can also compute this partial sum in sub-linear time:
`\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))`
However, a more practical and efficient approach uses the following recursive formula:
Let:
`Φ_m(n) = \sum_{k=1}^n J_m(k)`
which can be expressed recursively as:
`Φ_m(n) = F_m(n) - \sum_{k=2}^n Φ_m(floor(n/k))`
and computed in sub-linear time as:
`Φ_m(n) = F_m(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) Φ_m(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) Φ_m(k) * (floor(n/k) - floor(n/(k+1)))`
where `F_n(x)` are Faulhaber's polynomials.
This approach requires precomputing `Φ_m(k)` up to about `k = floor(n^(2/3))` and doing a simple lookup when `k` is less than the precomputed bound.
An asymptotic approximation for partial sums of the Jordan's totient function, is given by the following formula:
`\sum_{k=1}^n J_m(k) ~ n^(m+1) / ((m+1) * zeta(m+1))`
or, alternatively:
`\sum_{k=1}^n J_m(k) ~ \frac{F_m(n)}{zeta(m+1)}`
where `zeta(s)` is the Riemann zeta function, for `m >= 1`.
Additionally, the following symmetric identity is satisfied for all 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))`
where `F_n(x)` are Faulhaber's polynomials.
More generally, we also consider the partial sums of `J_m(k) * k^p`, for `m >= 1` and `p >= 0`:
`R(n) = \sum_{k=1}^n J_m(k) * k^p`
for which we have the following beautiful recursive formula:
`R(n) = F_(m+p)(n) - \sum_{k=2}^n k^p * R(floor(n/k))`
which can be computed in sub-linear time as:
`R(n) = F_(m+p)(n) - \sum_{k=2}^floor(n / (1+floor(sqrt(n)))) k^p * R(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) R(k) * (F_p(floor(n/k)) - F_p(floor(n/(k+1))))`
with the asymptotic formula:
`\sum_{k=1}^n J_m(k) * k^p ~ \frac{F_(m+p)(n)}{\zeta(m+1)}`
where `J_m(k)` is the Jordan totient function, `zeta(s)` is the Riemann zeta function and `F_n(x)` are Faulhaber's polynomials, for `m >= 1` and `p >= 0`.
# Partial sums of the Dedekind psi function
The generalized Dedekind psi function is defined as:`ψ_m(n) = \sum_{d|n} |μ(d)| (n/d)^m = (J_(2 m)(n)) / (J_m(n))`
where `μ(k)` is the Möbius function and `J_m(n)` is the Jordan totient function, and is multiplicative with:
`ψ_m(p^k) = p^(k m) + p^((k - 1) m)`
The Möbius transform of the Dedekind psi function can be expressed in terms of the Jordan totient function `J_m(k)`:
`\sum_{d|n} μ(d) ψ_m(n/d) = \sum_{d|n} |μ(d)| J_m(n/d)`
The inverse Möbius transform of the Dedekind psi function satisfies the following identity:
`\sum_{d|n} ψ_m(d) = \sum_{d|n} 2^(ω(d)) (n/d)^m`
which gives us the identity:
`ψ_m(n) = \sum_{d|n} 2^(ω(d)) J_m(n/d)`
where `J_m(k)` is the Jordan totient function and `ω(n)` is the prime omega 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 \sum_{d|k} |μ(d)| (k/d)^m = \sum_{k=1}^n |μ(k)| * F_m(floor(n/k))`
with `F_n(x)` is Faulhaber's formula and `|μ(k)|` is the absolute value of the Möbius function, defined as:
`|μ(k)| = 1` iff `k` is a square-free integer
`|μ(k)| = 0` otherwise
Next we define `S(n)` as the sum over the absolute value of the Möbius function `|μ(k)|` for `1 <= k <= n`, which represents the number of square-free integers between `1` and `n`. Importantly, `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 sub-linear 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))`
By using the Dirichlet hyperbola method, we can rewrite the formula as:
`\sum_{k=1}^n ψ_m(k) = \sum_{k=1}^floor(sqrt(n)) (|μ(k)| * F_m(floor(n/k)) + k^m * S(floor(n/k))) - F_m(floor(sqrt(n)))*S(floor(sqrt(n)))`
where `F_n(x)` are the Faulhaber polynomials.
Additionally, if we define `E_m(n)` and `D_m(n)` as:
`E_m(n) = \sum_{k=1}^n ψ_m(k)`
`D_m(n) = \sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n \sum_{j=1}^floor(n/k) ψ_m(j) = \sum_{k=1}^n E_m(floor(n/k))`
then we can create the following recursive formula to express partial sums of the Dedekind psi function `ψ_m(n)`:
`E_m(n) = D_m(n) - \sum_{k=2}^n E_m(floor(n/k))`
which can be computed in sublinear time if we can compute `D_m(n)` efficiently:
`E_m(n) = D_m(n) - \sum_{k=2}^floor(n/(1+floor(sqrt(n)))) E_m(floor(n/k)) - \sum_{k=1}^floor(sqrt(n)) E_m(k) * (floor(n/k) - floor(n/(k+1)))`
An asymptotic approximation to 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 all 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 inverse Möbius transform of the Dedekind psi function
The inverse Möbius transform of the generalized Dedekind psi function `ψ_m(n)`, is defined as:`\sum_{d|n} ψ_m(d) = \sum_{d|n} 2^(ω(d)) * (n/d)^m`
where `ω(n)` is the prime omega function.
For any fixed integer `m >= 0`, we're interested in computing the following partial sum (A306775):
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n \sum_{d|k} 2^(ω(d)) * (k/d)^m`
which can be expressed in terms of the Faulhaber polynomials `F_n(x)`, as:
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n ψ_m(k) * floor(n/k) = \sum_{k=1}^n 2^(ω(k)) * F_m(floor(n/k))`
which is also equivalent with:
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^n k^m * \sum_{j=1}^floor(n/k) 2^(ω(j))`
By defining `R(n)` as the partial sums of `2^(ω(k))` (see: A064608):
`R(n) = \sum_{k=1}^n 2^(ω(k))`
we can create the following sublinear formula, assuming that we can compute `R(n)` efficiently:
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) = \sum_{k=1}^floor(sqrt(n)) (2^(ω(k)) * F_m(floor(n/k)) + k^m * R(floor(n/k))) - F_m(floor(sqrt(n))) * R(floor(sqrt(n)))`
An asymptotic approximation to this partial sum is given by the following formula:
`\sum_{k=1}^n \sum_{d|k} ψ_m(d) ~ (n^(m+1) * zeta(m+1)^2) / ((m+1) * zeta(2*(m+1)))`
where `zeta(s)` is the Riemann zeta function, for `m >= 1`.
# 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))`
or, alternatively:
`A_m(n) = \sum_{k=1}^floor(sqrt(n)) k^m * \pi(floor(n/k)) + \sum_{p}^floor(sqrt(n)) F_m(floor(n/p)) - F_m(floor(sqrt(n))) * \pi(floor(sqrt(n)))`
where `F_n(x)` are Faulhaber's polynomials, `\pi(x)` is the prime-counting function and `p` runs over the prime numbers.
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)` are Faulhaber's polynomials and `P(s)` is defined as:
`P(s) = \sum_{p} 1/p^s`
with `p` running over the prime numbers.
# Partial sums of the prime big omega function
The prime big omega function, denoted as `Ω(n)`, is defined to be the number of prime factors of `n` counted with multiplicity.We're interested in computing the following partial sum:
`\sum_{k=1}^n Ω(k)`
which is equivalent with:
`\sum_{k=1}^n f(k) * floor(n/k)`
where `f(k)` is defined as:
`f(k) = 1` if `k` is a prime power with exponent `>= 1`
`f(k) = 0` otherwise
Let:
`V(n) = \sum_{k=1}^n f(k)`
This formula can be computed in sub-linear time, because we can use the prime-counting function `\pi(n)` to efficiently compute partial sums of `f(n)`, which represent the number of prime powers `<= n`:
`V(n) = \sum_{k=1}^(floor(log_2(n))) \pi(floor(n^(1/k)))`
allowing us to create the following sub-linear formula:
`\sum_{k=1}^n f(k) * floor(n/k) = \sum_{k=1}^(floor(sqrt(n))) (V(floor(n/k)) - V(floor(n / (k+1)))) * k + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * floor(n/k)`
with the full generalization:
`G_m(n) = \sum_{k=1}^n Ω_m(k) = \sum_{k=1}^n f(k) * F_m(floor(n/k))`
computed efficiently as:
`G_m(n) = \sum_{k=1}^floor(sqrt(n)) (V(floor(n/k)) - V(floor(n / (k+1)))) * F_m(k) + \sum_{k=1}^(floor(n / (floor(sqrt(n)) + 1))) f(k) * F_m(floor(n/k))`
or, alternatively:
`G_m(n) = \sum_{k=1}^floor(sqrt(n)) (k^m * V(floor(n/k)) + f(k) * F_m(floor(n/k))) - F_m(floor(sqrt(n))) * V(floor(sqrt(n)))`
where `F_n(x)` are Faulhaber's polynomials.
Evaluating `G_1(10^n)`, we observe an asymptotic behavior:
`G_1(10^1) = 30`
`G_1(10^2) = 2815`
`G_1(10^3) = 276337`
`G_1(10^4) = 27591490`
`G_1(10^5) = 2758525172`
`G_1(10^6) = 275847515154`
`G_1(10^7) = 27584671195911`
`G_1(10^8) = 2758466558498626`
`G_1(10^9) = 275846649393437566`
`G_1(10^10) = 27584664891073330599`
`G_1(10^11) = 2758466488352698209587`
which can be described as:
`G_1(n) ~ (n (n+1))/2 * 0.55169329765699918...` (see: A154945)
For `m=2` we have:
`G_2(10^1) = 82`
`G_2(10^2) = 66799`
`G_2(10^3) = 64901405`
`G_2(10^4) = 64727468210`
`G_2(10^5) = 64708096890744`
`G_2(10^6) = 64706281936598588`
`G_2(10^7) = 64706077322294843451`
`G_2(10^8) = 64706058761567362618628`
`G_2(10^9) = 64706056807390376400359474`
`G_2(10^10) = 64706056632561375736945155965`
`G_2(10^11) = 64706056612919470606889256184409`
with the asymptotic behavior:
`G_2(n) ~ (n (n + 1) (2 n + 1))/6 * 0.19411816983263379...` (see: A286229)
In general, for `m >= 1`, the asymptotic behavior of `G_m(n)` can be described in terms of the prime zeta function, as:
`G_m(n) ~ F_m(n) * Q(m+1)`
where `F_m(n)` are Faulhaber's polynomials and `Q(m)` is defined as:
`Q(m) = \sum_{k=1}^\infty \sum_{p} 1/p^(m*k) = \sum_{p} 1/(p^m - 1)`
where `p` runs over the prime numbers.
# Partial sums of Dirichlet convolutions of `(-1)^(ω(n))` and `(-1)^(Ω(n))`
In this section, we consider the following partial sums:`\sum_{k=1}^n \sum_{d|k} (-1)^(ω(d)) * (k/d)^m = \sum_{k=1}^n (-1)^(ω(k)) * F_m(floor(n/k))`
and
`\sum_{k=1}^n \sum_{d|k} (-1)^(Ω(d)) * (k/d)^m = \sum_{k=1}^n (-1)^(Ω(k)) * F_m(floor(n/k))`
where `ω(n)` and `Ω(n)` are the prime omega functions, and `F_m(n)` are the Faulhaber polynomials.
Several well-known identities include:
`(-1)^(ω(n)) = μ(rad(n))`
`(-1)^(Ω(n)) = λ(n)`
where `rad(n)` is the radical function, `μ(n)` is the Möbius function and `λ(n)` is the Liouville function.
For a square-free integer `n > 1`, we have:
`\sum_{d|n} (-1)^(ω(d)) = 0`
More generally, if `r^k` is a perfect power, where `r` is square-free and `k >= 2`, then:
`\sum_{d|r^k} (-1)^(ω(d)) = μ(r) * (k-1)^(ω(r))`
For `m=0`, we have the easy-to-compute identity:
`\sum_{k=1}^n (-1)^(Ω(k)) * floor(n/k) = floor(sqrt(n))`
For `m >= 1`, the partial sums can be described asymptotically as:
` \sum_{k=1}^n (-1)^(ω(k)) * F_m(floor(n/k)) ~ F_m(n) * Q(m+1)`
`\sum_{k=1}^n (-1)^(Ω(k)) * F_m(floor(n/k)) ~ F_m(n) * \frac{zeta(2 (m+1))}{zeta(m+1)}`
where `zeta(s)` is the Riemann zeta function, and `Q(s)` is defined as:
`Q(s) = \prod_{p} (1 - 1/(p^s - 1))`
where `p` runs 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:
`ω_m(n) = \sum_{p^k | n} n^m / p^m`
`Ω_m(n) = \sum_{p^k|n} \sum_{j=1}^k n^m / p^(j * m)`
Where for `ω_m(n)` we have the following OEIS sequences:
`ω_0(n)` = A001221(n)
`ω_1(n)` = A069359(n)
`ω_2(n)` = A322078(n)
While for `Ω_m(n)` we have:
`Ω_0(n)` = A001222(n)
`Ω_1(n)` = A095112(n)
`Ω_2(n)` = A322664(n)
The generalized 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)` are Faulhaber's polynomials.
For taking partial sums of the generalized prime omega function `ω_m(k)`, we have:
`\sum_{k=1}^n ω_m(k) = \sum_{p <= n} F_m(floor(n/p))`
where `p` runs over the prime numbers `<= n`.
Very similar, for `Ω_m(k)` we have:
`\sum_{k=1}^n Ω_m(k) = \sum_{p^k <= n} F_m(floor(n/p^k))`
where `p^k` runs over the prime powers `<= n`, with `k >= 1`.
Additionally, we can also consider the following formula, which has some interesting properties:
`a_m(n) = \sum_{p^k | n} (k * n^m) / p^m`
where `p^k` are the prime power factors of `n`, where `n = \prod_{e=1} p_e^(k_e)`.
Several sequences include:
`a_0(n)` = A001222(n)
`a_1(n)` = A003415(n)
`a_2(n)` = A322079(n)
For `m >= 1`, partial sums of this function can be described asymptotically as:
`\sum_{k=1}^n a_m(k) ~ F_m(n) * \sum_{p} 1/(p^m * (p-1))` ; (see: A136141, A152441)
where `F_n(x)` are Faulhaber's polynomials and `p` runs over the prime numbers.
# OEIS sequences
Let's consider the following general partial sum, where `f(k)` is an arithmetic function:`R_m(n) = \sum_{k=1}^n f(k) * F_m(floor(n/k))`
where `F_n(x)` are Faulhaber's polynomials.
Below we have a table of several interesting OEIS sequences for `f(k)`, which produce `R_0(n)` and `R_1(n)`:
+------------+------------+------------+
| `f(k)` | `R_0(n)` | `R_1(n)` |
+------------+------------+------------+
| A000012(k) | A006218(n) | A024916(n) |
| A000027(k) | A024916(n) | A143127(n) |
| A010051(k) | A013939(n) | A322068(n) |
| A008966(k) | A064608(n) | A173290(n) |
| A008683(k) | A000012(n) | A002088(n) |
| A000290(k) | A064602(n) | A143128(n) |
| A000010(k) | A000217(n) | A272718(n) |
| A034444(k) | A061503(n) | A306775(n) |
+------------+------------+------------+
# Implementations in Sidef
This section includes several Sidef scripts, implementing the formulas described in this post:- Generalized prime omega function
- Generalized prime big omega function
- Partial sums of the Möbius function
- Partial sums of the Liouville function
- Partial sums of the sigma function
- Partial sums of the gcd-sum function
- Partial sums of the gcd-sum function fast
- Partial sums of Euler's totient function
- Partial sums of Euler's totient function (recursive)
- Partial sums of Jordan's totient function
- Partial sums of Jordan's totient function times k^m
- Partial sums of Dedekind's psi function
- Partial sums of Dedekind's psi function recursive
- Partial sums of inverse Möbius transform of Dedekind's psi function
- Partial sums of the prime omega function
- Partial sums of the prime bigomega function
- Partial sums of 2^bigomega(k) function
- Partial sums of the least prime factor of n
- Partial sums of the greatest prime factor of n
- Partial sum of the alternating sum of divisors
# Implementations 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 Möbius function
- Partial sums of the Liouville function
- Partial sums of the gcd-sum function
- Partial sums of the gcd-sum function fast
- Partial sums of the gcd-sum function faster
- Partial sums of Euler's totient function
- Partial sums of Euler's totient function faster
- Partial sums of Euler's totient function times k^m
- Partial sums of Jordan's totient function
- Partial sums of Jordan's totient function faster
- Partial sums of Jordan's totient function times k^m
- Partial sums of Dedekind's psi function
- Partial sums of inverse Möbius transform of Dedekind's psi function
- Partial sums of the prime omega function
- Partial sums of the prime bigomega function
- Partial sums of n over k-almost prime divisors
- Partial sums of the least prime factor of n
- Partial sums of the greatest prime factor of n