-
Notifications
You must be signed in to change notification settings - Fork 33
/
partial_sums_of_prime_omega_function.pl
executable file
·155 lines (123 loc) · 4.36 KB
/
partial_sums_of_prime_omega_function.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
#!/usr/bin/perl
# Author: Daniel "Trizen" Șuteu
# Date: 24 November 2018
# https://github.com/trizen
# A new algorithm for computing the partial-sums of the generalized prime omega function `ω_m(k)`, for `1 <= k <= n`:
# A_m(n) = Sum_{k=1..n} ω_m(k)
#
# where:
# ω_m(n) = n^m * Sum_{p|n} 1/p^m
# Based on the formula:
# Sum_{k=1..n} ω_m(k) = Sum_{p prime <= n} F_m(floor(n/p))
#
# where F_n(x) is Faulhaber's formula.
# Example for `m=0`:
# 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
# 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
# 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
# OEIS sequences:
# https://oeis.org/A013939 -- Partial sums of sequence A001221 (number of distinct primes dividing n).
# https://oeis.org/A064182 -- Sum_{k <= 10^n} number of distinct primes dividing k.
# See also:
# 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
use 5.020;
use strict;
use warnings;
use experimental qw(signatures);
use Math::AnyNum qw(faulhaber_sum ipow);
use ntheory qw(forprimes prime_count sqrtint is_prime);
sub prime_omega_partial_sum ($n, $m) { # O(sqrt(n)) complexity
my $total = 0;
my $s = sqrtint($n);
my $u = int($n / ($s + 1));
for my $k (1 .. $s) {
$total += faulhaber_sum($k, $m) * prime_count(int($n/($k+1))+1, int($n/$k));
}
forprimes {
$total += faulhaber_sum(int($n/$_), $m);
} $u;
return $total;
}
sub prime_omega_partial_sum_2 ($n, $m) { # O(sqrt(n)) complexity
my $total = 0;
my $s = sqrtint($n);
for my $k (1 .. $s) {
$total += ipow($k, $m) * prime_count(int($n/$k));
$total += faulhaber_sum(int($n/$k), $m) if is_prime($k);
}
$total -= faulhaber_sum($s, $m) * prime_count($s);
return $total;
}
sub prime_omega_partial_sum_test ($n, $m) { # just for testing
my $total = 0;
forprimes {
$total += faulhaber_sum(int($n/$_), $m);
} $n;
return $total;
}
for my $m (0 .. 10) {
my $n = int rand 100000;
my $t1 = prime_omega_partial_sum($n, $m);
my $t2 = prime_omega_partial_sum_2($n, $m);
my $t3 = prime_omega_partial_sum_test($n, $m);
die "error: $t1 != $t2" if ($t1 != $t2);
die "error: $t1 != $t3" if ($t1 != $t3);
say "Sum_{k=1..$n} omega_$m(k) = $t1";
}
__END__
Sum_{k=1..93178} omega_0(k) = 247630
Sum_{k=1..60545} omega_1(k) = 828906439
Sum_{k=1..61222} omega_2(k) = 13368082621946
Sum_{k=1..58175} omega_3(k) = 220463446471253532
Sum_{k=1..26576} omega_4(k) = 94816277435320229002
Sum_{k=1..17978} omega_5(k) = 96085844643312478233603
Sum_{k=1..99336} omega_6(k) = 112956550182103434253591001302255
Sum_{k=1..15217} omega_7(k) = 1459563487599016502195229269710
Sum_{k=1..62565} omega_8(k) = 3271462737352430519765722633491562894793
Sum_{k=1..91318} omega_9(k) = 4007044838270388920307792726568428120477189405
Sum_{k=1..28834} omega_10(k) = 514524955177931497535073881648700561462698676