sidef-scripts/Math/partial_sums_of_generalized_gcd-sum_function.sf at master · trizen/sidef-scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
#!/usr/bin/ruby

# 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

func partial_sums_of_gcd_sum_function(n, m=1) {

    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] = jordan_totient(k, m)))
    }

    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(m)

        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.ipow(m) * euler_phi_partial_sum(t)) + (euler_phi_lookup[k] * t.faulhaber(m))
    })

    var T = s.faulhaber(m)
    var C = euler_phi_partial_sum(s)

    return (A - T*C)
}

func partial_sums_of_gcd_sum_function_dirichlet(n, m=1) {
    dirichlet_hyperbola(n, { .totient(m) }, { .ipow(m) }, { .totient_sum(m) }, { .faulhaber(m) })
}

func sigma_ktm_partial_sum(n, m) {      # O(sqrt(n)) complexity

    var total = 0

    var s = n.isqrt

    for k in (1 .. s) {
        total += 2*(k**m * faulhaber_sum(idiv(n,k), m))
    }

    total -= faulhaber_sum(s, m)**2

    return total
}

func partial_sums_of_gcd_sum_function_dirichlet_2(n,m=1) {

    # See also:
    #   https://oeis.org/A372931

    dirichlet_hyperbola(n, { .mu }, { .ipow(m) * .tau }, { .mertens }, { sigma_ktm_partial_sum(_, m) })
}

say 20.of { partial_sums_of_gcd_sum_function(_) }
say 20.of { partial_sums_of_gcd_sum_function_dirichlet(_) }
say 20.of { partial_sums_of_gcd_sum_function_dirichlet_2(_) }

say ''

say 20.of { partial_sums_of_gcd_sum_function(_,2) }
say 20.of { partial_sums_of_gcd_sum_function_dirichlet(_,2) }
say 20.of { partial_sums_of_gcd_sum_function_dirichlet_2(_,2) }

__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]
[0, 1, 4, 9, 17, 26, 41, 54, 74, 95, 122, 143, 183, 208, 247, 292, 340, 373, 436, 473]

[0, 1, 8, 25, 65, 114, 233, 330, 538, 763, 1106, 1347, 2027, 2364, 3043, 3876, 4900, 5477, 7052, 7773]
[0, 1, 8, 25, 65, 114, 233, 330, 538, 763, 1106, 1347, 2027, 2364, 3043, 3876, 4900, 5477, 7052, 7773]
[0, 1, 8, 25, 65, 114, 233, 330, 538, 763, 1106, 1347, 2027, 2364, 3043, 3876, 4900, 5477, 7052, 7773]