123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596 |
- #!/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
|