123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114 |
- #!/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
|