-
Notifications
You must be signed in to change notification settings - Fork 0
/
sum_of_number_of_unitary_divisors.sf
105 lines (76 loc) · 3.53 KB
/
sum_of_number_of_unitary_divisors.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
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
#!/usr/bin/ruby
# Daniel "Trizen" Șuteu
# Date: 10 January 2019
# https://github.com/trizen
# Two fast algorithms for computing the sum of number of unitary divisors from 1 to n.
# a(n) = Sum_{k=1..n} usigma_0(k)
# Based on the formula:
# a(n) = Sum_{k=1..n} moebius(k)^2 * floor(n/k)
# See also:
# https://oeis.org/A034444 -- Partial sums of A034444: sum of number of unitary divisors from 1 to n.
# https://oeis.org/A180361 -- Sum of number of unitary divisors (A034444) from 1 to 10^n
# https://oeis.org/A268732 -- Sum of the numbers of divisors of gcd(x,y) with x*y <= n.
# Asymptotic formula:
# a(n) ~ n*log(n)/zeta(2) + O(n)
# Better asymptotic formula:
# a(n) ~ (n/zeta(2))*(log(n) + 2*γ - 1 - c) + O(sqrt(n) * log(n))
#
# where γ is the Euler-Mascheroni constant and c = 2*zeta'(2)/zeta(2) = -1.1399219861890656127997287200...
func asymptotic_formula(n) {
# c = 2*Zeta'(2)/Zeta(2) = (12 * Zeta'(2))/π^2 = 2 (-12 log(A) + γ + log(2) + log(π))
const c = -1.13992198618906561279972872003946000480696456161386195911639472087583455473348121357
# Asymptotic formula based on Merten's theorem (1874) (see: https://oeis.org/A064608)
(n/zeta(2)) * (log(n) + 2*Num.EulerGamma - 1 - c)
}
func usigma0_partial_sum_1 (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 usigma0_partial_sum_2 (n) { # based on formula by Jerome Raulin (https://oeis.org/A064608)
var total = 0
n.isqrt.each_squarefree {|k|
var v = moebius(k)
var t = 2*sum(1..isqrt(idiv(n, k*k)), {|j|
idiv(n, j*k*k)
})
total += v*(t - isqrt(idiv(n, k*k))**2)
}
return total
}
func usigma0_partial_sum_3(n) { # O(sqrt(n)) complexity, using Dirichlet's hyperbola method
n.dirichlet_sum({1}, {.mu.abs}, {_}, {.squarefree_count})
}
say 20.of(usigma0_partial_sum_1)
say 20.of(usigma0_partial_sum_2)
say 20.of(usigma0_partial_sum_3)
for k in (0..7) {
var n = 10**k
var t = usigma0_partial_sum_1(n)
var u = asymptotic_formula(n)
printf("a(10^%s) = %10s ~ %-15s -> %s\n", k, t, u.round(-2), t/u)
}
__END__
[0, 1, 3, 5, 7, 9, 13, 15, 17, 19, 23, 25, 29, 31, 35, 39, 41, 43, 47, 49]
[0, 1, 3, 5, 7, 9, 13, 15, 17, 19, 23, 25, 29, 31, 35, 39, 41, 43, 47, 49]
a(10^0) = 1 ~ 0.79 -> 1.27085398285349342897812915198984638968899591751
a(10^1) = 23 ~ 21.87 -> 1.05182461403816051734935994402113331145060974294
a(10^2) = 359 ~ 358.65 -> 1.00098140095602073835866744824992972185806123685
a(10^3) = 4987 ~ 4986.28 -> 1.00014357239778054254970740667091143421188177813
a(10^4) = 63869 ~ 63860.88 -> 1.00012715302552355451250212258735392366329621935
a(10^5) = 778581 ~ 778589.19 -> 0.999989484576929013867264739526374966823956960403
a(10^6) = 9185685 ~ 9185695.75 -> 0.99999882923368455522780513812504287278271814501
a(10^7) = 105854997 ~ 105854996.37 -> 1.00000000598372061072117962943109677794267023891
a(10^8) = 1198530315 ~ 1198530351.90 -> 0.999999969211002320383540850995519903094748492418
a(10^9) = 13385107495 ~ 13385107401.37 -> 1.00000000699496540213133746406895764726726792391
a(10^10) = 147849112851 ~ 147849112837.28 -> 1.00000000009281141854332921757852421030396550125