1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495 |
- #!/usr/bin/ruby
- # Author: Trizen
- # Date: 11 May 2024
- # https://github.com/trizen
- # A sublinear formula for computing the partial sums alternating sum of divisors function (A206369), using Dirichlet's hyperbola method.
- # Formula:
- #
- # a(n) = Sum_{k=1..n} Sum_{d|k} (-1)^bigomega(k/d) * d
- # = Sum_{k=1..n} Sum_{d|k, d is a perfect square} phi(k/d)
- # = (1/2) * Sum_{k=1..n} (-1)^bigomega(k) * floor(n/k) * floor(n/k + 1)
- # = (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) = Sum_{k=1..n} phi(k), where `phi` is the Euler totient function.
- # Example:
- # a(10^1) = 35
- # a(10^2) = 3311
- # a(10^3) = 329283
- # a(10^4) = 32900704
- # a(10^5) = 3289890783
- # a(10^6) = 328986877961
- # a(10^7) = 32898683779520
- # a(10^8) = 3289868149892131
- # a(10^9) = 328986813691207320
- # a(10^10) = 32898681337806276406
- # Asymptotic formula due to Tóth, 2013:
- # a(n) ~ (Pi^2/30) * n^2.
- # OEIS sequences:
- # https://oeis.org/A370905 -- Partial sums of the alternating sum of divisors function (A206369).
- # https://oeis.org/A002088 -- Sum of totient function: a(n) = Sum_{k=1..n} phi(k), cf. A000010.
- # See also:
- # https://en.wikipedia.org/wiki/Dirichlet_hyperbola_method
- # https://trizenx.blogspot.com/2018/11/partial-sums-of-arithmetical-functions.html
- func euler_totient_partial_sum (n) { # using Dirichlet's hyperbola method
- var total = 0
- var s = n.isqrt
- for k in (1..s) {
- total += k*mertens(idiv(n,k))
- }
- s.each_squarefree {|k|
- total += moebius(k)*faulhaber(idiv(n,k), 1)
- }
- total -= mertens(s)*faulhaber(s, 1)
- return total
- }
- func partial_sums_of_asod_1(n) { # sublinear complexity
- func S(n) {
- return euler_totient_partial_sum(n)
- }
- sum(1..n.iroot(4), {|k|
- S(floor(n / k**2))
- }) + sum(1..n.isqrt, {|k|
- phi(k) * isqrt(floor(n/k))
- }) - S(n.isqrt)*n.iroot(4)
- }
- func partial_sums_of_asod_2(n) { # sublinear time (slow-ish)
- sum(1..n.isqrt, {|k|
- k*liouville_sum(idiv(n,k)) + liouville(k)*faulhaber_sum(idiv(n,k), 1)
- }) - liouville_sum(n.isqrt)*faulhaber_sum(n.isqrt, 1)
- }
- func partial_sums_of_asod_test(n) { # just for testing
- sum(1..n, {|k| liouville(k) * faulhaber(idiv(n,k), 1) })
- }
- for m in (0..10) {
- var n = 1000.irand
- var t1 = partial_sums_of_asod_1(n)
- var t2 = partial_sums_of_asod_2(n)
- var t3 = partial_sums_of_asod_test(n)
- assert_eq(t1, t2, [n, t1, t2])
- assert_eq(t1, t3, [n, t1, t3])
- say "Sum_{k=1..#{n}} Sum_{d|k} λ(k/d) * d = #{t2}"
- }
|