Skip to content

Commit

Permalink
- Added the Number n.almost_prime_divisors(k) method.
Browse files Browse the repository at this point in the history
Returns the k-almost prime divisors of n.

Example:

    say 5040.almost_prime_divisors(7)   #=> [720, 1008, 1680, 2520]

When k is omitted, an array of arrays with the k-almost prime divisors of n, for each k in the range 0..bigomega(n), is returned:

    say 120.almost_prime_divisors       #=> [[1], [2, 3, 5], [4, 6, 10, 15], [8, 12, 20, 30], [24, 40, 60], [120]]
  • Loading branch information
trizen committed Mar 10, 2021
1 parent 64c932e commit a8cc27e
Show file tree
Hide file tree
Showing 2 changed files with 160 additions and 0 deletions.
143 changes: 143 additions & 0 deletions lib/Sidef/Types/Number/Number.pm
Original file line number Diff line number Diff line change
Expand Up @@ -15622,6 +15622,149 @@ package Sidef::Types::Number::Number {
return \@almost_primes;
}

sub almost_prime_divisors {
my ($n, $k) = @_;

if (!defined($k)) {
my @list;

foreach my $k (0 .. $n->bigomega) {
push @list, $n->almost_prime_divisors(_set_int($k));
}

return Sidef::Types::Array::Array->new(\@list);
}

_valid(\$k);

$k = _any2ui($$k) // return Sidef::Types::Array::Array->new;
my $z = _any2mpz($$n) // return Sidef::Types::Array::Array->new;

if ($k == 0) {
return Sidef::Types::Array::Array->new($ONE);
}

my @list;

if (HAS_PRIME_UTIL and Math::GMPz::Rmpz_fits_ulong_p($z)) {

$z = Math::GMPz::Rmpz_get_ui($z);

my @factor_exp = Math::Prime::Util::factor_exp($z);
my %valuations = map { @$_ } @factor_exp;
my @factors = map { $_->[0] } @factor_exp;

if ($k == 1) {
return Sidef::Types::Array::Array->new([map { _set_int($_) } @factors]);
}

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

my $L =
defined(&Math::Prime::Util::divint)
? Math::Prime::Util::divint($z, $m)
: Math::Prime::Util::GMP::divint($z, $m);

if ($k == 1) {

foreach my $q (@factors) {

$q < $p and next;
$q > $L and last;

Math::Prime::Util::valuation($m, $q) < $valuations{$q} or next;

push @list, $m * $q;
}

return;
}

my $s = Math::Prime::Util::rootint($L, $k);

foreach my $q (@factors) {

$q < $p and next;
$q > $s and last;

Math::Prime::Util::valuation($m, $q) < $valuations{$q} or next;

__SUB__->($m * $q, $q, $k - 1);
}
}
->(1, $factors[0], $k);

return Sidef::Types::Array::Array->new([map { _set_int($_) } sort { $a <=> $b } @list]);
}

my @factor_exp = _factor_exp(Math::GMPz::Rmpz_get_str($z, 10));
my %valuations = map { @$_ } @factor_exp;
my @factors = map { ($_ < ULONG_MAX) ? $_ : ${_set_int($_)} } map { $_->[0] } @factor_exp;

if ($k == 1) {
return Sidef::Types::Array::Array->new([map { _set_int($_) } @factors]);
}

my $t = Math::GMPz::Rmpz_init();

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

Math::GMPz::Rmpz_div($t, $z, $m);

if ($k == 1) {

my $L =
Math::GMPz::Rmpz_fits_ulong_p($t)
? Math::GMPz::Rmpz_get_ui($t)
: Math::GMPz::Rmpz_init_set($t);

foreach my $q (@factors) {

$q < $p and next;
$q > $L and last;

my $v = ref($q) ? Math::GMPz::Rmpz_remove($t, $m, $q) : do {
Math::GMPz::Rmpz_set_ui($t, $q);
Math::GMPz::Rmpz_remove($t, $m, $t);
};

$v < $valuations{$q} or next;

push @list, $m * $q;
}

return;
}

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

my $s =
Math::GMPz::Rmpz_fits_ulong_p($t)
? Math::GMPz::Rmpz_get_ui($t)
: Math::GMPz::Rmpz_init_set($t);

foreach my $q (@factors) {

$q < $p and next;
$q > $s and last;

my $v = ref($q) ? Math::GMPz::Rmpz_remove($t, $m, $q) : do {
Math::GMPz::Rmpz_set_ui($t, $q);
Math::GMPz::Rmpz_remove($t, $m, $t);
};

$v < $valuations{$q} or next;

__SUB__->($m * $q, $q, $k - 1);
}
}
->(Math::GMPz::Rmpz_init_set_ui(1), $factors[0], $k);

Sidef::Types::Array::Array->new([map { _set_int($_) } sort { $a <=> $b } @list]);
}

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

Expand Down
17 changes: 17 additions & 0 deletions lib/Sidef/Types/Number/Number.pod
Original file line number Diff line number Diff line change
Expand Up @@ -663,6 +663,23 @@ This is done by first running a primality pre-test on all values and returning e

=cut

=head2 almost_prime_divisors

n.almost_prime_divisors
n.almost_prime_divisors(k)

Returns the k-almost prime divisors of n.

Example:

say 5040.almost_prime_divisors(7) #=> [720, 1008, 1680, 2520]

When C<k> is omitted, an array of arrays with the k-almost prime divisors of n, for each k in the range C<0..bigomega(n)>, is returned:

say 120.almost_prime_divisors #=> [[1], [2, 3, 5], [4, 6, 10, 15], [8, 12, 20, 30], [24, 40, 60], [120]]

=cut

=head2 almost_primes

k.almost_primes(n)
Expand Down

0 comments on commit a8cc27e

Please sign in to comment.