- Added the Number lucas_factor(n, j=1, tries=100). method.

Effective in factoring Carmichael numbers, Fermat pseudoprimes, Lucas pseudoprimes and Lucas-Carmichael numbers. Example: say lucas_factor(2425361208749736840354501506901183117777758034612345610725789878400467) Output: [19760381716819645293083, 151496259828950613913643, 810175650389605457016443] ...taking ~0.01 seconds.
trizen · Aug 27, 2020 · 2d91cf5 · 2d91cf5
1 parent d435c38
commit 2d91cf5
Show file tree

Hide file tree

Showing 2 changed files with 161 additions and 25 deletions.
diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -12490,33 +12490,40 @@ package Sidef::Types::Number::Number {
         );
     }
 
-    sub miller_factor {
+    sub _miller_factor {
         my ($n, $tries) = @_;
 
-        _valid(\$tries) if defined($tries);
-        $n = _any2mpz($$n) // return Sidef::Types::Array::Array->new();
+        # $n is a Math::GMPz object holding a positive value
+        # $tries is a positive native integer or `undef`
+
+        if ($HAS_PRIME_UTIL and Math::GMPz::Rmpz_fits_ulong_p($n)) {
+            return Math::Prime::Util::factor(Math::GMPz::Rmpz_get_ui($n));
+        }
+
+        my $D = Math::GMPz::Rmpz_init();    # n-1
+        Math::GMPz::Rmpz_sub_ui($D, $n, 1);
 
-        my $D = $n - 1;
-        my $s = Math::GMPz::Rmpz_scan1($D, 0);
+        my $s = Math::GMPz::Rmpz_scan1($D, 0) || return ($n);
         my $r = $s - 1;
-        my $d = $D >> $s;
 
-        if (defined($tries)) {
-            $tries = _any2ui($$tries) // 100;
-        }
-        else {
-            $tries = 100;
-            if ($s > 20 and $tries > 10) {
-                $tries = 10;
-            }
-        }
+        my $d = Math::GMPz::Rmpz_init();    # D >> s
+        Math::GMPz::Rmpz_div_2exp($d, $D, $s);
+
+        $tries //= 1 + CORE::int(200 / $s);
 
         my $x = Math::GMPz::Rmpz_init();
         my $g = Math::GMPz::Rmpz_init();
 
+        state $t = Math::GMPz::Rmpz_init_nobless();
+
         for (1 .. $tries) {
 
-            my $p = ($HAS_PRIME_UTIL ? Math::Prime::Util::random_prime(1e7) : Math::Prime::Util::GMP::random_prime(1e7));
+            my $p = (
+                     $HAS_PRIME_UTIL
+                     ? Math::Prime::Util::random_prime(1e7)
+                     : Math::Prime::Util::GMP::random_prime(1e7)
+                    );
+
             Math::GMPz::Rmpz_powm($x, Math::GMPz::Rmpz_init_set_ui($p), $d, $n);
 
             foreach my $k (0 .. $r) {
@@ -12525,28 +12532,147 @@ package Sidef::Types::Number::Number {
                 last if (Math::GMPz::Rmpz_cmp($x, $D) == 0);
 
                 foreach my $i (1, -1) {
-                    Math::GMPz::Rmpz_gcd($g, $x + $i, $n);
-                    if (    Math::GMPz::Rmpz_cmp_ui($g, 1) > 0
-                        and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
+
+                    ($i < 0)
+                      ? Math::GMPz::Rmpz_sub_ui($t, $x, -$i)
+                      : Math::GMPz::Rmpz_add_ui($t, $x, +$i);
+
+                    Math::GMPz::Rmpz_gcd($g, $t, $n);
+
+                    if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0 and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
+
                         Math::GMPz::Rmpz_divexact($x, $n, $g);
 
-                        my @g_factors = Math::GMPz::Rmpz_probab_prime_p($g, 5) ? (bless \$g) : @{__SUB__->(bless \$g)};
-                        my @x_factors = Math::GMPz::Rmpz_probab_prime_p($x, 5) ? (bless \$x) : @{__SUB__->(bless \$x)};
+                        my @g_factors = (Math::GMPz::Rmpz_probab_prime_p($g, 5) ? $g : __SUB__->($g));
+                        my @x_factors = (Math::GMPz::Rmpz_probab_prime_p($x, 5) ? $x : __SUB__->($x));
 
-                        return Sidef::Types::Array::Array->new(
-                                                           [sort { Math::GMPz::Rmpz_cmp($$a, $$b) } (@g_factors, @x_factors)]);
+                        return (@g_factors, @x_factors);
                     }
                 }
 
                 Math::GMPz::Rmpz_powm_ui($x, $x, 2, $n);
             }
         }
 
-        Sidef::Types::Array::Array->new([bless(\$n)]);
+        return ($n);
+    }
+
+    sub miller_factor {
+        my ($n, $tries) = @_;
+
+        $n = _any2mpz($$n) // return Sidef::Types::Array::Array->new();
+        Math::GMPz::Rmpz_sgn($n) > 0 or return Sidef::Types::Array::Array->new();
+
+        if (defined($tries)) {
+            _valid(\$tries);
+            $tries = _any2ui($$tries) // undef;
+        }
+
+        my @factors = map { ref($_) ? $_ : Math::GMPz::Rmpz_init_set_ui($_) } _miller_factor($n, $tries);
+        @factors = sort { Math::GMPz::Rmpz_cmp($a, $b) } @factors;
+        @factors = map  { bless \$_ } @factors;
+
+        Sidef::Types::Array::Array->new(\@factors);
     }
 
     *miller_rabin_factor = \&miller_factor;
 
+    sub _lucas_factor {
+        my ($n, $j, $tries) = @_;
+
+        # $n is a Math::GMPz object holding a positive value
+        # $j is a signed native integer or `undef`
+        # $tries is an unsigned native integer or `undef`
+
+        if ($HAS_PRIME_UTIL and Math::GMPz::Rmpz_fits_ulong_p($n)) {
+            return Math::Prime::Util::factor(Math::GMPz::Rmpz_get_ui($n));
+        }
+
+        $j //= 1;
+
+        my $D = Math::GMPz::Rmpz_init();    # n + j
+
+        ($j < 0)
+          ? Math::GMPz::Rmpz_sub_ui($D, $n, -$j)
+          : Math::GMPz::Rmpz_add_ui($D, $n, +$j);
+
+        my $s = Math::GMPz::Rmpz_scan1($D, 0) || return ($n);
+        my $r = $s;
+
+        my $d = Math::GMPz::Rmpz_init();    # D >> s
+        Math::GMPz::Rmpz_div_2exp($d, $D, $s);
+
+        $tries //= 1 + CORE::int(100 / $s);
+
+        my $x = Math::GMPz::Rmpz_init();
+        my $g = Math::GMPz::Rmpz_init();
+
+        my $PQ_limit = 1e6;
+
+        if (Math::GMPz::Rmpz_cmp_ui($n, 1e6) <= 0) {
+            $PQ_limit = Math::GMPz::Rmpz_get_ui($n) - 1;
+        }
+
+        foreach my $i (1 .. $tries) {
+
+            my $P = 1 + CORE::int(CORE::rand($PQ_limit));
+            my $Q = 1 + CORE::int(CORE::rand($PQ_limit));
+
+            $Q *= -1 if (CORE::rand(1) < 0.5);
+
+            next if Math::Prime::Util::GMP::is_square($P * $P - 4 * $Q);
+
+            Math::GMPz::Rmpz_set($x, $d);
+
+            foreach my $k (0 .. $r) {
+
+                my ($U, $V) = Math::Prime::Util::GMP::lucas_sequence($n, $P, $Q, $x);
+
+                foreach my $t ($U, $V, Math::Prime::Util::GMP::subint($V, $P)) {
+                    Math::GMPz::Rmpz_set_str($g, $t, 10);
+                    Math::GMPz::Rmpz_gcd($g, $g, $n);
+                    if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0 and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
+                        Math::GMPz::Rmpz_divexact($x, $n, $g);
+#<<<
+                        my @g_factors = (Math::GMPz::Rmpz_probab_prime_p($g, 5) ? $g : __SUB__->($g, [-1, +1]->[CORE::int(CORE::rand(2))]));
+                        my @x_factors = (Math::GMPz::Rmpz_probab_prime_p($x, 5) ? $x : __SUB__->($x, [-1, +1]->[CORE::int(CORE::rand(2))]));
+#>>>
+                        return (@g_factors, @x_factors);
+                    }
+                }
+
+                Math::GMPz::Rmpz_mul_2exp($x, $x, 1);
+            }
+        }
+
+        return ($n);
+    }
+
+    sub lucas_factor {
+        my ($n, $j, $tries) = @_;
+
+        $n = _any2mpz($$n) // return Sidef::Types::Array::Array->new();
+        Math::GMPz::Rmpz_sgn($n) > 0 or return Sidef::Types::Array::Array->new();
+
+        if (defined($j)) {
+            _valid(\$j);
+            $j = _any2si($$j) // undef;
+        }
+
+        if (defined($tries)) {
+            _valid(\$tries);
+            $tries = _any2ui($$tries) // undef;
+        }
+
+        my @factors = map { ref($_) ? $_ : Math::GMPz::Rmpz_init_set_ui($_) } _lucas_factor($n, $tries);
+        @factors = sort { Math::GMPz::Rmpz_cmp($a, $b) } @factors;
+        @factors = map  { bless \$_ } @factors;
+
+        Sidef::Types::Array::Array->new(\@factors);
+    }
+
+    *lucas_miller_factor = \&lucas_factor;
+
     sub fermat_factor {
         my ($n, $k) = @_;
         _valid(\$k) if defined($k);
@@ -14519,7 +14645,7 @@ package Sidef::Types::Number::Number {
 
         Math::Prime::Util::GMP::is_carmichael($n) || return Sidef::Types::Bool::Bool::FALSE;
 
-        my @factors = map { $$_ } @{(bless \$n)->miller_factor};
+        my @factors = map { ref($_) ? Math::GMPz::Rmpz_get_str($_, 10) : $_ } _miller_factor($n);
         my @primes;
 
         foreach my $k (@factors) {

diff --git a/lib/Sidef/Types/Number/Number.pod b/lib/Sidef/Types/Number/Number.pod
@@ -3156,6 +3156,16 @@ Return the
 
 =cut
 
+=head2 lucas_factor
+
+Number.lucas_factor() -> I<Obj>
+
+Return the
+
+Aliases: I<lucas_miller_factor>
+
+=cut
+
 =head2 lucas_mod
 
 Number.lucas_mod() -> I<Obj>