/
partial_sums_of_exponential_prime_omega_functions.pl
194 lines (156 loc) · 4.83 KB
/
partial_sums_of_exponential_prime_omega_functions.pl
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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
#!/usr/bin/perl
# Daniel "Trizen" Șuteu
# Date: 16 March 2021
# https://github.com/trizen
# Compute partial sums of the following three functions in sublinear time:
# S1(n) = Sum_{k=1..n} v^bigomega(k)
# S2(n) = Sum_{k=1..n} v^omega(k)
# S3(n) = Sum_{k=1..n} v^omega(k) * mu(k)^2
use 5.020;
use warnings;
use ntheory qw(:all);
use experimental qw(signatures);
sub squarefree_almost_prime_count ($k, $n) {
if ($k == 0) {
return (($n <= 0) ? 0 : 1);
}
if ($k == 1) {
return prime_count($n);
}
my $count = 0;
sub ($m, $p, $k, $j = 0) {
my $s = rootint(divint($n, $m), $k);
if ($k == 2) {
foreach my $q (@{primes($p, $s)}) {
++$j;
if (modint($m, $q) != 0) {
$count += prime_count(divint($n, mulint($m, $q))) - $j;
}
}
return;
}
foreach my $p (@{primes($p, $s)}) {
if (modint($m, $p) != 0) {
__SUB__->(mulint($m, $p), $p, $k - 1, $j);
}
++$j;
}
}->(1, 2, $k);
return $count;
}
sub S1 ($n, $v = 2) { # Sum_{k=1..n} v^bigomega(k)
vecsum(map { mulint(powint($v, $_), almost_prime_count($_, $n)) } 0 .. logint($n, 2));
}
sub S2 ($n, $v = 2) { # Sum_{k=1..n} v^omega(k)
vecsum(map { mulint(powint($v, $_), omega_prime_count($_, $n)) } 0 .. logint($n, 2));
}
sub S3 ($n, $v = 2) { # Sum_{k=1..n} v^omega(k) * mu(k)^2
vecsum(map { mulint(powint($v, $_), squarefree_almost_prime_count($_, $n)) } 0 .. logint($n, 2));
}
say join ', ', map { S1($_) } 1 .. 20; #=> A069205: [1, 3, 5, 9, 11, 15, 17, 25, 29, 33, 35, 43, 45, 49, 53, 69, 71, 79, 81]
say join ', ', map { S2($_) } 1 .. 20; #=> A064608: [1, 3, 5, 7, 9, 13, 15, 17, 19, 23, 25, 29, 31, 35, 39, 41, 43, 47, 49]
say join ', ', map { S3($_) } 1 .. 20; #=> A069201: [1, 3, 5, 5, 7, 11, 13, 13, 13, 17, 19, 19, 21, 25, 29, 29, 31, 31, 33]
say '';
say join ', ', map { S1($_, -1) } 1 .. 20; #=> A002819: [1, 0, -1, 0, -1, 0, -1, -2, -1, 0, -1, -2, -3, -2, -1, 0, -1, -2, -3, -4]
say join ', ', map { S2($_, -1) } 1 .. 20; #=> A174863: [1, 0, -1, -2, -3, -2, -3, -4, -5, -4, -5, -4, -5, -4, -3, -4, -5, -4, -5, -4]
say join ', ', map { S3($_, -1) } 1 .. 20; #=> A002321: [1, 0, -1, -1, -2, -1, -2, -2, -2, -1, -2, -2, -3, -2, -1, -1, -2, -2, -3, -3]
__END__
# A069205(n) = Sum_{k=1..n} 2^bigomega(k)
A069205(10^1) = 33
A069205(10^2) = 811
A069205(10^3) = 15301
A069205(10^4) = 260615
A069205(10^5) = 3942969
A069205(10^6) = 55282297
A069205(10^7) = 746263855
A069205(10^8) = 9613563919
A069205(10^9) = 120954854741
A069205(10^10) = 1491898574939
A069205(10^11) = 17944730372827
A069205(10^12) = 212986333467973
A069205(10^13) = 2498962573520227
A069205(10^14) = 28874142998632109
# A002819(n) = Sum_{k=1..n} (-1)^bigomega(k)
# See also: A090410
A002819(10^1) = 0
A002819(10^2) = -2
A002819(10^3) = -14
A002819(10^4) = -94
A002819(10^5) = -288
A002819(10^6) = -530
A002819(10^7) = -842
A002819(10^8) = -3884
A002819(10^9) = -25216
A002819(10^10) = -116026
A002819(10^11) = -342224
A002819(10^12) = -522626
A002819(10^13) = -966578
A002819(10^14) = -7424752
# A064608(n) = Sum_{k=1..n} 2^omega(k)
# See also: A180361
A064608(10^1) = 23
A064608(10^2) = 359
A064608(10^3) = 4987
A064608(10^4) = 63869
A064608(10^5) = 778581
A064608(10^6) = 9185685
A064608(10^7) = 105854997
A064608(10^8) = 1198530315
A064608(10^9) = 13385107495
A064608(10^10) = 147849112851
A064608(10^11) = 1618471517571
A064608(10^12) = 17584519050293
# A174863(n) = Sum_{k=1..n} (-1)^omega(k)
A174863(10^1) = -4
A174863(10^2) = 14
A174863(10^3) = 64
A174863(10^4) = -16
A174863(10^5) = -720
A174863(10^6) = -1908
A174863(10^7) = -1650
A174863(10^8) = 10734
A174863(10^9) = 53740
A174863(10^10) = 108654
A174863(10^11) = 195702
A174863(10^12) = 27158
# A069201(n) = Sum_{k=1..n} mu(k)^2 * 2^omega(k)
A069201(10^1) = 17
A069201(10^2) = 211
A069201(10^3) = 2825
A069201(10^4) = 34891
A069201(10^5) = 414813
A069201(10^6) = 4808081
A069201(10^7) = 54684335
A069201(10^8) = 612868643
A069201(10^9) = 6788951097
A069201(10^10) = 74492096539
A069201(10^11) = 810947010335
A069201(10^12) = 8769730440341
# A002321(n) = Sum_{k=1..n} (-1)^omega(k) * mu(k)^2 = Sum_{k=1..n} mu(k)
# See also: A084237
A002321(10^1) = -1
A002321(10^2) = 1
A002321(10^3) = 2
A002321(10^4) = -23
A002321(10^5) = -48
A002321(10^6) = 212
A002321(10^7) = 1037
A002321(10^8) = 1928
A002321(10^9) = -222
A002321(10^10) = -33722
A002321(10^11) = -87856
A002321(10^12) = 62366
# A013928(n) = Sum_{k=1..n} mu(k)^2
# See also: A071172
A013928(10^1) = 7
A013928(10^2) = 61
A013928(10^3) = 608
A013928(10^4) = 6083
A013928(10^5) = 60794
A013928(10^6) = 607926
A013928(10^7) = 6079291
A013928(10^8) = 60792694
A013928(10^9) = 607927124
A013928(10^10) = 6079270942
A013928(10^11) = 60792710280
A013928(10^12) = 607927102274