Skip to content

Commit

Permalink
- Added the Number k.squarefree_almost_prime_count(a,b).
Browse files Browse the repository at this point in the history 
Returns the number of squarefree k-almost primes in the range a..b.

Example:

    say 1.squarefree_almost_prime_count(1000)          # count of primes <= 1000
    say 2.squarefree_almost_prime_count(50, 1000)      # count of squarefree semiprimes in range 50..1000

Also aliased as `squarefree_pi_k(k, a, b)`.
  • Loading branch information
@trizen
trizen committed Mar 15, 2021
1 parent 0bedc21 commit a236a25
Show file tree
Hide file tree
Showing 3 changed files with 148 additions and 26 deletions.
135 changes: 114 additions & 21 deletions lib/Sidef/Types/Number/Number.pm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10405,8 +10405,7 @@ package Sidef::Types::Number::Number {

return ZERO if ($y < $x);

my $count = _prime_count($x, $y);
_set_int($count);
_set_int(_prime_count($x, $y));
}

*primepi = \&prime_count;
Expand Down Expand Up @@ -10661,14 +10660,14 @@ package Sidef::Types::Number::Number {
sub {
my ($m, $p, $k, $j) = @_;

Math::GMPz::Rmpz_div($t, $n, $m);
my $s = do {
Math::GMPz::Rmpz_div($t, $n, $m);
Math::GMPz::Rmpz_root($t, $t, $k);
Math::GMPz::Rmpz_get_ui($t);
};

if ($k == 2) {

Math::GMPz::Rmpz_sqrt($t, $t);

my $s = Math::GMPz::Rmpz_get_str($t, 10);

for (
my $q = $p ;
$q <= $s ;
Expand Down Expand Up @@ -10699,9 +10698,15 @@ package Sidef::Types::Number::Number {
return;
}

Math::GMPz::Rmpz_root($t, $t, $k);

foreach my $q (Math::Prime::Util::GMP::sieve_primes($p, Math::GMPz::Rmpz_get_str($t, 10))) {
for (
my $q = $p ;
$q <= $s ;
$q = (
HAS_PRIME_UTIL
? Math::Prime::Util::next_prime($q)
: Math::Prime::Util::GMP::next_prime($q)
)
) {
__SUB__->($m * $q, $q, $k - 1, $j++);
}
}
Expand All @@ -10713,6 +10718,103 @@ package Sidef::Types::Number::Number {
*pi_k = \&almost_prime_count;
*almost_primepi = \&almost_prime_count;

sub squarefree_almost_prime_count {
my ($k, $from, $to) = @_;

_valid(\$from);

if (defined($to)) {
_valid(\$to);
return ZERO if $to->lt($from);
return $k->squarefree_almost_prime_count($to)->sub($k->squarefree_almost_prime_count($from->dec));
}

$k = _any2ui($$k) // return ZERO;

if ($k == 0) {
return ONE;
}
elsif ($k == 1) {
return $_[1]->prime_count;
}

my $v = _big2uistr($from) // return ZERO;

state $t = Math::GMPz::Rmpz_init_nobless();
state $u = Math::GMPz::Rmpz_init_nobless();
state $v = Math::GMPz::Rmpz_init_nobless();

my $count = Math::GMPz::Rmpz_init_set_ui(0);

my $n = _any2mpz($$from) // return ZERO;

sub {
my ($m, $p, $k, $j) = @_;

my $s = do {
Math::GMPz::Rmpz_div($t, $n, $m);
Math::GMPz::Rmpz_root($t, $t, $k);
Math::GMPz::Rmpz_get_ui($t);
};

if ($k == 2) {

for (
my $q = $p ;
$q <= $s ;
$q = (
HAS_PRIME_UTIL
? Math::Prime::Util::next_prime($q)
: Math::Prime::Util::GMP::next_prime($q)
)
) {

++$j;

if (!Math::GMPz::Rmpz_divisible_ui_p($m, $q)) {

Math::GMPz::Rmpz_mul_ui($t, $m, $q);
Math::GMPz::Rmpz_div($t, $n, $t);

my $pi = _prime_count(Math::GMPz::Rmpz_get_str($t, 10));

if ($pi < ULONG_MAX) {
Math::GMPz::Rmpz_add_ui($count, $count, $pi - $j);
}
else {
Math::GMPz::Rmpz_set_str($t, "$pi", 10);
Math::GMPz::Rmpz_sub_ui($t, $t, $j);
Math::GMPz::Rmpz_add($count, $count, $t);
}
}
}

return;
}

for (
;
$p <= $s ;
$p = (
HAS_PRIME_UTIL
? Math::Prime::Util::next_prime($p)
: Math::Prime::Util::GMP::next_prime($p)
)
) {
if (!Math::GMPz::Rmpz_divisible_ui_p($m, $p)) {
__SUB__->($m * $p, $p, $k - 1, $j);
}
++$j;
}
}
->(Math::GMPz::Rmpz_init_set_ui(1), 2, $k, 0);

bless \$count;
}

*squarefree_pi_k = \&squarefree_almost_prime_count;
*squarefree_almost_primepi = \&squarefree_almost_prime_count;

sub omega_prime_count {
my ($k, $from, $to) = @_;

Expand Down Expand Up @@ -12828,8 +12930,7 @@ package Sidef::Types::Number::Number {
return _set_int("$r");
}

my $r = Math::Prime::Util::GMP::vecsum(Math::Prime::Util::GMP::sieve_primes($x, $y));
_set_int($r);
_set_int(Math::Prime::Util::GMP::vecsum(Math::Prime::Util::GMP::sieve_primes($x, $y)));
}

*prime_sum = \&sum_primes;
Expand Down Expand Up @@ -16146,15 +16247,7 @@ package Sidef::Types::Number::Number {

my $s = Math::Prime::Util::GMP::rootint(Math::Prime::Util::GMP::divint($B, $m), $k);

for (
;
$p <= $s ;
$p = (
HAS_PRIME_UTIL
? Math::Prime::Util::next_prime($p)
: Math::Prime::Util::GMP::next_prime($p)
)
) {
foreach my $p (Math::Prime::Util::GMP::sieve_primes($p, $s)) {

if ($squarefree and Math::GMPz::Rmpz_divisible_ui_p($m, $p)) {
next;
Expand Down
20 changes: 18 additions & 2 deletions lib/Sidef/Types/Number/Number.pod
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4915,8 +4915,8 @@ Returns the number of k-almost primes <= n, or in the range C<a..b>.

Example:

say 5.almost_prime_count(1e7) #=> 1349779
say 5.almost_prime_count(1e6, 1e7) #=> 1225314
say 1.almost_prime_count(100) # count of primes <= 100
say 2.almost_prime_count(50, 100) # count of semiprimes in range 50..100

Aliases: I<almost_primepi>, I<almost_prime_count>

Expand Down Expand Up @@ -6078,6 +6078,22 @@ Example:

=cut

=head2 squarefree_pi_k

k.squarefree_almost_prime_count(n)
k.squarefree_almost_prime_count(a,b)

Returns the number of squarefree k-almost primes <= n, or in the range C<a..b>.

Example:

say 1.squarefree_almost_prime_count(1000) # count of primes <= 1000
say 2.squarefree_almost_prime_count(50, 1000) # count of squarefree semiprimes in range 50..1000

Aliases: I<squarefree_almost_primepi>, I<squarefree_almost_prime_count>

=cut

=head2 squarefree_sigma

squarefree_sigma(n,k=1)
Expand Down
19 changes: 16 additions & 3 deletions scripts/Tests/number_methods.sf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -280,11 +280,13 @@ for k in (1..4) {

assert_eq(k.omega_prime_count(n), 1..n -> count { .is_omega_prime(k) })
assert_eq(k.almost_prime_count(n), 1..n -> count { .is_almost_prime(k) })
assert_eq(k.squarefree_almost_prime_count(n), 1..n -> count { .is_squarefree && .is_almost_prime(k) })

n = irand(100, 1000)

assert_eq(k.omega_prime_count(n), k.omega_primes(n).len)
assert_eq(k.almost_prime_count(n), k.almost_primes(n).len)
assert_eq(k.squarefree_almost_prime_count(n), k.squarefree_almost_primes(n).len)
}

assert_eq(2.almost_primes(50, 100).len, 2.almost_primepi(50, 100))
Expand All @@ -299,9 +301,20 @@ assert_eq(2.omega_primes(50, 106).len, 2.omega_prime_count(50, 106))
assert_eq(3.omega_primes(50, 106).len, 3.omega_prime_count(50, 106))
assert_eq(3.omega_primes(49, 105).len, 3.omega_prime_count(49, 105))

assert_eq(1.omega_primes(49, 105), range(49, 105).grep { .is_omega_prime(1) })
assert_eq(2.omega_primes(49, 105), range(49, 105).grep { .is_omega_prime(2) })
assert_eq(3.omega_primes(49, 105), range(49, 105).grep { .is_omega_prime(3) })
assert_eq(2.squarefree_almost_primes(50, 100).len, 2.squarefree_almost_prime_count(50, 100))
assert_eq(2.squarefree_almost_primes(10, 106).len, 2.squarefree_almost_prime_count(10, 106))
assert_eq(2.squarefree_almost_primes(50, 106).len, 2.squarefree_almost_prime_count(50, 106))
assert_eq(3.squarefree_almost_primes(50, 106).len, 3.squarefree_almost_prime_count(50, 106))
assert_eq(3.squarefree_almost_primes(49, 105).len, 3.squarefree_almost_prime_count(49, 105))

for k in (1..4) {

var a = irand(50, 100)
var b = irand(a, 150)

assert_eq(k.omega_primes(a,b), range(a,b).grep { .is_omega_prime(k) })
assert_eq(k.squarefree_almost_primes(a,b), range(a,b).grep { .is_squarefree && .is_almost_prime(k) })
}

assert_eq(2.powerful(50, 100).len, 2.powerful_count(50, 100))
assert_eq(2.powerful(10, 106).len, 2.powerful_count(10, 106))
Expand Down

0 comments on commit a236a25

Please sign in to comment.