- Added the Number n.omega_prime_divisors(k) method.

Returns the k-omega prime divisors of n. Example: say 5040.omega_prime_divisors(4) #=> [210, 420, 630, 840, 1260, 1680, 2520, 5040] When `k` is omitted, an array of arrays with the k-omega prime divisors of n, for each k in the range 0..omega(n), is returned: say 120.omega_prime_divisors #=> [[1], [2, 3, 4, 5, 8], [6, 10, 12, 15, 20, 24, 40], [30, 60, 120]]
trizen · Mar 21, 2021 · 6737ed3 · 6737ed3
1 parent ef48f44
commit 6737ed3
Show file tree

Hide file tree

Showing 3 changed files with 112 additions and 54 deletions.
diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -10953,6 +10953,7 @@ package Sidef::Types::Number::Number {
                         __SUB__->($w, $p, $k - 1, $j, $s);
                     }
                 }
+
                 ++$j;
             }
           }
@@ -16018,6 +16019,85 @@ package Sidef::Types::Number::Number {
         return \@omega_primes;
     }
 
+    sub omega_prime_divisors {
+        my ($n, $k) = @_;
+
+        if (!defined($k)) {
+            my @list;
+
+            foreach my $k (0 .. $n->omega) {
+                push @list, $n->omega_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]);
+        }
+
+        if ($k == 1) {
+            return $_[0]->prime_power_divisors;
+        }
+
+        my @factor_exp = _factor_exp(Math::GMPz::Rmpz_get_str($z, 10));
+
+        if ($k > scalar(@factor_exp)) {
+            return Sidef::Types::Array::Array->new();
+        }
+
+        my %valuations = map { @$_ } @factor_exp;
+        my @factors    = map { ($_ < ULONG_MAX) ? $_ : ${_set_int($_)} } map { $_->[0] } @factor_exp;
+
+        my $t = Math::GMPz::Rmpz_init();
+
+        my @list;
+
+        sub {
+            my ($m, $p, $k) = @_;
+
+            Math::GMPz::Rmpz_div($t, $z, $m);
+            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;
+
+                if (ref($q) ? Math::GMPz::Rmpz_divisible_p($m, $q) : Math::GMPz::Rmpz_divisible_ui_p($m, $q)) {
+                    next;
+                }
+
+                my $v = $m * $q;
+
+                foreach my $j (1 .. $valuations{$q}) {
+
+                    if ($k == 1) {
+                        push @list, $v;
+                    }
+                    else {
+                        __SUB__->($v, $q, $k - 1);
+                    }
+
+                    $v *= $q;
+                }
+            }
+          }
+          ->(Math::GMPz::Rmpz_init_set_ui(1), $factors[0], $k);
+
+        Sidef::Types::Array::Array->new([map { _set_int($_) } sort { $a <=> $b } @list]);
+    }
+
     sub omega_primes {
         my ($k, $from, $to) = @_;
 
@@ -16306,60 +16386,6 @@ package Sidef::Types::Number::Number {
             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;
@@ -16368,8 +16394,20 @@ package Sidef::Types::Number::Number {
             return Sidef::Types::Array::Array->new([map { _set_int($_) } @factors]);
         }
 
+        my $bigomega = 0;
+
+        foreach my $pp (@factor_exp) {
+            $bigomega += $pp->[1];
+        }
+
+        if ($k > $bigomega) {
+            return Sidef::Types::Array::Array->new();
+        }
+
         my $t = Math::GMPz::Rmpz_init();
 
+        my @list;
+
         sub {
             my ($m, $p, $k) = @_;
 

diff --git a/lib/Sidef/Types/Number/Number.pod b/lib/Sidef/Types/Number/Number.pod
@@ -4815,6 +4815,23 @@ Example:
 
 =cut
 
+=head2 omega_prime_divisors
+
+    n.omega_prime_divisors
+    n.omega_prime_divisors(k)
+
+Returns the k-omega prime divisors of n.
+
+Example:
+
+    say 5040.omega_prime_divisors(4)    #=> [210, 420, 630, 840, 1260, 1680, 2520, 5040]
+
+When C<k> is omitted, an array of arrays with the k-omega prime divisors of n, for each k in the range C<0..omega(n)>, is returned:
+
+    say 120.omega_prime_divisors        #=> [[1], [2, 3, 4, 5, 8], [6, 10, 12, 15, 20, 24, 40], [30, 60, 120]]
+
+=cut
+
 =head2 omega_primes
 
     k.omega_primes(n)

diff --git a/scripts/Tests/number_methods.sf b/scripts/Tests/number_methods.sf
@@ -339,11 +339,13 @@ for k in (1..7) {
 # Almost prime divisors
 
 assert_eq(5040.almost_prime_divisors.flat.sort, 5040.divisors)
+assert_eq(5040.omega_prime_divisors.flat.sort, 5040.divisors)
 
 with (10!) {|n|
     var divisors = n.divisors
     for k in (0 .. (1+n.bigomega)) {
         assert_eq(n.almost_prime_divisors(k), divisors.grep { .is_almost_prime(k) })
+        assert_eq(n.omega_prime_divisors(k), divisors.grep { .is_omega_prime(k) })
     }
 }
 
@@ -354,6 +356,7 @@ with (21!) {|n|
     assert_eq(almost_prime_divisors(n, 37), n.divisors.grep { .is_almost_prime(37) })
     assert_eq(almost_prime_divisors(n, n.bigomega), [n])
     assert_eq(almost_prime_divisors(n, n.bigomega+1), [])
+    assert_eq(omega_prime_divisors(n, n.omega+1), [])
 }
 
 do {