danaj / Math-Prime-Util

Perl (XS) module implementing prime number utilities, including sieves
Other
47 stars 21 forks source link

[suggestion] Algorithm for counting k-almost primes <= n. #39

Closed trizen closed 4 years ago

trizen commented 4 years ago

It would be nice to have a new method, called almost_prime_count(n,k), which counts the number of k-almost primes <= n.

Algorithm in Perl:

use 5.020;
use ntheory qw(:all);
use experimental qw(signatures);

sub almost_prime_count ($n, $k) {

    if ($k == 1) {
        return prime_count($n);
    }

    if ($k == 2) {
        return semiprime_count($n);
    }

    my $count = 0;

    sub ($m, $p, $r) {

        my $s = rootint(divint($n, $m), $r);

        if ($r == 2) {
            my $j = prime_count($p) - 2;

            forprimes {
                $count += (prime_count(divint($n, $m * $_)) - ++$j);
            } $p, $s;

            return;
        }

        for (my $q = $p ; $q <= $s ; $q = next_prime($q)) {
            __SUB__->($m * $q, $q, $r - 1);
        }
    }->(1, 2, $k);

    return $count;
}

foreach my $k (1 .. 10) {
    printf("Count of %2d-almost primes <= 10^n: %s\n", $k, join(', ', map { almost_prime_count(powint(10, $_), $k) } 0 .. 10));
}

Note:

See also: https://en.wikipedia.org/wiki/Almost_prime

danaj commented 4 years ago

Great idea.

I need to get the current code released, especially the GMP code. It's been way too long between pushes to CPAN.

danaj commented 4 years ago

I added is_almost_prime(n,k) as a start, to the native int code. It needs optimized PP and some GMP code.

danaj commented 4 years ago

I just added almost_prime_count(n,k) which seems reasonably efficient. I'm testing an nth_almost_prime(n,k) which just reverses using the count, so not the most efficient but looks handy. Some of the OEIS sequences end at 10^14 which this should be able to add a couple terms to. After more testing.

If you have any good results on upper bounds, or better ways to get the nth value that would be great. I haven't found anything obvious.

I still have to do bigint coding.