sidef-scripts/Math/partial_sums_of_prime_bigomega_function.sf

#!/usr/bin/ruby

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

# A nice algorithm in terms of the prime-counting function for computing partial sums of the generalized bigomega(n) function:
#   B_m(n) = Sum_{k=1..n} Ω_m(k)

# For `m=0`, we have:
#   B_0(n) = bigomega(n!)

# OEIS related sequences:
#   https://oeis.org/A025528
#   https://oeis.org/A022559
#   https://oeis.org/A071811
#   https://oeis.org/A154945  (0.55169329765699918...)
#   https://oeis.org/A286229  (0.19411816983263379...)

# Example for `B_0(n)`:
#    B_0(10^1) = 15
#    B_0(10^2) = 239
#    B_0(10^3) = 2877
#    B_0(10^4) = 31985
#    B_0(10^5) = 343614
#    B_0(10^6) = 3626619
#    B_0(10^7) = 37861249
#    B_0(10^8) = 392351272
#    B_0(10^9) = 4044220058
#    B_0(10^10) = 41518796555
#    B_0(10^11) = 424904645958

# Example for `B_1(n)`:
#   B_1(10^1) = 30
#   B_1(10^2) = 2815
#   B_1(10^3) = 276337
#   B_1(10^4) = 27591490
#   B_1(10^5) = 2758525172
#   B_1(10^6) = 275847515154
#   B_1(10^7) = 27584671195911
#   B_1(10^8) = 2758466558498626
#   B_1(10^9) = 275846649393437566
#   B_1(10^10) = 27584664891073330599
#   B_1(10^11) = 2758466488352698209587

# Example for `B_2(n)`:
#   B_2(10^1) = 82
#   B_2(10^2) = 66799
#   B_2(10^3) = 64901405
#   B_2(10^4) = 64727468210
#   B_2(10^5) = 64708096890744
#   B_2(10^6) = 64706281936598588
#   B_2(10^7) = 64706077322294843451
#   B_2(10^8) = 64706058761567362618628
#   B_2(10^9) = 64706056807390376400359474
#   B_2(10^10) = 64706056632561375736945155965
#   B_2(10^11) = 64706056612919470606889256184409

# Asymptotic formulas:
#   B_1(n) ~ 0.55169329765699918... * n*(n+1)/2
#   B_2(n) ~ 0.19411816983263379... * n*(n+1)*(2*n+1)/6

# In general, for `m>=1`, we have the following asymptotic formula:
#   B_m(n) ~ (Sum_{k>=1} primezeta((m+1)*k)) * F_m(n)
#
# where F_n(x) is Faulhaber's formula and primezeta(s) is the prime zeta function.

# The prime zeta function is defined as:
#   primezeta(s) = Sum_{p prime >= 2} 1/p^s

# See also:
#   https://oeis.org/A013939
#   https://oeis.org/A064182
#   https://en.wikipedia.org/wiki/Prime_omega_function
#   https://en.wikipedia.org/wiki/Prime-counting_function
#   https://trizenx.blogspot.com/2018/11/partial-sums-of-arithmetical-functions.html

func prime_bigomega_partial_sum(n, m) {     # O(sqrt(n)) complexity

    var s = n.isqrt
    var u = floor(n/(s+1))

    var total = 0

    for k in (1..s) {
        total += faulhaber_sum(k, m)*(prime_power_count(floor(n/(k+1))+1, floor(n/k)))
    }

    for k in (1..u) {
        total += faulhaber_sum(floor(n/k), m) if k.is_prime_power
    }

    return total
}

func prime_bigomega_partial_sum_2(n, m) {

    var s = n.isqrt
    var total = 0

    for k in (1..s) {
        total += (k**m * prime_power_count(floor(n/k)))
        total += faulhaber_sum(floor(n/k), m) if k.is_prime_power
    }

    total -= prime_power_count(s)*faulhaber_sum(s, m)

    return total
}

func prime_bigomega_partial_sum_test (n, m) {    # just for testing
    var total = 0

    for k in (1..n) {
        total += faulhaber_sum(floor(n/k), m) if k.is_prime_power
    }

    return total
}

for m in (0 .. 10) {

    var n = 10000.irand

    var t1 = prime_bigomega_partial_sum(n, m)
    var t2 = prime_bigomega_partial_sum_2(n, m)
    var t3 = prime_bigomega_partial_sum_test(n, m)

    assert_eq(t1, t2)
    assert_eq(t1, t3)

    say "Sum_{k=1..#{n}} bigomega_#{m}(k) = #{t1}"
}

__END__
Sum_{k=1..6682} bigomega_0(k) = 21048
Sum_{k=1..3447} bigomega_1(k) = 3277974
Sum_{k=1..9287} bigomega_2(k) = 51825111802
Sum_{k=1..9418} bigomega_3(k) = 159996950760258
Sum_{k=1..7808} bigomega_4(k) = 213588687883972941
Sum_{k=1..1349} bigomega_5(k) = 17394189051711781
Sum_{k=1..6124} bigomega_6(k) = 385556244990136505278441
Sum_{k=1..8317} bigomega_7(k) = 11667206972644664046687556070
Sum_{k=1..5282} bigomega_8(k) = 715286720988907445715200810310
Sum_{k=1..7094} bigomega_9(k) = 32146805184233019106400204605763207
Sum_{k=1..1863} bigomega_10(k) = 42157952745357988027790444565086