new file: Math/partial_sums_of_gcd-sum_function.sf

Author: trizen

Math/partial_sums_of_gcd-sum_function.sf

@@ -0,0 +1,104 @@
+#!/usr/bin/ruby
+
+# Daniel "Trizen" Șuteu
+# Date: 04 February 2019
+# https://github.com/trizen
+
+# A sublinear algorithm for computing the partial sums of the gcd-sum function, using Dirichlet's hyperbola method.
+
+# 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)
+
+# 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
+
+# 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
+
+func partial_sums_of_gcd_sum_function(n) {
+    var s = n.isqrt
+
+    var mertens_lookup   = [0,1]
+    var euler_sum_lookup = [0,1]
+
+    var lookup_size = n.iroot(3)**2     # O(n^(2/3))
+
+    for i in (1 .. lookup_size) {
+        mertens_lookup[i]   = (mertens_lookup[i-1]   + i.moebius)
+        euler_sum_lookup[i] = (euler_sum_lookup[i-1] + i.euler_phi)
+    }
+
+    func moebius_partial_sum(n) {
+
+        if (n <= lookup_size) {
+            return mertens_lookup[n]
+        }
+
+        n.mertens
+    }
+
+    func euler_phi_partial_sum(n) {
+
+        if (n <= lookup_size) {
+            return euler_sum_lookup[n]
+        }
+
+        var s = n.isqrt
+
+        var A = sum(1..s, {|a|
+            (a * moebius_partial_sum(floor(n/a))) + (moebius(a) * faulhaber(floor(n/a), 1))
+        })
+
+        var C = (moebius_partial_sum(s) * faulhaber(s, 1))
+
+        return (A - C)
+    }
+
+    var A = sum(1..s, {|a|
+        (a * euler_phi_partial_sum(floor(n/a))) + (euler_phi(a) * faulhaber(floor(n/a), 1))
+    })
+
+    var C = (euler_phi_partial_sum(s) * faulhaber(s, 1))
+
+    return (A - C)
+}
+
+func g(n) { n }
+func h(n) { euler_phi(n) }
+
+func test_sum(n, g, h) {
+    sum(1..n, {|k|
+        k.divisors.sum {|d|
+            g(d) * h(k/d)
+        }
+    })
+}
+
+say 20.of { test_sum(_, g, h) }
+say 20.of { partial_sums_of_gcd_sum_function(_) }
+
+__END__
+[0, 1, 4, 9, 17, 26, 41, 54, 74, 95, 122, 143, 183, 208, 247, 292, 340, 373, 436, 473]
+[0, 1, 4, 9, 17, 26, 41, 54, 74, 95, 122, 143, 183, 208, 247, 292, 340, 373, 436, 473]

README.md

@@ -244,6 +244,7 @@ A simple collection of Sidef scripts.
     * [Ore's harmonic numbers](./Math/ore_s_harmonic_numbers.sf)
     * [Partial sums of dedekind psi function](./Math/partial_sums_of_dedekind_psi_function.sf)
     * [Partial sums of euler totient function](./Math/partial_sums_of_euler_totient_function.sf)
+    * [Partial sums of gcd-sum function](./Math/partial_sums_of_gcd-sum_function.sf)
     * [Partial sums of jordan totient function](./Math/partial_sums_of_jordan_totient_function.sf)
     * [Partial sums of prime bigomega function](./Math/partial_sums_of_prime_bigomega_function.sf)
     * [Partial sums of prime omega function](./Math/partial_sums_of_prime_omega_function.sf)