new file: Math/partial_sums_of_generalized_gcd-sum_function.pl

Author: trizen

Math/partial_sums_of_generalized_gcd-sum_function.pl

@@ -0,0 +1,131 @@
+#!/usr/bin/perl
+
+# Daniel "Trizen" Șuteu
+# Date: 25 May 2025
+# https://github.com/trizen
+
+# A sublinear algorithm for computing the partial sums of the generalized gcd-sum function, using Dirichlet's hyperbola method.
+
+# Generalized Pillai's function:
+#   pillai(n,k) = Sum_{d|n} mu(n/d) * d^k * tau(d)
+
+# Multiplicative formula for Sum_{1 <= x_1, x_2, ..., x_k <= n} gcd(x_1, x_2, ..., x_k, n)^k:
+#   a(p^e) = (e - e/p^k + 1) * p^(k*e) = p^((e - 1) * k) * (p^k + e*(p^k - 1))
+
+# The partial sums of the gcd-sum function is defined as:
+#
+#   a(n) = Sum_{k=1..n} Sum_{d|k} d*phi(k/d)
+#
+# where phi(k) is the Euler totient function.
+
+# Also equivalent with:
+#   a(n) = Sum_{j=1..n} Sum_{i=1..j} gcd(i, j)
+
+# Based on the formula:
+#   a(n) = (1/2)*Sum_{k=1..n} phi(k) * floor(n/k) * floor(1+n/k)
+
+# Generalized formula:
+#   a(n,k) = Sum_{x=1..n} J_k(x) * F_k(floor(n/x))
+# where F_k(n) are the Faulhaber polynomials: F_k(n) = Sum_{x=1..n} x^k.
+
+# Example:
+#   a(10^1) = 122
+#   a(10^2) = 18065
+#   a(10^3) = 2475190
+#   a(10^4) = 317257140
+#   a(10^5) = 38717197452
+#   a(10^6) = 4571629173912
+#   a(10^7) = 527148712519016
+#   a(10^8) = 59713873168012716
+#   a(10^9) = 6671288261316915052
+
+#   a(10^1, 2) = 1106
+#   a(10^2, 2) = 1598361
+#   a(10^3, 2) = 2193987154
+#   a(10^4, 2) = 2828894776292
+#   a(10^5, 2) = 3466053625977000
+#   a(10^6, 2) = 4104546122851466704
+#   a(10^7, 2) = 4742992578252739471520
+#   a(10^8, 2) = 5381500783126483704718848
+#   a(10^9, 2) = 6020011093886996189443484608
+
+# OEIS sequences:
+#   https://oeis.org/A272718 -- Partial sums of gcd-sum sequence A018804.
+#   https://oeis.org/A018804 -- Pillai's arithmetical function: Sum_{k=1..n} gcd(k, n).
+
+# See also:
+#   https://en.wikipedia.org/wiki/Dirichlet_hyperbola_method
+#   https://trizenx.blogspot.com/2018/11/partial-sums-of-arithmetical-functions.html
+
+use 5.020;
+use strict;
+use warnings;
+
+use experimental qw(signatures);
+use Math::AnyNum qw(faulhaber_sum ipow);
+use ntheory      qw(jordan_totient sqrtint rootint);
+
+sub partial_sums_of_gcd_sum_function($n, $m) {
+
+    my $s                  = sqrtint($n);
+    my @totient_sum_lookup = (0);
+
+    my $lookup_size    = 2 + 2 * rootint($n, 3)**2;
+    my @jordan_totient = (0);
+
+    foreach my $x (1 .. $lookup_size) {
+        push @jordan_totient, jordan_totient($m, $x);
+    }
+
+    foreach my $i (1 .. $lookup_size) {
+        $totient_sum_lookup[$i] = $totient_sum_lookup[$i - 1] + $jordan_totient[$i];
+    }
+
+    my %seen;
+
+    my sub totient_partial_sum($n) {
+
+        if ($n <= $lookup_size) {
+            return $totient_sum_lookup[$n];
+        }
+
+        if (exists $seen{$n}) {
+            return $seen{$n};
+        }
+
+        my $s = sqrtint($n);
+        my $T = ${faulhaber_sum($n, $m)};
+
+        foreach my $k (2 .. int($n / ($s + 1))) {
+            $T -= __SUB__->(int($n / $k));
+        }
+
+        foreach my $k (1 .. $s) {
+            $T -= (int($n / $k) - int($n / ($k + 1))) * $totient_sum_lookup[$k];
+        }
+
+        $seen{$n} = $T;
+    }
+
+    my $A = 0;
+
+    foreach my $k (1 .. $s) {
+        my $t = int($n / $k);
+        $A += ${ipow($k, $m)} * totient_partial_sum($t) + $jordan_totient[$k] * ${faulhaber_sum($t, $m)};
+    }
+
+    my $T = ${faulhaber_sum($s, $m)};
+    my $C = totient_partial_sum($s);
+
+    return ($A - $T * $C);
+}
+
+foreach my $n (1 .. 8) {    # takes less than 1 second
+    say "a(10^$n, 1) = ", partial_sums_of_gcd_sum_function(10**$n, 1);
+}
+
+say '';
+
+foreach my $n (1 .. 8) {    # takes less than 1 second
+    say "a(10^$n, 2) = ", partial_sums_of_gcd_sum_function(10**$n, 2);
+}

README.md

@@ -1013,6 +1013,7 @@ A nice collection of day-to-day Perl scripts.
     * [Partial sums of gcd-sum function](./Math/partial_sums_of_gcd-sum_function.pl)
     * [Partial sums of gcd-sum function fast](./Math/partial_sums_of_gcd-sum_function_fast.pl)
     * [Partial sums of gcd-sum function faster](./Math/partial_sums_of_gcd-sum_function_faster.pl)
+    * [Partial sums of generalized gcd-sum function](./Math/partial_sums_of_generalized_gcd-sum_function.pl)
     * [Partial sums of gpf](./Math/partial_sums_of_gpf.pl)
     * [Partial sums of inverse moebius transform of dedekind function](./Math/partial_sums_of_inverse_moebius_transform_of_dedekind_function.pl)
     * [Partial sums of jordan totient function](./Math/partial_sums_of_jordan_totient_function.pl)