-
Notifications
You must be signed in to change notification settings - Fork 0
/
partial_sums_of_lcm_count_function.sf
72 lines (56 loc) · 1.63 KB
/
partial_sums_of_lcm_count_function.sf
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
#!/usr/bin/ruby
# Let f(n) be the number of couples (x,y) with x and y positive integers, x ≤ y and the least common multiple of x and y equal to n.
# Let a(n) = A007875(n), with a(1) = 1, for n > 1 (due to Vladeta Jovovic, Jan 25 2002):
# a(n) = (1/2)*Sum_{d|n} abs(mu(d))
# = 2^(omega(n)-1)
# = usigma_0(n)/2
# This gives us f(n) as:
# f(n) = Sum_{d|n} a(d)
# This script implements a sub-linear formula for computing partial sums of f(n):
# S(n) = Sum_{k=1..n} f(k)
# = Sum_{k=1..n} Sum_{d|k} a(d)
# = Sum_{k=1..n} a(k) * floor(n/k)
# See also:
# https://oeis.org/A007875
# https://oeis.org/A064608
# https://oeis.org/A182082
# Problem from:
# https://projecteuler.net/problem=379
# Several values for S(10^n):
# S(10^1) = 29
# S(10^2) = 647
# S(10^3) = 11751
# S(10^4) = 186991
# S(10^5) = 2725630
# S(10^6) = 37429395
# S(10^7) = 492143953
# S(10^8) = 6261116500
# S(10^9) = 77619512018
# S(10^10) = 942394656385
# S(10^11) = 11247100884096
func usigma0_partial_sum (n) { # O(sqrt(n)) complexity
var total = 0
var s = n.isqrt
var u = idiv(n, s+1)
var prev = squarefree_count(n)
for k in (1..s) {
var curr = squarefree_count(idiv(n, k+1))
total += (prev - curr)*k
prev = curr
}
u.each_squarefree {|k|
total += idiv(n, k)
}
return total
}
func S(n) {
dirichlet_sum(n,
{|k| (k == 1) ? 1 : usigma0(k)/2 }, # a(n)
{ 1 },
{|k| (usigma0_partial_sum(k)-1)/2 }, # partial sums of a(n)
{ _ }
)
}
for k in (1..5) {
say "S(10^#{k}) = #{S(10**k)}"
}