sum_of_remainders.sf 2.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990
  1. #!/usr/bin/ruby
  2. # Daniel "Trizen" Șuteu
  3. # Date: 10 March 2021
  4. # https://github.com/trizen
  5. # Let's consider the following function:
  6. # a(n,v) = Sum_{k=1..n} (v mod k)
  7. # The goal is to compute a(n,v) in sublinear time with respect to v.
  8. # Formula:
  9. # a(n,v) = n*v - A024916(v) + Sum_{k=n+1..v} k*floor(v/k).
  10. # Formula derived from:
  11. # a(n,v) = Sum_{k=1..n} (v - k*floor(v/k))
  12. # = n*v - Sum_{k=1..n} k*floor(v/k)
  13. # = n*v - Sum_{k=1..v} k*floor(v/k) + Sum_{k=n+1..v} k*floor(v/k)
  14. # Related problem:
  15. # Is there a sublinear formula for computing: Sum_{1<=k<=n, gcd(k,n)=1} k*floor(n/k) ?
  16. # See also:
  17. # https://oeis.org/A099726 -- Sum of remainders of the n-th prime mod k, for k = 1,2,3,...,n.
  18. # https://oeis.org/A340976 -- Sum_{1 < k < n} sigma(n) mod k, where sigma = A000203.
  19. # 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.
  20. func T(n) { # Sum_{k=1..n} k = n-th triangular number
  21. n.faulhaber(1)
  22. }
  23. func S(n) { # A024916(n) = Sum_{k=1..n} sigma(k) = Sum_{k=1..n} k*floor(n/k)
  24. #with (n.isqrt) { |s| sum(1..s, {|k| T(idiv(n,k)) + k*idiv(n,k) }) - T(s)*s }
  25. if (n < 0) {
  26. return (T(n.abs) + __FUNC__(n.abs-1))
  27. }
  28. n.dirichlet_sum({1}, {_}, {_}, {.faulhaber(1)})
  29. }
  30. func g(a,b) { # g(a,b) = Sum_{k=a..b} k*floor(b/k)
  31. if (b < 0) {
  32. return (T(b.abs) - T(a.abs-1) + __FUNC__(a, b.abs-1))
  33. }
  34. var total = 0
  35. while (a <= b) {
  36. var t = idiv(b, a)
  37. var u = idiv(b, t)
  38. total += t*(T(u) - T(a-1))
  39. a = u+1
  40. }
  41. return total
  42. }
  43. func sum_remainders(n, v) { # sub-linear formula
  44. sgn(v) * (n*v.abs - S(v) + g(n+1, v))
  45. }
  46. say {|n| sum_remainders(n, n.prime) }.map(1..20) #=> A099726
  47. say {|n| sum_remainders(n-1, n.sigma) }.map(1..20) #=> A340976
  48. for k in (1..8) {
  49. say ("A099726(10^#{k}) = ", sum_remainders(10**k, prime(10**k)))
  50. }
  51. # Positive v
  52. assert_eq(
  53. 20.of {|n| 20.of {|v| sum(1..n, {|k| v % k }) } },
  54. 20.of {|n| 20.of {|v| sum_remainders(n,v) } }
  55. )
  56. # Negative v
  57. assert_eq(
  58. 20.of {|n| 20.of {|v| sum(1..n, {|k| -v % k }) } },
  59. 20.of {|n| 20.of {|v| sum_remainders(n,-v) } }
  60. )
  61. __END__
  62. A099726(10^1) = 30
  63. A099726(10^2) = 2443
  64. A099726(10^3) = 248372
  65. A099726(10^4) = 25372801
  66. A099726(10^5) = 2437160078
  67. A099726(10^6) = 252670261459
  68. A099726(10^7) = 24690625139657
  69. A099726(10^8) = 2516604108737704