12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879 |
- #!/usr/bin/ruby
- # Sum of remainders of the n-th prime mod k, for k = 1,2,3,...,n.
- # https://oeis.org/A099726
- # Formula:
- # a(n) = n*p - A024916(p) + Sum_{k=n+1..p} k*floor(p/k), where p = prime(n).
- # Some large terms:
- # a(10^1) = 30
- # a(10^2) = 2443
- # a(10^3) = 248372
- # a(10^4) = 25372801
- # a(10^5) = 2437160078
- # a(10^6) = 252670261459
- # a(10^7) = 24690625139657
- # a(10^8) = 2516604108737704
- # a(10^9) = 249876964098609078
- # a(10^10) = 24994462548503343285
- # a(10^11) = 2513170619960257982344
- func T(n) { # Sum_{k=1..n} k = n-th triangular number
- n.faulhaber(1)
- }
- func S(n) { # Sum_{k=1..n} sigma(k) = Sum_{k=1..n} k*floor(n/k)
- #var s = n.isqrt
- #sum(1..s, {|k| T(idiv(n,k)) + k*idiv(n,k) }) - T(s)*s
- n.dirichlet_sum({1}, {_}, {_}, {.faulhaber(1)})
- }
- func g(a,b) {
- 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 a(n) { # sub-linear formula
- var p = n.prime
- n*p - S(p) + g(n+1, p)
- }
- func a_linear(n) { # just for testing
- var p = n.prime
- sum(1..n, {|k| p % k })
- }
- assert_eq(a.map(1..100), a_linear.map(1..100))
- say a.map(1..58) #=> [0, 1, 3, 5, 7, 7, 14, 18, 28, 30, 31, 26, 38, 45, 63, 71, 93, 75, 96, 115, 101, 142, 161, 167, 152, 159, 203, 224, 219, 222, 216, 250, 263, 296, 341, 320, 319, 349, 433, 427, 496, 419, 487, 481, 538, 537, 495, 631, 635, 676, 697, 777, 665, 820, 784, 874, 929, 856]
- for k in (1..9) {
- say "a(10^#{k}) = #{a(10**k)}"
- }
- __END__
- # PARI/GP program (linear time)
- a(n) = my(p=prime(n)); sum(k=1, n, p%k);
- # PARI/GP program (sublinear time)
- T(n) = n*(n+1)/2;
- S(n) = my(s=sqrtint(n)); sum(k=1, s, T(n\k) + k*(n\k)) - s*T(s); \\ A024916
- g(a,b) = my(s=0); while(a <= b, my(t=b\a); my(u=b\t); s += t*(T(u) - T(a-1)); a = u+1); s;
- a(n) = my(p=prime(n)); n*p - S(p) + g(n+1, p);
- # PARI/GP program (sublinear time) -- shorter, but slower
- T(n) = n*(n+1)/2;
- g(a,b) = my(s=0); while(a <= b, my(t=b\a); my(u=b\t); s += t*(T(u) - T(a-1)); a = u+1); s;
- a(n) = my(p=prime(n)); n*p - g(1, p) + g(n+1, p);
|