-
Notifications
You must be signed in to change notification settings - Fork 0
/
sum_of_remainders.sf
89 lines (69 loc) · 2.39 KB
/
sum_of_remainders.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
#!/usr/bin/ruby
# Daniel "Trizen" Șuteu
# Date: 10 March 2021
# https://github.com/trizen
# Let's consider the following function:
# a(n,v) = Sum_{k=1..n} (v mod k)
# The goal is to compute a(n,v) in sublinear time with respect to v.
# Formula:
# a(n,v) = n*v - A024916(v) + Sum_{k=n+1..v} k*floor(v/k).
# Formula derived from:
# a(n,v) = Sum_{k=1..n} (v - k*floor(v/k))
# = n*v - Sum_{k=1..n} k*floor(v/k)
# = n*v - Sum_{k=1..v} k*floor(v/k) + Sum_{k=n+1..v} k*floor(v/k)
# Related problem:
# Is there a sublinear formula for computing: Sum_{1<=k<=n, gcd(k,n)=1} k*floor(n/k) ?
# See also:
# https://oeis.org/A099726 -- Sum of remainders of the n-th prime mod k, for k = 1,2,3,...,n.
# https://oeis.org/A340976 -- Sum_{1 < k < n} sigma(n) mod k, where sigma = A000203.
# https://oeis.org/A340180 -- a(n) = Sum_{x in C(n)} (sigma(n) mod x), where C(n) is the set of numbers < n coprime to n, and sigma = A000203.
func T(n) { # Sum_{k=1..n} k = n-th triangular number
n.faulhaber(1)
}
func S(n) { # A024916(n) = Sum_{k=1..n} sigma(k) = Sum_{k=1..n} k*floor(n/k)
#with (n.isqrt) { |s| sum(1..s, {|k| T(idiv(n,k)) + k*idiv(n,k) }) - T(s)*s }
if (n < 0) {
return (T(n.abs) + __FUNC__(n.abs-1))
}
n.dirichlet_sum({1}, {_}, {_}, {.faulhaber(1)})
}
func g(a,b) { # g(a,b) = Sum_{k=a..b} k*floor(b/k)
if (b < 0) {
return (T(b.abs) - T(a.abs-1) + __FUNC__(a, b.abs-1))
}
var total = 0
while (a <= b) {
var t = idiv(b, a)
var u = idiv(b, t)
total += t*(T(u) - T(a-1))
a = u+1
}
return total
}
func sum_remainders(n, v) { # sub-linear formula
sgn(v) * (n*v.abs - S(v) + g(n+1, v))
}
say {|n| sum_remainders(n, n.prime) }.map(1..20) #=> A099726
say {|n| sum_remainders(n-1, n.sigma) }.map(1..20) #=> A340976
for k in (1..8) {
say ("A099726(10^#{k}) = ", sum_remainders(10**k, prime(10**k)))
}
# Positive v
assert_eq(
20.of {|n| 20.of {|v| sum(1..n, {|k| v % k }) } },
20.of {|n| 20.of {|v| sum_remainders(n,v) } }
)
# Negative v
assert_eq(
20.of {|n| 20.of {|v| sum(1..n, {|k| -v % k }) } },
20.of {|n| 20.of {|v| sum_remainders(n,-v) } }
)
__END__
A099726(10^1) = 30
A099726(10^2) = 2443
A099726(10^3) = 248372
A099726(10^4) = 25372801
A099726(10^5) = 2437160078
A099726(10^6) = 252670261459
A099726(10^7) = 24690625139657
A099726(10^8) = 2516604108737704