sidef-scripts/Math/partial_sums_of_2_to_the_bigomega_of_n.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 22 May 2020
# Edit: 11 July 2020
# https://github.com/trizen

# Fast formulas for computing: Sum_{k=1..n} 2^bigomega(k).

# OEIS:
#   https://oeis.org/A069205 -- a(n) = Sum_{k=1..n} 2^bigomega(k).

# Let:
#   S(n) = Sum_{k=1..n} 2^bigomega(k)

# Is it possible to compute S(n) in sublinear time? The answer is "yes".

# Formula #1:
#
#   S(n) = 1 + Sum_{k=1..floor(log_2(n))} 2^k * pi_k(n)
#
# where pi_k(n) is the number of k-almost primes <= n.

# Formula #2:
#
#   S(n) = Sum_{2-powerful k <= n} 2^(bigomega(k) - 2*omega(k)) * D(floor(n/k))
#
# where D(n) = Sum_{k=1..n} floor(n/k), computed in O(sqrt(n)) steps as:
#
#   D(n) = 2*Sum_{k=1..floor(sqrt(n))} floor(n/k) - floor(sqrt(n))^2

# See also:
#   https://projecteuler.net/problem=708

func S(n) {
    1 + sum(1..n.ilog2, {|k|
        2**k * k.almost_primepi(n)
    })
}

func T(n) {
    2.powerful(n).sum {|k|
        2**(bigomega(k) - 2*omega(k)) * dirichlet_hyperbola(
            floor(n/k), { 1 }, { 1 }, { _ }, { _ }
        )
    }
}

assert_eq(
    1..100->map{|k| 2**bigomega(k) }.accumulate,
    1..100->map(S),
)

assert_eq(
    1..100->map(S),
    1..100->map(T),
)

for k in (1..6) {

    var a = S(10**k)
    var b = T(10**k)

    assert_eq(a,b)

    say "S(10^#{k}) = #{a}"
}

__END__
S(10^1) = 33
S(10^2) = 811
S(10^3) = 15301
S(10^4) = 260615
S(10^5) = 3942969
S(10^6) = 55282297
S(10^7) = 746263855
S(10^8) = 9613563919
S(10^9) = 120954854741