- Faster factorization of large integers with many small prime factors.

Example: say factor(lcm(1..38861)).len #=> 4175 The above code now takes only ~0.08 seconds. Before it took ~8 seconds. Related to: danaj/Math-Prime-Util-GMP#27
trizen · Feb 1, 2022 · 4b9567b · 4b9567b
1 parent 827ae37
commit 4b9567b
Show file tree

Hide file tree

Showing 2 changed files with 88 additions and 45 deletions.
diff --git a/lib/Sidef.pm b/lib/Sidef.pm
@@ -88,6 +88,7 @@ package Sidef {
 
     sub get_sidef_config_dir {
         my ($self) = @_;
+
         $self->{sidef_config_dir} //= $ENV{SIDEF_CONFIG_DIR}
           || File::Spec->catdir(
                                 $ENV{XDG_CONFIG_DIR}
@@ -107,8 +108,8 @@ package Sidef {
 
         if (not -d $self->{sidef_config_dir}) {
             require File::Path;
-            File::Path::make_path($self->{sidef_config_dir})
-                or warn "[WARNING] Can't create directory <<$self->{sidef_config_dir}>>: $!";
+            eval { File::Path::make_path($self->{sidef_config_dir}) }
+              or warn "[WARNING] Can't create directory <<$self->{sidef_config_dir}>>: $!";
         }
 
         return $self->{sidef_config_dir};

diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -818,6 +818,34 @@ package Sidef::Types::Number::Number {
         Math::GMPz::Rmpz_get_str($x, 10);
     }
 
+    sub _factor {
+        my ($n) = @_;
+
+        if (ref($n) eq 'Math::GMPz') {
+            if (HAS_PRIME_UTIL and Math::GMPz::Rmpz_fits_ulong_p($n)) {
+                $n = Math::GMPz::Rmpz_get_ui($n);
+            }
+            else {
+                $n = Math::GMPz::Rmpz_get_str($n, 10);
+            }
+        }
+
+        my @factors;
+
+        if (length($n) > 2000) {
+            ($n, @factors) = _adaptive_trial_factor($n);
+        }
+
+        (
+         @factors,
+         (
+          (HAS_PRIME_UTIL and $n < ULONG_MAX)
+          ? Math::Prime::Util::factor($n)
+          : Math::Prime::Util::GMP::factor($n)
+         )
+        );
+    }
+
     # Prime factorization in [p,k] form, where k is the multiplicity of p.
     sub _factor_exp {
         my ($n) = @_;
@@ -835,7 +863,7 @@ package Sidef::Types::Number::Number {
             return Math::Prime::Util::factor_exp($n);
         }
 
-        my @factors = Math::Prime::Util::GMP::factor($n);
+        my @factors = _factor($n);
 
         my $prev_value = shift(@factors) // return;
         my @factor_exp = [$prev_value, 1];
@@ -853,23 +881,6 @@ package Sidef::Types::Number::Number {
         @factor_exp;
     }
 
-    sub _factor {
-        my ($n) = @_;
-
-        if (ref($n) eq 'Math::GMPz') {
-            if (HAS_PRIME_UTIL and Math::GMPz::Rmpz_fits_ulong_p($n)) {
-                $n = Math::GMPz::Rmpz_get_ui($n);
-            }
-            else {
-                $n = Math::GMPz::Rmpz_get_str($n, 10);
-            }
-        }
-
-        (HAS_PRIME_UTIL and $n < ULONG_MAX)
-          ? Math::Prime::Util::factor($n)
-          : Math::Prime::Util::GMP::factor($n);
-    }
-
     sub _divisors {
         my ($n) = @_;
 
@@ -1009,25 +1020,66 @@ package Sidef::Types::Number::Number {
     }
 
     sub _adaptive_trial_factor {
-        my ($n) = @_;
+        my ($n, $L, $R) = @_;
 
-        my @limits = (5e4);
-        my $size   = Math::GMPz::Rmpz_sizeinbase($n, 10);
+        $L //= 5e4;
+        $R //= 5e6;
 
-        if ($size > 55) { push @limits, 5e5 }
-        if ($size > 65) { push @limits, 5e6 }
+        if (ref($n) eq 'Math::GMPz') {
+            $n = Math::GMPz::Rmpz_init_set($n);    # copy
+        }
+        else {
+            $n = Math::GMPz::Rmpz_init_set_str($n, 10);
+        }
+
+        my @factors;
+
+        my $P = _cached_primorial($L);
+
+        my $g = Math::GMPz::Rmpz_init();
+        my $t = Math::GMPz::Rmpz_init();
+
+        my $F = 2;
+
+        while (1) {
+
+            Math::GMPz::Rmpz_gcd($g, $P, $n);
+
+            # Early stop when n seems to no longer have small factors
+            if (Math::GMPz::Rmpz_cmp_ui($g, 1) == 0) {
+                last;
+            }
 
-        my @prime_factors;
+            # Factorize n over primes in P
+            foreach my $p (Math::Prime::Util::GMP::sieve_primes($F, $L)) {
+                if (Math::GMPz::Rmpz_divisible_ui_p($g, $p)) {
 
-        foreach my $trial_limit (@limits) {
-            my ($r, @factors) = _primorial_trial_factor($n, $trial_limit);
-            @factors || last;
-            push @prime_factors, @factors;
-            $n = $r;
-            last if Math::GMPz::Rmpz_fits_ulong_p($n);
+                    Math::GMPz::Rmpz_set_ui($t, $p);
+                    my $valuation = Math::GMPz::Rmpz_remove($n, $n, $t);
+                    push @factors, ($p) x $valuation;
+
+                    # Stop the loop early when no more primes divide `u` (optional)
+                    Math::GMPz::Rmpz_divexact_ui($g, $g, $p);
+                    last if (Math::GMPz::Rmpz_cmp_ui($g, 1) == 0);
+                }
+            }
+
+            # Early stop when n has been fully factored
+            if (Math::GMPz::Rmpz_cmp_ui($n, 1) == 0) {
+                last;
+            }
+
+            # Early stop when the trial range has been exhausted
+            if ($L > $R) {
+                last;
+            }
+
+            $F = $L;
+            $L <<= 1;
+            $P = _cached_primorial($L);
         }
 
-        return ($n, @prime_factors);
+        return ($n, @factors);
     }
 
     #
@@ -14782,7 +14834,7 @@ package Sidef::Types::Number::Number {
             last if (($j >= 7) && ($size <= 150));    # 45 digits
         }
 
-        my @f = Math::Prime::Util::GMP::factor(Math::GMPz::Rmpz_get_str($n, 10));
+        my @f = _factor($n);
         _set_int($f[0]);
     }
 
@@ -14801,18 +14853,8 @@ package Sidef::Types::Number::Number {
     }
 
     sub factor {
-
         my $n = &_big2pistr // return Sidef::Types::Array::Array->new();
-        my @factors;
-
-        if (HAS_PRIME_UTIL and $n < ULONG_MAX) {
-            @factors = Math::Prime::Util::factor($n);
-        }
-        else {
-            @factors = Math::Prime::Util::GMP::factor($n);
-        }
-
-        Sidef::Types::Array::Array->new([map { _set_int($_) } @factors]);
+        Sidef::Types::Array::Array->new([map { _set_int($_) } _factor($n)]);
     }
 
     *factors = \&factor;