sidef-scripts/Math/sum_of_prime_powers.sf

#!/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)
}