123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102 |
- #!/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]
|