Really large prime sum

Solution using Miller Rabin algorithm for primality testing:

 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
44
45
46
47
48
49
50
51

#include <iostream>
#include <cmath>
using namespace std;
typedef __int128 lli;

lli powMod(lli b, lli e, lli m){
    lli ans = 1;
    b %= m;
    while(e){
        if(e & 1) ans = ans * b % m;
        e >>= 1;
        if(e) b = b * b %m;
    }
    return ans;
}

bool test(lli n){
    if(n < 2) return false;
    if(n == 2) return true;
    lli d = n - 1, s = 0;
    while(!(d & 1)){
        d >>= 1;
        ++s;
    }
    for(int i = 0; i < 16; ++i){
        lli a = 1 + rand() % (n - 1);
        lli m = powMod(a, d, n);
        if(m == 1 || m == n - 1) goto exit;
        for(int k = 0; k < s - 1; ++k){
            m = m * m % n;
            if(m == n - 1) goto exit;
        }
        return false;
        exit:;
    }
    return true;
}

int main(){
    ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    long long L, mod = 10007;
    cin >> L;
    long long sum = 0;
    for(long long i = 2, u = sqrt(L); i <= u; i+=2){
        long long num = i*i+1;
        if(test(num)) sum += num%mod;
        if(sum>=mod) sum -= mod;
    }
    cout << sum << "\n";
    return 0;
}

Testing $10^{12}$ numbers for primality and then adding them up is just not feasible so we need to be smart about how we are testing. The first thing we need to do is realize that we only need to test $10^{6}$ numbers by testing all values of $n^2 + 1$ for $n$ in the interval $[1, 10^6]$ . The second thing is to realize that if you are testing if a number $n$ divides $m$ we only need to test the values of $n$ on the interval $[0, \sqrt{m}]$ . This will make the problem at least feasible, but we can cut off a lot more time in it. Other than 2 there won't be any other even prime; for an integer $n > 1$ , $n^2 + 1$ can only be prime if $n$ is even, so we only need to try the even values of $n$ which shrinks our search by half. The same thing for the divisors of a number, if we have an odd number $m$ we only need to try if odd numbers divide it, so the code looks like the code below. Another good solution that I did not implement but you are welcome to try is using the Pollard's rho algorithm.

 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

#include <bits/stdc++.h>

using namespace std;
typedef unsigned long long ull;

bool isPrime(const ull &x){
    ull sqrt_x = round(sqrt(x));
    for(ull i=3; i<= sqrt_x; i+=2){
        if(x%i == 0)
            return false;
    }
    return true;
}

int main(){
    ull sum = 0;

    for(ull i=2; i < 1000000; i = i+2){
        if(isPrime(i*i + 1)){
            sum += (i*i + 1);
            sum = sum % (10007);
        }
        if(i%10000 == 0)
            printf("%llu%%\n", (100*i)/1000000);
    }
    printf("%llu\n", sum);
    return 0;
}

$(\text{Plus}\text{@@}\text{Table}\left[\text{If}\left[\text{PrimeQ}\left[i^2+1\right],i^2+1,\text{Nothing}\right],\{i,2,1000000\}\right] \bmod 10007) \Longrightarrow 3601$