-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpartial_sums_of_prime_bigomega_function.sf
149 lines (117 loc) · 3.96 KB
/
partial_sums_of_prime_bigomega_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
#!/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 += (ipow(k,m) * prime_power_count(floor(n/k)))
}
each_prime_power(1, s, {|k|
total += faulhaber_sum(floor(n/k), m)
})
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