- Added the Number gcud(...) method.

Returns the greatest common unitary divisor of a list of integers. Example: say gcud(9!, 5040) #=> 35 say gcud(9!, 5040, 120) #=> 5 - Also added the Array `.gcud` and `.gcud_by { ... }` methods. say [9!, 5040, 120].gcud #=> 5 say [9!, 5040, 120].gcud_by { _**2 } #=> 25
trizen · Nov 26, 2021 · 2e1635e · 2e1635e
1 parent a59e5e3
commit 2e1635e
Show file tree

Hide file tree

Showing 5 changed files with 77 additions and 37 deletions.
diff --git a/lib/Sidef/Types/Array/Array.pm b/lib/Sidef/Types/Array/Array.pm
@@ -1099,6 +1099,12 @@ package Sidef::Types::Array::Array {
         $self->_reduce_by('gcd', Sidef::Types::Number::Number::ZERO, sub { $block->run($_[1]) });
     }
 
+    sub gcud_by {
+        my ($self, $block) = @_;
+        $block //= Sidef::Types::Block::Block::IDENTITY;
+        $self->_reduce_by('gcud', Sidef::Types::Number::Number::ZERO, sub { $block->run($_[1]) });
+    }
+
     sub gcd {
         my ($self, $block) = @_;
 
@@ -1109,6 +1115,16 @@ package Sidef::Types::Array::Array {
         Sidef::Types::Number::Number::gcd(@$self);
     }
 
+    sub gcud {
+        my ($self, $block) = @_;
+
+        if (defined($block)) {
+            goto &gcud_by;
+        }
+
+        Sidef::Types::Number::Number::gcud(@$self);
+    }
+
     sub lcm_by {
         my ($self, $block) = @_;
         $block //= Sidef::Types::Block::Block::IDENTITY;

diff --git a/lib/Sidef/Types/Array/Array.pod b/lib/Sidef/Types/Array/Array.pod
@@ -912,6 +912,22 @@ Returns the
 
 =cut
 
+=head2 gcud
+
+    Array.gcud()
+
+Returns the
+
+=cut
+
+=head2 gcud_by
+
+    Array.gcud_by()
+
+Returns the
+
+=cut
+
 =head2 getopt
 
     Array.getopt()

diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -12677,7 +12677,7 @@ package Sidef::Types::Number::Number {
 
             foreach my $z (@terms) {
                 Math::GMPz::Rmpz_gcd($r, $r, $z);
-                Math::GMPz::Rmpz_cmp_ui($r, 1) || last;
+                last if (Math::GMPz::Rmpz_cmp_ui($r, 1) == 0);
             }
 
             return bless \$r;
@@ -12692,6 +12692,39 @@ package Sidef::Types::Number::Number {
         bless \$r;
     }
 
+    sub gcud {    # greatest common unitary divisor (OEIS: A165430)
+        my (@vals) = @_;
+        _valid(\(@vals));
+
+        @vals || return ZERO;    # By convention, gcd of an empty set is 0.
+        @vals == 1 and return $vals[0];
+
+        my @terms = map { _any2mpz($$_) // goto &nan } @vals;
+        my $g     = Math::GMPz::Rmpz_init_set($terms[0]);
+
+        foreach my $i (1 .. $#terms) {
+            Math::GMPz::Rmpz_gcd($g, $g, $terms[$i]);
+            if (Math::GMPz::Rmpz_cmp_ui($g, 1) == 0) {
+                return bless \$g;
+            }
+        }
+
+        state $t = Math::GMPz::Rmpz_init_nobless();
+
+        foreach my $n (@terms) {
+            next if (Math::GMPz::Rmpz_sgn($n) == 0);
+            while (1) {
+                Math::GMPz::Rmpz_divexact($t, $n, $g);
+                Math::GMPz::Rmpz_gcd($t, $t, $g);
+                last if (Math::GMPz::Rmpz_cmp_ui($t, 1) == 0);
+                Math::GMPz::Rmpz_divexact($g, $g, $t);
+            }
+            last if (Math::GMPz::Rmpz_cmp_ui($g, 1) == 0);
+        }
+
+        bless \$g;
+    }
+
     sub gcdext {
         my ($n, $k) = @_;
 
@@ -15616,7 +15649,7 @@ package Sidef::Types::Number::Number {
 
     *exponential_divisors = \&edivisors;
 
-    sub idivisors {
+    sub idivisors {    # OEIS: A077609
         my ($n) = @_;
 
         $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
@@ -2543,7 +2543,7 @@ Returns the n-th term in Stern's diatomic series (or Stern-Brocot sequence).
 
     gcd(...)
 
-The greatest common divisor of a list of integers.
+Returns the greatest common divisor of a list of integers.
 
 =cut
 
@@ -2559,6 +2559,14 @@ The value of C<d> is always non-negative.
 
 =cut
 
+=head2 gcud
+
+    gcud(...)
+
+Returns the greatest common unitary divisor of a list of integers.
+
+=cut
+
 =head2 geometric_sum
 
     geometric_sum(n,r)
@@ -2747,7 +2755,7 @@ When C<a> and C<b> are integers, this is equivalent with:
 
 Returns an array with the infinitary divisors (or i-divisors) of n.
 
-    say 36.idivisors        #=> [1, 4, 9, 36]
+    say 96.idivisors        #=> [1, 2, 3, 6, 16, 32, 48, 96]
 
 Aliases: I<infinitary_divisors>
 

diff --git a/scripts/Tests/divisor_functions.sf b/scripts/Tests/divisor_functions.sf
@@ -106,14 +106,6 @@ do {
 }
 
 do {    # bi-unitary divisors
-    func udivs(n) {
-        n.divisors.grep {|x| gcd(x, n/x) == 1 }
-    }
-
-    func gcud(n, m) {
-        udivs(n) & udivs(m) -> max
-    }
-
     func biudivs(n) {
         n.divisors.grep {|x| gcud(x,n/x) == 1 }
     }
@@ -201,31 +193,6 @@ do {    # non-unitary divisors
 
 do {    # infinitary divisors
 
-    func isidiv(d, n) {
-
-        if (n % d != 0) {
-            return false
-        }
-
-        n.prime_divisors.each {|p|
-
-            var n2 = valuation(n, p).digits(2)
-            var d2 = valuation(d, p).digits(2)
-
-            for j in (0 .. d2.end) {
-                if ((n2[j] == 0) && (d2[j] != 0)) {
-                    return false
-                }
-            }
-        }
-
-        return true;
-    }
-
-    func idivisors(n) {
-        n.divisors.grep{ isidiv(_, n) }
-    }
-
     func a(n, k=1) {
 
         return 0 if (n == 0)