- Added the Number k.powerfree_usigma(n, j=1) and `k.powerfree_usig…

…ma0(n)` method. Returns the sum and the count of the k-powerfree divisors of n. When computing the sum, each divisor is raised to the power j. Example: say 20.of { 2.powerfree_usigma(_) } # OEIS: A092261 say 20.of { 2.powerfree_usigma(_, 2) } # OEIS: A091306 say 20.of { 2.powerfree_usigma0(_) } # OEIS: A056671
trizen · Nov 25, 2021 · b13a1b9 · b13a1b9
1 parent ab740cf
commit b13a1b9
Show file tree

Hide file tree

Showing 3 changed files with 126 additions and 4 deletions.
diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -15714,6 +15714,8 @@ package Sidef::Types::Number::Number {
             my ($k, $factor_exp) = @_;
 
             my @d = ($ONE);
+            my $r = Math::GMPz::Rmpz_init();
+
             foreach my $pe (grep { $_->[1] >= $k } @$factor_exp) {
 
                 my ($p, $e) = @$pe;
@@ -15726,10 +15728,12 @@ package Sidef::Types::Number::Number {
 
                 my @t;
                 for (my $i = $k ; $i <= $e ; $i += $k) {
+
+                    Math::GMPz::Rmpz_pow_ui($r, $p, $i);
+
                     foreach my $d (@d) {
                         my $z = Math::GMPz::Rmpz_init();
-                        Math::GMPz::Rmpz_pow_ui($z, $p, $i);
-                        Math::GMPz::Rmpz_mul($z, $z, $d);
+                        Math::GMPz::Rmpz_mul($z, $r, $d);
                         push @t, $z;
                     }
                 }
@@ -17441,6 +17445,71 @@ package Sidef::Types::Number::Number {
         bless \$s;
     }
 
+    sub powerfree_usigma0 {
+        my ($k, $n) = @_;
+
+        # Multiplicative with:
+        #   a(p^e) = 2          # for e < k
+        #   a(p^e) = 1          # for e >= k
+
+        $k = defined($k) ? do { _valid(\$k); _any2ui($$k)   // goto &nan } : 1;
+        $n = defined($n) ? do { _valid(\$n); _big2uistr($n) // goto &nan } : (goto &nan);
+
+        $k > 0 or return ZERO;
+
+        my @factor_exp = _factor_exp($n);
+        @factor_exp and $factor_exp[0][0] eq '0' and return ZERO;
+
+        my $r = Math::GMPz::Rmpz_init();
+        Math::GMPz::Rmpz_setbit($r, scalar grep { $_->[1] < $k } @factor_exp);
+        bless \$r;
+    }
+
+    sub powerfree_usigma {
+        my ($k, $n, $j) = @_;
+
+        # Multiplicative with:
+        #   a(p^e) = p^(e*j) + 1      # for e < k
+        #   a(p^e) = 1                # for e >= k
+
+        $k = defined($k) ? do { _valid(\$k); _any2ui($$k) // goto &nan } : 1;
+        $j = defined($j) ? do { _valid(\$j); _any2ui($$j) // goto &nan } : 1;
+
+        $k > 0 or return ZERO;
+
+        if ($j == 0) {
+            goto &powerfree_usigma0;
+        }
+
+        $n = defined($n) ? do { _valid(\$n); _big2uistr($n) // goto &nan } : (goto &nan);
+
+        my @factor_exp = _factor_exp($n);
+        @factor_exp and $factor_exp[0][0] eq '0' and return ZERO;
+
+        my $t = Math::GMPz::Rmpz_init();
+        my $s = Math::GMPz::Rmpz_init_set_ui(1);
+
+        foreach my $pe (@factor_exp) {
+
+            my ($p, $e) = @$pe;
+
+            $e < $k or next;
+
+            if ($p < ULONG_MAX) {
+                Math::GMPz::Rmpz_ui_pow_ui($t, $p, $e * $j);
+            }
+            else {
+                Math::GMPz::Rmpz_set_str($t, $p, 10);
+                Math::GMPz::Rmpz_pow_ui($t, $t, $e * $j);
+            }
+
+            Math::GMPz::Rmpz_add_ui($t, $t, 1);
+            Math::GMPz::Rmpz_mul($s, $s, $t);
+        }
+
+        bless \$s;
+    }
+
     sub square_usigma0 {
 
         # Multiplicative with:

diff --git a/lib/Sidef/Types/Number/Number.pod b/lib/Sidef/Types/Number/Number.pod
@@ -5486,6 +5486,31 @@ Aliases: I<powerfree_unitary_divisors>, I<unitary_powerfree_divisors>
 
 =cut
 
+=head2 powerfree_usigma
+
+    k.powerfree_usigma(n, j=1)
+
+Returns the sum of the k-powerfree unitary divisors of n, each divisors raised to the power j.
+
+    say 20.of { 2.powerfree_usigma(_) }        # OEIS: A092261
+    say 20.of { 2.powerfree_usigma(_, 2) }     # OEIS: A091306
+
+Equivalent with (but much faster):
+
+    n.udivisors.grep { .is_powerfree(k) }.sum {|d| d**j }
+
+=cut
+
+=head2 powerfree_usigma0
+
+    k.powerfree_usigma0(n)
+
+Returns the number of the k-powerfree unitary divisors of n.
+
+    say 20.of { 2.powerfree_usigma0(_) }       # OEIS: A056671
+
+=cut
+
 =head2 powerful
 
     k.powerful(n)

diff --git a/scripts/Tests/divisor_functions.sf b/scripts/Tests/divisor_functions.sf
@@ -24,7 +24,17 @@ for k in (0..5) {
     }
 }
 
-with (10!) {|n|
+var h2 = Hash(
+    powerfree_udivisors => 'powerfree_usigma',
+)
+
+for k in (0..5), j in (0..5) {
+    h2.each {|a,b|
+        assert_eq(k.(a)(n).sum {|d| d**j }, k.(b)(n,j), "#{a} != #{b} for n=#{n}, k=#{k}, j=#{j}")
+    }
+}
+
+with (irand(5, 10)!) {|n|
     var D = n.divisors
     assert_eq(n.power_divisors,    D.grep{.is_power})
     assert_eq(n.power_divisors(1), D.grep{.is_power(1)})
@@ -35,7 +45,7 @@ with (10!) {|n|
     assert_eq(n.power_divisors(9), D.grep{.is_power(9)})
 }
 
-with (21!) {|n|
+with (irand(10, 21)!) {|n|
     var D = n.udivisors
     assert_eq(n.power_udivisors,    D.grep{.is_power})
     assert_eq(n.power_udivisors(1), D.grep{.is_power(1)})
@@ -46,6 +56,17 @@ with (21!) {|n|
     assert_eq(n.power_udivisors(9), D.grep{.is_power(9)})
 }
 
+with (irand(10, 21)!) {|n|
+    var D = n.udivisors
+    assert_eq(1.powerfree_udivisors(n), D.grep{.is_powerfree(1)})
+    assert_eq(2.powerfree_udivisors(n), D.grep{.is_powerfree(2)})
+    assert_eq(3.powerfree_udivisors(n), D.grep{.is_powerfree(3)})
+    assert_eq(4.powerfree_udivisors(n), D.grep{.is_powerfree(4)})
+    assert_eq(5.powerfree_udivisors(n), D.grep{.is_powerfree(5)})
+    assert_eq(6.powerfree_udivisors(n), D.grep{.is_powerfree(6)})
+    assert_eq(7.powerfree_udivisors(n), D.grep{.is_powerfree(7)})
+}
+
 for n in (0..50) {
     var D1 = n.divisors
     var D2 = n.udivisors
@@ -59,6 +80,13 @@ for n in (0..50) {
     assert_eq(n.power_udivisors(1), D2.grep { .is_power(1) })
     assert_eq(n.power_udivisors(2), D2.grep { .is_power(2) })
     assert_eq(n.power_udivisors(3), D2.grep { .is_power(3) })
+
+    assert_eq(0.powerfree_udivisors(n), D2.grep { .is_powerfree(0) })
+    assert_eq(1.powerfree_udivisors(n), D2.grep { .is_powerfree(1) })
+    assert_eq(2.powerfree_udivisors(n), D2.grep { .is_powerfree(2) })
+    assert_eq(3.powerfree_udivisors(n), D2.grep { .is_powerfree(3) })
+    assert_eq(4.powerfree_udivisors(n), D2.grep { .is_powerfree(4) })
+    assert_eq(5.powerfree_udivisors(n), D2.grep { .is_powerfree(5) })
 }
 
 assert_eq(5040.sigma(-1), 5040.divisors.sum {|d| d**(-1) })