- Added the Number cyclotomic_factor(n, bases...) method.

Factorization method, based on cyclotomic polynomials, trying to find algebraic factors of `n`. say cyclotomic_factor(2**120 + 1) say cyclotomic_factor((10**258 - 1)/9 - 10**(258/2) - 1, 10) Output: [257, 65281, 4278255361, 18518800563924107521] [10, 11, 11, 91, 101, 10001, 100000001, 10000000000000001, 100000000000000000000000000000001, 909090909090909090909090909090909090909091, 10000000000000000000000000000000000000000000000000000000000000001, 1098901098901098901098901098901098901098900989010989010989010989010989010989010989011] An optional list of bases can be given to restrict the search only to the given bases.
trizen · May 23, 2022 · 39ba86f · 39ba86f
1 parent 9c4c74f
commit 39ba86f
Show file tree

Hide file tree

Showing 3 changed files with 113 additions and 0 deletions.
diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -11445,6 +11445,99 @@ package Sidef::Types::Number::Number {
         bless \$prod;
     }
 
+    sub cyclotomic_factor {
+        my ($n, @bases) = @_;
+
+        $n = _any2mpz($$n) // return Sidef::Types::Array::Array->new;
+
+        Math::GMPz::Rmpz_cmp_ui($n, 1) > 0
+          or return Sidef::Types::Array::Array->new;
+
+        if (@bases) {
+            _valid(\(@bases));
+            @bases = grep { defined($_) } map { _any2mpz($$_) } @bases;
+        }
+        else {
+            @bases = map { ${_set_int($_)} } (2 .. __ilog__($n, 2));
+        }
+
+        my $cyclotomicmod = sub {
+            my ($n, $x, $m) = @_;
+
+            my @factor_exp = _factor_exp($n);
+
+            # Generate the squarefree divisors of n, along
+            # with the number of prime factors of each divisor
+            my @sd;
+            foreach my $pe (@factor_exp) {
+                my ($p) = @$pe;
+
+                $p =
+                  ($p < ULONG_MAX)
+                  ? Math::GMPz::Rmpz_init_set_ui($p)
+                  : Math::GMPz::Rmpz_init_set_str("$p", 10);
+
+                push @sd, map { [$_->[0] * $p, $_->[1] + 1] } @sd;
+                push @sd, [$p, 1];
+            }
+
+            push @sd, [$ONE, 0];
+
+            my $prod = Math::GMPz::Rmpz_init_set_ui(1);
+
+            foreach my $pair (@sd) {
+                my ($d, $c) = @$pair;
+
+                my $base = Math::GMPz::Rmpz_init();
+                my $exp  = CORE::int($n / $d);
+
+                Math::GMPz::Rmpz_powm_ui($base, $x, $exp, $m);    # x^(n/d) mod m
+                Math::GMPz::Rmpz_sub_ui($base, $base, 1);
+
+                if ($c % 2 == 1) {
+                    Math::GMPz::Rmpz_invert($base, $base, $m) || return $base;
+                }
+
+                Math::GMPz::Rmpz_mul($prod, $prod, $base);
+                Math::GMPz::Rmpz_mod($prod, $prod, $m);
+            }
+
+            $prod;
+        };
+
+        $n = Math::GMPz::Rmpz_init_set($n);    # copy
+
+        my @factors;
+        state $g = Math::GMPz::Rmpz_init_nobless();
+
+      OUTER: foreach my $x (@bases) {
+            my $limit = 1 + __ilog__($n, $x);
+
+            foreach my $k (3 .. $limit) {
+                my $c = $cyclotomicmod->($k, $x, $n);
+
+                Math::GMPz::Rmpz_gcd($g, $n, $c);
+                if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0 and Math::GMPz::Rmpz_cmp($g, $n) < 0) {
+
+                    my $valuation = Math::GMPz::Rmpz_remove($n, $n, $g);
+                    push(@factors, (Math::GMPz::Rmpz_init_set($g)) x $valuation);
+
+                    if (Math::GMPz::Rmpz_cmp_ui($n, 1) == 0 or _is_prob_prime($n)) {
+                        last OUTER;
+                    }
+                }
+            }
+        }
+
+        if (Math::GMPz::Rmpz_cmp_ui($n, 1) > 0) {
+            push @factors, $n;
+        }
+
+        @factors = sort { Math::GMPz::Rmpz_cmp($a, $b) } @factors;
+        @factors = map  { bless \$_ } @factors;
+        Sidef::Types::Array::Array->new(\@factors);
+    }
+
     sub powerfree_sum {
         my ($k, $from, $to) = @_;
 

diff --git a/lib/Sidef/Types/Number/Number.pod b/lib/Sidef/Types/Number/Number.pod
@@ -1793,6 +1793,22 @@ Aliases: I<cyclotomic_polynomial>
 
 =cut
 
+=head2 cyclotomic_factor
+
+    cyclotomic_factor(n)
+    cyclotomic_factor(n, bases...)
+
+Factorization method, based on cyclotomic polynomials, trying to find algebraic factors of C<n>.
+
+    say cyclotomic_factor(2**120 + 1)
+    say cyclotomic_factor((10**258 - 1)/9 - 10**(258/2) - 1, 10)
+
+An optional list of bases can be given to restrict the search only to the given bases.
+
+The product of the factors, will give back C<n>. However, some factors may be composite.
+
+=cut
+
 =head2 cyclotomicmod
 
     cyclotomicmod(n, x, m)

diff --git a/scripts/Tests/number_methods.sf b/scripts/Tests/number_methods.sf
@@ -400,6 +400,7 @@ do {
     {|n| assert_eq(squfof_factor(n).prod, n||1) }.each(0..30)
     {|n| assert_eq(qs_factor(n).prod, n||1) }.each(0..30)
     {|n| assert_eq(trial_factor(n).prod, n||1) }.each(0..30)
+    {|n| assert_eq(cyclotomic_factor(n).prod, n||1) }.each(0..30)
 }
 
 do {
@@ -408,4 +409,7 @@ do {
     assert_eq(20.of {.lnhyperfactorial} ~Z=~= 20.of{.hyperfactorial.ln} -> uniq, [true])
 }
 
+assert_eq(cyclotomic_factor(((10**258 - 1)/9 - 10**(258/2) - 1)),     %n[2, 5, 7, 11, 11, 13, 17, 73, 101, 137, 353, 449, 641, 1409, 69857, 5882353, 100000000000000000000000000000001, 909090909090909090909090909090909090909091, 10000000000000000000000000000000000000000000000000000000000000001, 1098901098901098901098901098901098901098900989010989010989010989010989010989010989011])
+assert_eq(cyclotomic_factor(((10**258 - 1)/9 - 10**(258/2) - 1), 10), %n[10, 11, 11, 91, 101, 10001, 100000001, 10000000000000001, 100000000000000000000000000000001, 909090909090909090909090909090909090909091, 10000000000000000000000000000000000000000000000000000000000000001, 1098901098901098901098901098901098901098900989010989010989010989010989010989010989011])
+
 say "** Tests passed!"