-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpartial_sums_of_prime_omega_function.sf
151 lines (121 loc) · 4.21 KB
/
partial_sums_of_prime_omega_function.sf
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#!/usr/bin/ruby
# Author: Daniel "Trizen" Șuteu
# Date: 25 November 2018
# https://github.com/trizen
# A new algorithm for computing the partial-sums of the prime omega function `ω(k)`, for `1 <= k <= n`:
# Sum_{k=1..n} ω(k)
# Based on the formula:
# Sum_{k=1..n} ω(k) = Sum_{p prime <= n} floor(n/p)
# Generalization:
# A_m(n) = Sum_{p prime <= n} F_m(floor(n/p))
#
# where F_m(x) is Faulhaber's formula for `Sum_{k=1..x} k^m`.
# See also:
# https://oeis.org/A013939
# https://oeis.org/A064182
# https://oeis.org/A069359
# https://oeis.org/A322068
# https://en.wikipedia.org/wiki/Prime_zeta_function
# https://en.wikipedia.org/wiki/Prime_omega_function
# https://en.wikipedia.org/wiki/Prime-counting_function
# https://trizenx.blogspot.com/2018/08/interesting-formulas-and-exercises-in.html
# Example for `m=0` (A064182):
# A_0(10^1) = 11
# A_0(10^2) = 171
# A_0(10^3) = 2126
# A_0(10^4) = 24300
# A_0(10^5) = 266400
# A_0(10^6) = 2853708
# A_0(10^7) = 30130317
# A_0(10^8) = 315037281
# A_0(10^9) = 3271067968
# A_0(10^10) = 33787242719
# A_0(10^11) = 347589015681
# A_0(10^12) = 3564432632541
# Example for `m=1`:
# A_1(10^1) = 25
# A_1(10^2) = 2298
# A_1(10^3) = 226342
# A_1(10^4) = 22616110
# A_1(10^5) = 2261266482
# A_1(10^6) = 226124236118
# A_1(10^7) = 22612374197143
# A_1(10^8) = 2261237139656553
# A_1(10^9) = 226123710243814636
# A_1(10^10) = 22612371006991736766
# A_1(10^11) = 2261237100241987653515
# A_1(10^12) = 226123710021083492369813
# A_1(10^13) = 22612371002056432695022703
# A_1(10^14) = 2261237100205367824451036203
# Example for `m=2`:
# A_2(10^1) = 75
# A_2(10^2) = 59962
# A_2(10^3) = 58403906
# A_2(10^4) = 58270913442
# A_2(10^5) = 58255785988898
# A_2(10^6) = 58254390385024132
# A_2(10^7) = 58254229074894448703
# A_2(10^8) = 58254214780225801032503
# A_2(10^9) = 58254213248247357411667320
# A_2(10^10) = 58254213116747777047390609694
# A_2(10^11) = 58254213101385832019517484266265
# A_2(10^12) = 58254213099991292350208499967189227
# A_2(10^13) = 58254213099830361065330294973944269431
# Asymptotic formulas:
# A_1(n) ~ 0.4522474200410654985065... * n*(n+1)/2 (see: https://oeis.org/A085548)
# A_2(n) ~ 0.1747626392994435364231... * n*(n+1)*(2*n+1)/6 (see: https://oeis.org/A085541)
# For `m >= 1`, `A_m(n)` can be described asymptotically in terms of the prime zeta function:
# A_m(n) ~ F_m(n) * P(m+1)
#
# where P(s) is defined as:
# P(s) = Sum_{p prime >= 2} 1/p^s
func prime_omega_partial_sum (n, m=0) { # O(sqrt(n)) complexity
var total = 0
var s = n.isqrt
var u = floor(n/(s+1))
for k in (1..s) {
total += (faulhaber_sum(k, m) * prime_count(floor(n/(k+1))+1, floor(n/k)))
}
for p in (primes(u)) {
total += faulhaber_sum(floor(n/p), m)
}
return total
}
func prime_omega_partial_sum_2 (n, m=0) { # O(sqrt(n)) complexity
var total = 0
var s = n.isqrt
for k in (1..s) {
total += (k**m * prime_count(floor(n/k)))
total += faulhaber_sum(floor(n/k), m) if k.is_prime
}
total -= prime_count(s)*faulhaber_sum(s, m)
return total
}
func prime_omega_partial_sum_test (n, m=0) { # just for testing
var total = 0
for p in (primes(n)) {
total += faulhaber_sum(floor(n/p), m)
}
return total
}
for m in (0 .. 10) {
var n = 10000.irand
var t1 = prime_omega_partial_sum(n, m)
var t2 = prime_omega_partial_sum_2(n, m)
var t3 = prime_omega_partial_sum_test(n, m)
assert_eq(t1, t2)
assert_eq(t1, t3)
say "Sum_{k=1..#{n}} omega_#{m}(k) = #{t1}"
}
__END__
Sum_{k=1..2526} omega_0(k) = 5708
Sum_{k=1..5691} omega_1(k) = 7324853
Sum_{k=1..1867} omega_2(k) = 379161614
Sum_{k=1..3539} omega_3(k) = 3018968330718
Sum_{k=1..3152} omega_4(k) = 2227743127912110
Sum_{k=1..3789} omega_5(k) = 8419989859247261517
Sum_{k=1..8477} omega_6(k) = 3722366146862937902377096
Sum_{k=1..5670} omega_7(k) = 543094671659948787399171188
Sum_{k=1..3126} omega_8(k) = 6366824727060827116574504619
Sum_{k=1..5202} omega_9(k) = 1444635343895056823620338084965462
Sum_{k=1..7271} omega_10(k) = 134840274708837491442524789498005113305