diff --git a/lib/Sidef/Math/Math.pm b/lib/Sidef/Math/Math.pm
index 753e214f..000d3acf 100644
--- a/lib/Sidef/Math/Math.pm
+++ b/lib/Sidef/Math/Math.pm
@@ -170,37 +170,7 @@ package Sidef::Math::Math {
 
     sub gcd_factors {
         my ($self, $n, $arr) = @_;
-
-        $n->is_pos
-          || return Sidef::Types::Array::Array->new;
-
-        my $orig_n = $n;
-
-        my @factors;
-        my $G = Sidef::Types::Array::Array->new([grep { $_->is_ntf($n) } map { $n->gcd($_) } @$arr])->sort->uniq;
-
-        foreach my $g (@$G) {
-
-            my $new_g = $g;
-
-            if ($new_g->divides($n)) {
-                ## ok
-            }
-            else {
-                $new_g = $n->gcd($g);
-                $new_g->is_ntf($n) || next;
-            }
-
-            my $v = $n->valuation($new_g);
-            $n = $n->remove($new_g);
-            push @factors, ($new_g) x CORE::int($v);
-        }
-
-        if (!$n->is_one) {
-            push @factors, $n;
-        }
-
-        Sidef::Types::Array::Array->new(\@factors)->sort;
+        Sidef::Types::Number::Number::gcd_factors($n, $arr);
     }
 
     sub smooth_numbers {
diff --git a/lib/Sidef/Math/Math.pod b/lib/Sidef/Math/Math.pod
index 57785f9a..462d3725 100644
--- a/lib/Sidef/Math/Math.pod
+++ b/lib/Sidef/Math/Math.pod
@@ -105,9 +105,9 @@ Returns the greatest common divisor (GCD) of a list of integers.
 
 =head2 gcd_factors
 
-    Math.gcd_factors(n, array)
+    Math.gcd_factors(n, [a, b, ...])
 
-Given a positive integer and a list of numbers, it tries to find non-trivial factors of n, checking each C<gcd(n, array[0])>, C<gcd(n, array[1])>, etc.
+Given a positive integer and an array of integers, it tries to find non-trivial factors of n, checking each C<gcd(n, array[0])>, C<gcd(n, array[1])>, etc.
 
     var n = 43*43*97*503
     var a = [19*43*97, 1, 13*41*43*101]
diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
index f928bb38..5c0e70d4 100644
--- a/lib/Sidef/Types/Number/Number.pm
+++ b/lib/Sidef/Types/Number/Number.pm
@@ -16266,6 +16266,55 @@ package Sidef::Types::Number::Number {
         _set_int($f[-1]);
     }
 
+    sub gcd_factors {
+        my ($n, $arr) = @_;
+
+        $n = _any2mpz($$n) // return Sidef::Types::Array::Array->new;
+
+        Math::GMPz::Rmpz_sgn($n) > 0
+          or return Sidef::Types::Array::Array->new;
+
+        my $z = Math::GMPz::Rmpz_init_set($n);    # copy
+        state $t = Math::GMPz::Rmpz_init_nobless();
+
+        my @gcds;
+        my %seen;
+
+        foreach my $k (@$arr) {
+            _valid(\$k);
+            my $m = _any2mpz($$k) // next;
+            Math::GMPz::Rmpz_gcd($t, $z, $m);
+            Math::GMPz::Rmpz_cmp_ui($t, 1) > 0 or next;
+            Math::GMPz::Rmpz_cmp($t, $z) < 0   or next;
+            if (!$seen{Math::GMPz::Rmpz_get_str($t, 10)}++) {
+                push @gcds, Math::GMPz::Rmpz_init_set($t);
+            }
+        }
+
+        @gcds = sort { Math::GMPz::Rmpz_cmp($a, $b) } @gcds;
+
+        my @factors;
+
+        foreach my $g (@gcds) {
+
+            Math::GMPz::Rmpz_gcd($t, $g, $z);
+            Math::GMPz::Rmpz_cmp_ui($t, 1) > 0 or next;
+            Math::GMPz::Rmpz_cmp($t, $z) < 0   or next;
+
+            my $v = Math::GMPz::Rmpz_remove($z, $z, $t);
+            push(@factors, (Math::GMPz::Rmpz_init_set($t)) x $v);
+        }
+
+        if (Math::GMPz::Rmpz_cmp_ui($z, 1) > 0) {
+            push @factors, $z;
+        }
+
+        @factors = sort { Math::GMPz::Rmpz_cmp($a, $b) } @factors;
+        @factors = map  { bless \$_ } @factors;
+
+        Sidef::Types::Array::Array->new(\@factors);
+    }
+
     sub factor {
         my ($n, $block) = @_;
 
@@ -16282,7 +16331,7 @@ package Sidef::Types::Number::Number {
                     }
                 )
             )->uniq;
-            return Sidef::Math::Math->gcd_factors($n, $f);
+            return $n->gcd_factors($f);
         }
 
         $n = _big2pistr($n) // return Sidef::Types::Array::Array->new();
diff --git a/lib/Sidef/Types/Number/Number.pod b/lib/Sidef/Types/Number/Number.pod
index f52fe04c..1dda1a42 100644
--- a/lib/Sidef/Types/Number/Number.pod
+++ b/lib/Sidef/Types/Number/Number.pod
@@ -2751,6 +2751,20 @@ The value of C<d> is always non-negative.
 
 =cut
 
+=head2 gcd_factors
+
+    gcd_factors(n, [a, b, ...])
+
+Given a positive integer and an array of integers, it tries to find non-trivial factors of n, checking each C<gcd(n, array[0])>, C<gcd(n, array[1])>, etc.
+
+    var n = 43*43*97*503
+    var a = [19*43*97, 1, 13*41*43*101]
+    say gcd_factors(n, a)                #=> [43, 43, 97, 503]
+
+The product of the factors gives back C<n>. However, some factors may be composite.
+
+=cut
+
 =head2 gcud
 
     gcud(...)