Hexy prime sum

1. Generate all primes less than $h_{2500} = 12497500$ using Sieve of Eratosthenes .
2. Build table prime[i] to store the largest prime smaller than $i$ .
3. Sum prime[k] for all $k=h_i$ for $i=2\dots 2500$ .

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

#include <iostream>
#include <bitset>
using namespace std;

const int MAXN = 12497505;

int prime[MAXN];
bitset<MAXN> pt;

int main() {

    // sieve
    pt.set(); pt[0] = pt[1] = 0;
    for (int i = 2; i < MAXN; i++) if (pt[i] == 1)
        for (int j = i + i; j < MAXN; j += i)
            pt[j] = 0;

    int curr = 0;
    for (int i = 2; i < MAXN; i++) {
        prime[i] = curr;

        if (pt[i] == 1)
            curr = i;
    }

    long long ans = 0;
    for (int i = 2; i <= 2500; i++) {
        int k = 2 * i * i - i;

        ans += prime[k];
    }

    cout << ans << endl;

    return 0;
}

An alternative approach is to use a fast primality test algorithm called Miller-Kabin test . The code in the wiki is quite inefficient. Below is an improved version mainly on faster computation of $a^d \mod n$ .

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

from random import randint

def pow(a,d,n):
    ans = 1
    while d != 0:
        if d & 1:    # d is odd
            ans = (ans * a) % n
        d = d / 2
        a = (a * a) % n
    return ans

def isPrime(n, k):
    if n < 2: return False
    if n < 4: return True

    # speedup test
    if n % 2 == 0: return False

    s = 0
    d = n - 1
    while d & 1 == 0:    # d is even
        s += 1
        d /= 2

    for i in range(k):
        a = randint(2, n-2)    # 2 <= a <= n-2
        x = pow(a,d,n)
        if x == 1: continue
        for j in range(s):
            if x == n-1: break
            x = (x * x) % n
        else:
            return False
    return True

ans = 0
for i in range(2,2501):
    num = 2 * i * i - i
    while not isPrime(num,100):
        num -= 1
    ans += num

print ans

On my machine, both algorithm takes roughly the same runtime (4~5 sec). However, as $n$ increases, I believe the first approach will be less efficient. In addition, the second approach use much less memory than the first one. Miller-Kabin for the win!

another general algorithm is to first generate all primes less than $h_n$ which can be done quickly using a sieve. Then you simply start a pointer at the first prime number, namely 2, and for each $h_n$ you move the pointer forward until you reach the largest prime less than $h_n$ at which point you increment the running sum by that prime number. Once you've stepped through all the values for n the running sum will contain the answer.

I use a raw approach to evaluate $\zeta(n)$ but it works since $n < 2501$

Here is my code -

1
2
3
4
5
6
7
8
9

def zeta(n):
    n -= 1 if n % 2 == 0 else 2
    while not all(n % i for i in range(3, int(n ** 0.5) + 1, 2)):
        n -= 2
    return n
z = 0
for x in xrange(2, 2501):
    z += zeta(2 * x * x - x)
print z

Here's my take on the problem:So basically I used the Miller-Rabin primality test,the code is available here

1
2
3
4
5
6
7
8
9

sumseries=0
for i in range(2,2501):
    N=i*(2*i -1)
    #if  N is not prime then keep reducing it till you reach the nearest prime number
    while   not millr_rabin(N):
        N=N-1
    sumseries=sumseries+N

print sumseries

My method took about 5 minutes, but that's because my primality test was fairly slow. I calculated $h=2n^2-n$ , then ran a For Loop that counted down from $h$ by 1 until it reached a prime number, and then I added the number to the running total.

My primality test involved using a function I had written previously that gave a list of all the divisors of the number, and if that list had a length of $2$ (1 and itself) then it concluded it was prime. I feel that with a faster primality test this method could be fairly quick.