sidef-scripts/Math/sum_of_primes.sf

#!/usr/bin/ruby

# Sublinear algorithm for computing the sum of primes <= n, each prime raised to a fixed power j >= 0.

# Algorithm from:
#   https://math.stackexchange.com/questions/1378286/find-the-sum-of-all-primes-smaller-than-a-big-number

# PARI/GP program (for j = 1):
#   a(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;

# Generalized PARI/GP program (for j >= 0):
#   a(n, j=2) = if(n <= 1, return(0)); my(r=sqrtint(n)); my(V=vector(r, k, n\k)); my(B=bernpol(j+1)); my(F(n)=(subst(B, x, n+1) - subst(B, 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]))); 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;

func sum_of_primes(n, j=1) {

    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(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
}

say sum_of_primes(1e6)         #=> 37550402023
say sum_of_primes(1e6, 2)      #=> 24693298341834533

assert_eq(
    30.of { sum_of_primes(_) }
    30.of { .sum_primes }
)

assert_eq(
    30.of { sum_of_primes(_, 2) }
    30.of { .primes.sum{.sqr} }
)

assert_eq(
    30.of { sum_of_primes(_, 0) }
    30.of { .prime_count }
)

__END__

# A failed attempt at creating a sublinear method for computing the sum of primes <= n.

# Based on the formulas:
#   a(n) = Sum_{k=1..n} Sum_{d|k} A008683(d) * A008472(k/d)
#   a(n) = Sum_{k=1..n} k*Sum_{d|k} mu(d) * omega(k/d)
#   a(n) = Sum_{k=1..n} floor(n/k) * Sum_{p prime | k} p*mu(k/p)

# Which can be computed in sublinear time as:
#   a(n) = Sum_{k=1..floor(sqrt(n))} (A008472(k)*A002321(floor(n/k)) + A008683(k)*A024924(floor(n/k))) - A002321(floor(sqrt(n)))*A024924(floor(sqrt(n)))
#   a(n) = Sum_{k=1..m} (A008472(k)*A002321(floor(n/k)) + A008683(k)*A024924(floor(n/k))) - A002321(m)*A024924(m), where m = floor(sqrt(n)).

# Where A024924(n) can be computed in sublinear time as (recursively, using the sum of primes function):
#   A024924(n) = Sum_{k=1..floor(sqrt(n))} (A061397(k)*floor(n/k) + A034387(floor(n/k))) - A034387(floor(sqrt(n)))*floor(sqrt(n))

# See also:
#   https://oeis.org/A024924
#   https://oeis.org/A034387

func sum_of_sopf(n) {
    dirichlet_sum(n,
        { .is_prime ? _ : 0 },
        { 1 },
        { .sum_primes },    # FIXME: remove the recursive definition
        { _ }
    )
}

func sum_of_primes(n) {
    dirichlet_sum(n,
        { .sopf },
        { .mu },
        (sum_of_sopf),
        { .mertens }
    )
}

func A137851_partial_sum(n) {
    dirichlet_sum(n,
        {.is_prime ? _ : 0},
        {.mu},
        {.sum_primes},
        {.mertens}
    )
}

func sum_of_primes_2(n) {
    dirichlet_sum(n,
        {|k| k.prime_divisors.sum {|p| p * mu(k/p) } },
        { 1 },
        (A137851_partial_sum),
        { _ },
    )
}

say sum_primes(1e5)         #=> 454396537
say sum_of_primes(1e5)      #=> 454396537
say sum_of_primes_2(1e5)    #=> 454396537