sidef-scripts/Math/partial_sums_of_sigma_function.sf

#!/usr/bin/ruby

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

# Formula in terms of Bernoulli polynomials for computing the partial-sums of:
#
#   a(n) = Sum_{k=1..n} sigma_j(k)
#
# for some fixed j >= 0, without using the sigma function.

# Formulas:
#   a(n) = Sum_{k=1..n} sigma_j(k)
#        = Sum_{k=1..n} Sum_{d|k} d^j
#        = Sum_{k=1..n} Sum_{b=1..floor(n/k)} b^j
#        = Sum_{k=1..n} (Bernoulli(j+1, 1+floor(n/k)) - Bernoulli(j+1, 0))/(j+1)
#        = floor((n+1)/2) + Sum_{j=1..floor(n/2)} (Bernoulli(j+1, 1+floor(n/k)) - Bernoulli(j+1, 0))/(j+1)
#
# where Bernoulli(n,x) are the Bernoulli polynomials.

# Alternative formulas:
#   a(n) = Sum_{k=1..n} sigma_j(k)
#        = Sum_{k=1..n} k^j * floor(n/k)
#        = (Bernoulli(j+1, n+1) - Bernoulli(j+1, 1+floor(n/2)))/(j+1) + Sum_{k=1..floor(n/2)} k^j * floor(n/k)

# 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, j) {
    sum(1..n, {|k| k.sigma(j) })
}

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

func bernoulli_sigma_partial_sum(n, j) {

    var t = bernoulli(j+1)

     floor((n+1)/2) + sum(1..floor(n/2), {|k|
        (bernoulli(j+1, 1+floor(n / k)) - t) / (j+1)
    })
}

func faulhaber_alternative_sigma_partial_sum(n, j) {
    faulhaber_sum(n, j) - faulhaber_sum(floor(n/2), j) + sum(1..floor(n/2), {|k|
        k**j * floor(n/k)
    })
}

func bernoulli_alternative_sigma_partial_sum(n, j) {
    (bernoulli(j+1, n+1) - bernoulli(j+1, 1+floor(n/2))) / (j+1) + sum(1..floor(n/2), {|k|
        k**j * floor(n/k)
    })
}

#
## Run some tests
#

for j in (1..10) {

    var n = 1000.irand
    var a = sigma_partial_sum(n, j)
    var b = bernoulli_sigma_partial_sum(n, j)
    var c = faulhaber_sigma_partial_sum(n, j)
    var d = faulhaber_alternative_sigma_partial_sum(n, j)
    var e = bernoulli_alternative_sigma_partial_sum(n, j)

    assert_eq(a, b)
    assert_eq(a, c)
    assert_eq(a, d)
    assert_eq(a, e)

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

__END__
Sum_{k=1..439} sigma_1(k) = 158390
Sum_{k=1..606} sigma_2(k) = 89428845
Sum_{k=1..906} sigma_3(k) = 182745921247
Sum_{k=1..158} sigma_4(k) = 20751543588
Sum_{k=1..591} sigma_5(k) = 7261290636306925
Sum_{k=1..840} sigma_6(k) = 42686645719301011492
Sum_{k=1..333} sigma_7(k) = 19205081250504734854
Sum_{k=1..327} sigma_8(k) = 4825111439384649236832
Sum_{k=1..472} sigma_9(k) = 55519692608089816684245494
Sum_{k=1..432} sigma_10(k) = 9008681557899400667788631360