123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108 |
- #!/usr/bin/ruby
- # Sublinear algorithm for computing the sum of prime powers <= n, based on the sublinear algorithm for computing the sum of primes <= n.
- # See also:
- # https://oeis.org/A074793
- # PARI/GP program:
- #`(
- f(n) = if(n <= 1, return(0)); my(r=sqrtint(n)); my(V=vector(r, k, n\k)); my(L=n\r-1); V=concat(V, vector(L, k, L-k+1)); my(T=vector(#V, k, V[k]*(V[k]+1)\2)); my(S=Map(matrix(#V, 2, x, y, if(y==1, V[x], T[x])))); forprime(p=2, r, my(sp=mapget(S, p-1), p2=p*p); for(k=1, #V, if(V[k] < p2, break); mapput(S, V[k], mapget(S, V[k]) - p*(mapget(S, V[k]\p) - sp)))); mapget(S, n)-1;
- g(n, j=2) = if(n <= 1, return(0)); my(r=sqrtint(n)); my(V=vector(r, k, n\k)); my(F(n, j)=(subst(bernpol(j+1), x, n+1) - subst(bernpol(j+1), x, 1)) / (j+1)); my(L=n\r-1); V=concat(V, vector(L, k, L-k+1)); my(T=vector(#V, k, F(V[k], j))); my(S=Map(matrix(#V, 2, x, y, if(y==1, V[x], T[x])))); forprime(p=2, r, my(sp=mapget(S, p-1), p2=p*p); for(k=1, #V, if(V[k] < p2, break); mapput(S, V[k], mapget(S, V[k]) - p^j*(mapget(S, V[k]\p) - sp)))); mapget(S, n)-1;
- a(n) = my(s=0); forprime(p=2, sqrtnint(n,3), s += (p^(logint(n, p)+1) - 1)/(p-1) - p*p - p - 1); f(n) + g(sqrtint(n)) + s;
- b(n) = my(s=0); forprime(p=2, sqrtint(n), s += (p^(logint(n, p)+1) - 1)/(p-1) - p - 1); f(n) + s;
- c(n) = sum(k=1, logint(n,2), g(sqrtnint(n,k), k));
- )
- func sum_of_primes(n, j=1) { # Sum_{p prime <= n} p^k
- return n.sum_primes if (j == 1) # optimization
- return 0 if (n <= 1)
- var r = n.isqrt
- var V = (1..r -> map {|k| idiv(n,k) })
- V << range(V.last-1, 1, -1).to_a...
- var S = Hash(V.map{|k| (k, faulhaber_sum(k,j)) }...)
- for p in (2..r) {
- S{p} > S{p-1} || next
- var sp = S{p-1}
- var p2 = p*p
- V.each {|v|
- break if (v < p2)
- S{v} -= ipow(p,j)*(S{idiv(v,p)} - sp)
- }
- }
- return S{n}-1
- }
- func sum_of_prime_powers(n) {
- # a(n) = Sum_{p prime <= n} p
- # b(n) = Sum_{p prime <= n^(1/2)} p^2
- # c(n) = Sum_{p prime <= n^(1/3)} f(p)
- # sum_of_prime_powers(n) = a(n) + b(n) + c(n)
- var ps1 = sum_of_primes(n)
- var ps2 = sum_of_primes(n.isqrt, 2)
- # f(p) = (Sum_{k=1..floor(log_p(n))} p^k) - p^2 - p
- # = (p^(1+floor(log_p(n))) - 1)/(p-1) - p^2 - p - 1
- var ps3 = n.icbrt.primes.sum {|p|
- (ipow(p, n.ilog(p)+1) - 1)/(p-1) - p*p - p - 1
- }
- ps1 + ps2 + ps3
- }
- func sum_of_prime_powers_2(n) {
- # a(n) = Sum_{p prime <= n} p
- # b(n) = Sum_{p prime <= n^(1/2)} f(p)
- # sum_of_prime_powers(n) = a(n) + b(n)
- var ps1 = sum_of_primes(n)
- # f(p) = (Sum_{k=1..floor(log_p(n))} p^k) - p
- # = (p^(1+floor(log_p(n))) - 1)/(p-1) - p - 1
- var ps2 = 0
- n.isqrt.each_prime {|p|
- ps2 += ((ipow(p, n.ilog(p)+1) - 1)/(p-1) - p - 1)
- }
- ps1 + ps2
- }
- func sum_of_prime_powers_3(n) {
- # a(n) = Sum_{k=1..floor(log_2(n))} Sum_{p prime <= n^(1/k)} p^k.
- sum(1..n.ilog2, {|k|
- sum_of_primes(n.iroot(k), k)
- })
- }
- say sum_of_prime_powers(1e5) #=> 457028152
- say sum_of_prime_powers(43**3) #=> 294752679
- for k in (1..100) {
- var n = 1e3.irand
- var x = sum_of_prime_powers(n)
- var y = sum_of_prime_powers_2(n)
- var z = sum_of_prime_powers_3(n)
- var w = n.prime_powers.sum
- assert_eq(x,y)
- assert_eq(x,z)
- assert_eq(x,w)
- }
|