-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcount_of_k-omega_primes.sf
101 lines (68 loc) · 2.24 KB
/
count_of_k-omega_primes.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
#!/usr/bin/ruby
# Daniel "Trizen" Șuteu
# Date: 14 March 2021
# https://github.com/trizen
# Count the number of k-omega primes <= n.
# Definition:
# k-omega primes are numbers n such that omega(n) = k.
# See also:
# https://en.wikipedia.org/wiki/Almost_prime
# https://en.wikipedia.org/wiki/Prime_omega_function
func omega_prime_count(n, k=1) {
if (k == 1) {
return prime_power_count(n)
}
var count = 0
func (m, p, k, j = 1, s = idiv(n,m).iroot(k)) {
if (k == 2) {
for (true; p <= s; ++j) {
var r = p.next_prime
for (var t = m*p ; t <= n; t *= p) {
var w = idiv(n, t)
break if (p > w)
count += (prime_count(w) - j)
for (var r2 = r; r2 <= w; r2.next_prime!) {
var u = t*r2*r2
break if (u > n)
for (true; u <= n; u *= r2) {
++count
}
}
}
p = r
}
return nil
}
for (true; p <= s; ++j) {
var r = p.next_prime
for (var t = m*p ; t <= n; t *= p) {
var s = idiv(n,t).iroot(k-1)
break if (s < r)
__FUNC__(t, r, k-1, j+1, s)
}
p = r
}
}(1, 2, k)
return count
}
# Run some tests
for k in (1..10) {
var upto = k.pn_primorial+1e5.irand
var x = omega_prime_count(upto, k)
var y = k.omega_primes(upto).len
var z = k.omega_prime_count(upto)
say "Testing: #{k} with n = #{upto} -> #{x}"
assert_eq(x, y)
assert_eq(x, z)
}
say ''
for k in (1..6) {
say ("Count of #{'%2d' % k}-omega primes for 10^n: ", 7.of {|n| omega_prime_count(10**n, k) })
}
__END__
Count of 1-omega primes for 10^n: [0, 7, 35, 193, 1280, 9700, 78734]
Count of 2-omega primes for 10^n: [0, 2, 56, 508, 4097, 33759, 288726]
Count of 3-omega primes for 10^n: [0, 0, 8, 275, 3695, 38844, 379720]
Count of 4-omega primes for 10^n: [0, 0, 0, 23, 894, 15855, 208034]
Count of 5-omega primes for 10^n: [0, 0, 0, 0, 33, 1816, 42492]
Count of 6-omega primes for 10^n: [0, 0, 0, 0, 0, 25, 2285]