partial_sums_of_gcd-sum_function.sf

#!/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 lookup_size = (2 * n.iroot(3)**2)     # O(n^(2/3))

    var mertens_lookup   = [0,1]
    var euler_sum_lookup = [0,1]

    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]   + moebius_lookup[i])
        euler_sum_lookup[i] = (euler_sum_lookup[i-1] + (euler_phi_lookup[i] = i.euler_phi))
    }

    var mertens_cache = Hash()

    func moebius_partial_sum (n) {

        if (n <= lookup_size) {
            return mertens_lookup[n]
        }

        mertens_cache{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, {|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, {|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))

    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]