- Optimized Number is_omega_prime(n,k) for much better performance …

…for large n and k.
trizen · Jul 10, 2022 · ac70392 · ac70392
1 parent a77b68a
commit ac70392
Show file tree

Hide file tree

Showing 3 changed files with 112 additions and 6 deletions.
diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -14965,9 +14965,7 @@ package Sidef::Types::Number::Number {
         Math::GMPz::Rmpz_sgn($n) > 0
           or return Sidef::Types::Bool::Bool::FALSE;
 
-        my $bigomega  = 0;
-        my $remainder = $n;
-        my $size      = Math::GMPz::Rmpz_sizeinbase($n, 2);
+        my $size = Math::GMPz::Rmpz_sizeinbase($n, 2);
 
         if ($k >= $size) {    # the smallest k-almost prime is 2^k
             return Sidef::Types::Bool::Bool::FALSE;
@@ -14981,6 +14979,9 @@ package Sidef::Types::Number::Number {
             );
         }
 
+        my $bigomega  = 0;
+        my $remainder = $n;
+
         if ($size > 100) {    # greater than 10^30
             foreach my $j (2 .. 8) {
 
@@ -15055,7 +15056,7 @@ package Sidef::Types::Number::Number {
                     }
 
                     my $prod = Sidef::Types::Number::Number::prod(@prime_factors);
-                    my $c    = (bless \$n)->div($prod);
+                    my $c    = (bless \$n)->idiv($prod);
 
                     $remainder = $$c;
                     $bigomega  = scalar(@prime_factors);
@@ -15128,9 +15129,103 @@ package Sidef::Types::Number::Number {
             );
         }
 
-        my @factors = _factor_exp($n);
+        my $omega     = 0;
+        my $remainder = $n;
+        my $size      = Math::GMPz::Rmpz_sizeinbase($n, 2);
+
+        if ($size > 100) {    # greater than 10^30
+            foreach my $j (2 .. 8) {
+
+                my ($r, @trial_factors) = _primorial_trial_factor($n, 10**$j);
+
+                $omega     = scalar(List::Util::uniq(@trial_factors));
+                $remainder = $r;
+
+                if (Math::GMPz::Rmpz_cmp_ui($r, 1) == 0) {
+                    if ($omega == $k) {
+                        return Sidef::Types::Bool::Bool::TRUE;
+                    }
+                    else {
+                        return Sidef::Types::Bool::Bool::FALSE;
+                    }
+                }
+
+                my $log = __ilog__($r, _next_prime(10**$j));
+
+                $omega + $log + 1 >= $k
+                  or return Sidef::Types::Bool::Bool::FALSE;
+
+                my $r_is_prime_power = Math::Prime::Util::GMP::is_prime_power(Math::GMPz::Rmpz_get_str($r, 10));
+
+                if ($r_is_prime_power) {
+                    if ($omega + 1 == $k) {
+                        return Sidef::Types::Bool::Bool::TRUE;
+                    }
+                    else {
+                        return Sidef::Types::Bool::Bool::FALSE;
+                    }
+                }
+
+                if ($omega + ($r_is_prime_power ? 1 : 2) > $k) {
+                    return Sidef::Types::Bool::Bool::FALSE;
+                }
+
+                last if (($j >= 5) && (Math::GMPz::Rmpz_sizeinbase($r, 2) <= 100));    # 30 digits
+                last if (($j >= 6) && (Math::GMPz::Rmpz_sizeinbase($r, 2) <= 133));    # 40 digits
+                last if (($j >= 7) && (Math::GMPz::Rmpz_sizeinbase($r, 2) <= 150));    # 45 digits
+
+                # Try to find special factors
+                if ($j == 8) {
+                    my @special_factors = @{(bless \$r)->special_factors};
+                    my @gcd_factors     = @{
+                        (bless \$n)->gcd_factors(
+                                     Sidef::Types::Array::Array->new([@special_factors, (map { _set_int($_) } @trial_factors)])
+                        )
+                    };
+
+                    my @prime_factors;
+
+                    foreach my $f (@gcd_factors) {
+                        if (_is_prob_prime($$f)) {
+                            push @prime_factors, $f;
+                        }
+                        elsif (Math::GMPz::Rmpz_sizeinbase($$f, 2) <= 150) {
+                            push @prime_factors, (map { _set_int($_) } _factor($$f));
+                        }
+                    }
+
+                    my $prod = Sidef::Types::Number::Number::prod(@prime_factors);
+                    my $c    = (bless \$n)->idiv($prod);
+
+                    $remainder = $$c;
+                    $omega     = scalar(List::Util::uniq(map { Math::GMPz::Rmpz_get_str($$_, 10) } @prime_factors));
+
+                    if (Math::GMPz::Rmpz_cmp_ui($remainder, 1) == 0) {
+                        if ($omega == $k) {
+                            return Sidef::Types::Bool::Bool::TRUE;
+                        }
+                        else {
+                            return Sidef::Types::Bool::Bool::FALSE;
+                        }
+                    }
+
+                    my $log = __ilog__($remainder, _next_prime(10**$j));
+
+                    $omega + $log + 1 >= $k
+                      or return Sidef::Types::Bool::Bool::FALSE;
+                }
+            }
+        }
+
+        if ($omega > $k) {
+            return Sidef::Types::Bool::Bool::FALSE;
+        }
+
+        my @factors = _factor_exp($remainder);
+
+        $omega += scalar(@factors);
 
-        (scalar(@factors) == $k)
+        ($omega == $k)
           ? Sidef::Types::Bool::Bool::TRUE
           : Sidef::Types::Bool::Bool::FALSE;
     }

diff --git a/scripts/Extended tests/is_almost_prime.sf b/scripts/Extended tests/is_almost_prime.sf
@@ -15,6 +15,8 @@ assert_eq(
     %n[2, 3, 3, 43, 7, 41, 23, 643, 17, 557, 251],
 )
 
+#assert_eq(A242786(21), 1151)
+
 func A241793(n) {
     for k in (1..Inf) {
         var b = bigomega(k)*n
@@ -29,6 +31,8 @@ assert_eq(
     %n[3, 34, 5, 15, 17, 55, 79, 5, 53, 23, 337, 13, 601, 79, 241, 41],
 )
 
+assert_eq(A241793(24), 79)
+
 func A280005(n) {
     for(var p = 2; true; p.next_prime!) {
         var v = (p**n + 1)

diff --git a/scripts/Extended tests/is_omega_prime.sf b/scripts/Extended tests/is_omega_prime.sf
@@ -15,4 +15,11 @@ assert_eq(
     %n[3, 4, 7, 8, 16, 11, 79, 44, 81, 91, 1024, 47],
 )
 
+assert_eq(A219019(14), 389)
+assert_eq(A219019(15), 256)
+assert_eq(A219019(16), 413)
+assert_eq(A219019(18), 373)
+assert_eq(A219019(20), 1000)
+#assert_eq(A219019(21), 4096)
+
 say "** Test passed!"