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