|
| 1 | +#!/usr/bin/perl |
| 2 | + |
| 3 | +# Daniel "Trizen" Șuteu |
| 4 | +# Date: 14 March 2021 |
| 5 | +# Edit: 25 March 2025 |
| 6 | +# https://github.com/trizen |
| 7 | + |
| 8 | +# Generate k-omega primes in range [a,b]. (not in sorted order) |
| 9 | + |
| 10 | +# Definition: |
| 11 | +# k-omega primes are numbers n such that omega(n) = k. |
| 12 | + |
| 13 | +# See also: |
| 14 | +# https://en.wikipedia.org/wiki/Almost_prime |
| 15 | +# https://en.wikipedia.org/wiki/Prime_omega_function |
| 16 | + |
| 17 | +use 5.020; |
| 18 | +use integer; |
| 19 | +use ntheory qw(:all); |
| 20 | +use experimental qw(signatures); |
| 21 | + |
| 22 | +sub omega_prime_numbers ($A, $B, $k, $callback) { |
| 23 | + |
| 24 | + $A = vecmax($A, pn_primorial($k)); |
| 25 | + |
| 26 | + sub ($m, $p, $k) { |
| 27 | + |
| 28 | + my $s = rootint($B / $m, $k); |
| 29 | + |
| 30 | + foreach my $q (@{primes($p, $s)}) { |
| 31 | + |
| 32 | + my $r = next_prime($q); |
| 33 | + |
| 34 | + for (my $v = $m * $q ; $v <= $B ; $v *= $q) { |
| 35 | + if ($k == 1) { |
| 36 | + $callback->($v) if ($v >= $A); |
| 37 | + } |
| 38 | + elsif ($v * $r <= $B) { |
| 39 | + __SUB__->($v, $r, $k - 1); |
| 40 | + } |
| 41 | + } |
| 42 | + } |
| 43 | + }->(1, 2, $k); |
| 44 | +} |
| 45 | + |
| 46 | +# Generate 5-omega primes in the range [3000, 10000] |
| 47 | + |
| 48 | +my $k = 5; |
| 49 | +my $from = 3000; |
| 50 | +my $upto = 10000; |
| 51 | + |
| 52 | +my @arr; |
| 53 | +omega_prime_numbers($from, $upto, $k, sub ($n) { push @arr, $n }); |
| 54 | + |
| 55 | +my @test = grep { prime_omega($_) == $k } $from .. $upto; # just for testing |
| 56 | +join(' ', sort { $a <=> $b } @arr) eq join(' ', @test) or die "Error: not equal!"; |
| 57 | + |
| 58 | +say join(', ', @arr); |
| 59 | + |
| 60 | +# Run some tests |
| 61 | + |
| 62 | +foreach my $k (1 .. 6) { |
| 63 | + |
| 64 | + my $from = pn_primorial($k) + int(rand(1e4)); |
| 65 | + my $upto = $from + int(rand(1e5)); |
| 66 | + |
| 67 | + say "Testing: $k with $from .. $upto"; |
| 68 | + |
| 69 | + my @arr; |
| 70 | + omega_prime_numbers($from, $upto, $k, sub ($n) { push @arr, $n }); |
| 71 | + |
| 72 | + my @test = grep { prime_omega($_) == $k } $from .. $upto; |
| 73 | + join(' ', sort { $a <=> $b } @arr) eq join(' ', @test) or die "Error: not equal!"; |
| 74 | +} |
0 commit comments