new file: Math/partial_sums_of_gcd-sum_function_fast.sf

Author: trizen

Math/partial_sums_of_gcd-sum_function.sf

@@ -40,23 +40,28 @@
 func partial_sums_of_gcd_sum_function(n) {
     var s = n.isqrt
 
+    var lookup_size = (2 * n.iroot(3)**2)     # O(n^(2/3))
+
     var mertens_lookup   = [0,1]
     var euler_sum_lookup = [0,1]
 
-    var lookup_size = n.iroot(3)**2     # O(n^(2/3))
+    var euler_phi_lookup = [0]
+    var moebius_lookup   = ::moebius(0, lookup_size)
 
     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)
+        mertens_lookup[i]   = (mertens_lookup[i-1]   + moebius_lookup[i])
+        euler_sum_lookup[i] = (euler_sum_lookup[i-1] + (euler_phi_lookup[i] = i.euler_phi))
     }
 
-    func moebius_partial_sum(n) {
+    var mertens_cache = Hash()
+
+    func moebius_partial_sum (n) {
 
         if (n <= lookup_size) {
             return mertens_lookup[n]
         }
 
-        n.mertens
+        mertens_cache{n} := n.mertens
     }
 
     func euler_phi_partial_sum(n) {
@@ -67,17 +72,19 @@ func partial_sums_of_gcd_sum_function(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 A = sum(1..s, {|k|
+            var t = floor(n/k)
+            (k * moebius_partial_sum(t)) + (moebius_lookup[k] * t.faulhaber(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 A = sum(1..s, {|k|
+        var t = floor(n/k)
+        (k * euler_phi_partial_sum(t)) + (euler_phi_lookup[k] * t.faulhaber(1))
     })
 
     var C = (euler_phi_partial_sum(s) * faulhaber(s, 1))

Math/partial_sums_of_gcd-sum_function_fast.sf

@@ -0,0 +1,92 @@
+#!/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 euler_sum_lookup = [0]
+
+    var lookup_size      = (2 + 2*n.iroot(3)**2)
+    var euler_phi_lookup = [0]
+
+    for k in (1 .. lookup_size) {
+        euler_sum_lookup[k] = (euler_sum_lookup[k-1] + (euler_phi_lookup[k] = k.euler_phi))
+    }
+
+    var seen = Hash()
+
+    func euler_phi_partial_sum(n) {
+
+        if (n <= lookup_size) {
+            return euler_sum_lookup[n]
+        }
+
+        if (seen.has(n)) {
+            return seen{n}
+        }
+
+        var s = n.isqrt
+        var T = n.faulhaber(1)
+
+        var A = sum(2..s, {|k|
+            __FUNC__(floor(n/k))
+        })
+
+        var B = sum(1 .. floor(n/s)-1, {|k|
+            (floor(n/k) - floor(n/(k+1))) * __FUNC__(k)
+        })
+
+        seen{n} = (T - A - B)
+    }
+
+    var A = sum(1..s, {|k|
+        var t = floor(n/k)
+        (k * euler_phi_partial_sum(t)) + (euler_phi_lookup[k] * t.faulhaber(1))
+    })
+
+    var T = s.faulhaber(1)
+    var C = euler_phi_partial_sum(s)
+
+    return (A - T*C)
+}
+
+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]

README.md

@@ -245,6 +245,7 @@ A simple collection of Sidef scripts.
     * [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 gcd-sum function fast](./Math/partial_sums_of_gcd-sum_function_fast.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)