sidef-scripts/Math/partial_sums_of_gcd-sum_function_fast.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 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]