modified: Math/sum_of_prime_powers.sf -- added PARI/GP program.

Author: trizen

Math/sum_of_k-almost_primes.sf

@@ -21,7 +21,7 @@
 func almost_prime_sum(n,k) {
 
     if (k == 1) {
-        return prime_sum(n)
+        return n.prime_sum
     }
 
     var total = 0

Math/sum_of_k-omega_primes.sf

@@ -16,7 +16,7 @@
 func omega_prime_sum(n, k=1) {
 
     if (k == 1) {
-        return prime_powers(n).sum
+        return n.prime_power_sum
     }
 
     var total = 0

Math/sum_of_prime_powers.sf

@@ -2,6 +2,19 @@
 
 # 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 - 1); f(n) + g(sqrtint(n)) - g(sqrtnint(n, 3)) + 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;
+
+)
+
 func sum_of_primes(n, j=1) {
 
     return 0 if (n <= 1)