new file: Math/partial_sums_of_sigma0_function.pl

modified:   Math/partial_sums_of_sigma_function.pl

Author: trizen

Math/partial_sums_of_sigma0_function.pl

@@ -0,0 +1,63 @@
+#!/usr/bin/perl
+
+# Author: Daniel "Trizen" Șuteu
+# Date: 09 November 2018
+# Edit: 30 March 2025
+# https://github.com/trizen
+
+# Algorithm with O(sqrt(n)) complexity for computing the partial-sums of the `sigma_0(k)` function:
+#   Sum_{k=1..n} sigma_0(k)
+
+# See also:
+#   https://oeis.org/A006218
+#   https://en.wikipedia.org/wiki/Divisor_function
+#   https://en.wikipedia.org/wiki/Faulhaber%27s_formula
+#   https://en.wikipedia.org/wiki/Bernoulli_polynomials
+#   https://trizenx.blogspot.com/2018/11/partial-sums-of-arithmetical-functions.html
+
+use 5.036;
+
+sub sigma0_partial_sum_faulhaber ($n) {
+
+    my $s   = int(sqrt($n));
+    my $sum = 0;
+
+    foreach my $k (1 .. $s) {
+        $sum += 2 * int($n / $k);
+    }
+
+    return ($sum - $s * $s);
+}
+
+sub sigma0_partial_sum_test ($n) {    # just for testing
+    my $sum = 0;
+    foreach my $k (1 .. $n) {
+        $sum += int($n / $k);
+    }
+    return $sum;
+}
+
+foreach my $m (0 .. 10) {
+
+    my $n = int(rand(1 << (2 * $m)));
+
+    my $t1 = sigma0_partial_sum_test($n);
+    my $t2 = sigma0_partial_sum_faulhaber($n);
+
+    say "Sum_{k=1..$n} sigma_0(k) = $t2";
+
+    die "error: $t1 != $t2" if ($t1 != $t2);
+}
+
+__END__
+Sum_{k=1..0} sigma_0(k) = 0
+Sum_{k=1..3} sigma_0(k) = 5
+Sum_{k=1..13} sigma_0(k) = 37
+Sum_{k=1..30} sigma_0(k) = 111
+Sum_{k=1..193} sigma_0(k) = 1049
+Sum_{k=1..51} sigma_0(k) = 211
+Sum_{k=1..2288} sigma_0(k) = 18059
+Sum_{k=1..15985} sigma_0(k) = 157208
+Sum_{k=1..10112} sigma_0(k) = 94818
+Sum_{k=1..152099} sigma_0(k) = 1838389
+Sum_{k=1..446108} sigma_0(k) = 5872025

Math/partial_sums_of_sigma_function.pl

@@ -2,6 +2,7 @@
 
 # Author: Daniel "Trizen" Șuteu
 # Date: 09 November 2018
+# Edit: 30 March 2025
 # https://github.com/trizen
 
 # A new generalized algorithm with O(sqrt(n)) complexity for computing the partial-sums of the `sigma_j(k)` function:
@@ -16,23 +17,19 @@
 #   https://en.wikipedia.org/wiki/Bernoulli_polynomials
 #   https://trizenx.blogspot.com/2018/11/partial-sums-of-arithmetical-functions.html
 
-use 5.020;
-use strict;
-use warnings;
-
-use ntheory qw(divisors);
-use experimental qw(signatures);
+use 5.036;
+use ntheory      qw(divisors);
 use Math::AnyNum qw(faulhaber_sum bernoulli sum isqrt ipow);
 
-sub sigma_partial_sum_faulhaber ($n, $m = 1) {       # using Faulhaber's formula
+sub sigma_partial_sum_faulhaber ($n, $m = 1) {    # using Faulhaber's formula
 
     my $s = isqrt($n);
     my $u = int($n / ($s + 1));
 
     my $sum = 0;
 
     foreach my $k (1 .. $s) {
-        $sum += $k * (faulhaber_sum(int($n/$k), $m) - faulhaber_sum(int($n/($k+1)), $m));
+        $sum += $k * (faulhaber_sum(int($n / $k), $m) - faulhaber_sum(int($n / ($k + 1)), $m));
     }
 
     foreach my $k (1 .. $u) {
@@ -42,15 +39,30 @@ ($n, $m = 1)
     return $sum;
 }
 
-sub sigma_partial_sum_bernoulli ($n, $m = 1) {       # using Bernoulli polynomials
+sub sigma_partial_sum_dirichlet ($n, $m = 1) {    # using the Dirichlet hyperbola method
+
+    my $total = 0;
+    my $s     = isqrt($n);
+
+    for my $k (1 .. $s) {
+        $total += faulhaber_sum(int($n / $k), $m);
+        $total += ipow($k, $m) * int($n / $k);
+    }
+
+    $total -= $s * faulhaber_sum($s, $m);
+
+    return $total;
+}
+
+sub sigma_partial_sum_bernoulli ($n, $m = 1) {    # using Bernoulli polynomials
 
     my $s = isqrt($n);
     my $u = int($n / ($s + 1));
 
     my $sum = 0;
 
     foreach my $k (1 .. $s) {
-        $sum += $k * (bernoulli($m+1, 1+int($n/$k)) - bernoulli($m+1, 1+int($n/($k+1)))) / ($m+1);
+        $sum += $k * (bernoulli($m + 1, 1 + int($n / $k)) - bernoulli($m + 1, 1 + int($n / ($k + 1)))) / ($m + 1);
     }
 
     foreach my $k (1 .. $u) {
@@ -61,7 +73,11 @@ ($n, $m = 1)
 }
 
 sub sigma_partial_sum_test ($n, $m = 1) {    # just for testing
-    sum(map { sum(map { ipow($_, $m) } divisors($_)) } 1..$n);
+    sum(
+        map {
+            sum(map { ipow($_, $m) } divisors($_))
+          } 1 .. $n
+       );
 }
 
 foreach my $m (0 .. 10) {
@@ -71,11 +87,13 @@ ($n, $m = 1)
     my $t1 = sigma_partial_sum_test($n, $m);
     my $t2 = sigma_partial_sum_faulhaber($n, $m);
     my $t3 = sigma_partial_sum_bernoulli($n, $m);
+    my $t4 = sigma_partial_sum_dirichlet($n, $m);
 
     say "Sum_{k=1..$n} sigma_$m(k) = $t2";
 
     die "error: $t1 != $t2" if ($t1 != $t2);
     die "error: $t1 != $t3" if ($t1 != $t3);
+    die "error: $t1 != $t4" if ($t1 != $t4);
 }
 
 __END__