generalized_partial_sums_of_sigma_function.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 08 November 2018
# https://github.com/trizen

# A new general formula for computing the partial-sums of:
#
#   `k^m * sigma_j(k)`, for `1 <= k <= n`
#
# with fixed positive integers `m` and `j`, without using the sigma function.

# Formulas:
#   a(n) = Sum_{k=1..n} k^m * sigma_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^(j+m)
#        = Sum_{k=1..n} k^(m+j) * Sum_{b=1..floor(n/k)} b^m
#        = Sum_{k=1..n} k^(m+j) * (Bernoulli(m+1, 1+floor(n/k)) - Bernoulli(m+1, 0))/(m+1)
#
# where Bernoulli(n,x) are the Bernoulli polynomials.

# See also:
#   https://en.wikipedia.org/wiki/Divisor_function
#   https://en.wikipedia.org/wiki/Faulhaber's_formula
#   https://en.wikipedia.org/wiki/Bernoulli_polynomials
#   https://trizenx.blogspot.com/2018/08/interesting-formulas-and-exercises-in.html

func sigma_partial_sum(n, m, j) {
    sum(1..n, {|k| k**m * k.sigma(j) })
}

func faulhaber_sigma_partial_sum(n, m, j) {
    sum(1..n, {|k|
         k**(m + j) * faulhaber_sum(floor(n / k), m)
    })
}

func bernoulli_sigma_partial_sum(n, m, j) {

    var b = bernoulli(m+1)

    sum(1..n, {|k|
        k**(m + j) * (bernoulli(m+1, 1+floor(n / k)) - b) / (m+1)
    })
}

#
## Run some tests
#

for m in (1..5) {

    for j in (1..3) {

        var n = 100.irand
        var a = sigma_partial_sum(n, m, j)
        var b = bernoulli_sigma_partial_sum(n, m, j)
        var c = faulhaber_sigma_partial_sum(n, m, j)

        assert_eq(a, b)
        assert_eq(a, c)

        say "Sum_{k=1..#{n}} k^#{m} * sigma_#{j}(k) = #{b}"
    }
}

__END__
Sum_{k=1..51} k^1 * sigma_1(k) = 74034
Sum_{k=1..50} k^1 * sigma_2(k) = 1959930
Sum_{k=1..32} k^1 * sigma_3(k) = 7871069
Sum_{k=1..66} k^2 * sigma_1(k) = 8111384
Sum_{k=1..72} k^2 * sigma_2(k) = 485826460
Sum_{k=1..62} k^2 * sigma_3(k) = 10762700836
Sum_{k=1..84} k^3 * sigma_1(k) = 1434314059
Sum_{k=1..68} k^3 * sigma_2(k) = 20740625414
Sum_{k=1..22} k^3 * sigma_3(k) = 453251592
Sum_{k=1..45} k^4 * sigma_1(k) = 2436788727
Sum_{k=1..34} k^4 * sigma_2(k) = 10012218107
Sum_{k=1..67} k^4 * sigma_3(k) = 58038503678200
Sum_{k=1..41} k^5 * sigma_1(k) = 47926340149
Sum_{k=1..92} k^5 * sigma_2(k) = 808713405614655
Sum_{k=1..93} k^5 * sigma_3(k) = 65544349176333209