Today I did a lot of things during the day. But many of them did not finish…. at least for now.

Hopefully, I was able to finish this problem today. Not much to say, it’s just a brute force approach.

Project Euler 263

The problem requires finding numbers with properties of primes in an arithmetic progression and being also practical:

Primes in arithmetic progression

Practical Number

So, let’s start thinking a bit about the requirement for a practical number:

Suppose that, $n$ is a practical number, then it’s possible to write each number from $1$ to $n$ as the sum of unique divisors of $n$.

1: It's the trivial case, because $n$ will always have divisor $1$.

2: It cannot be written in any other way than only $2$, so $n \equiv 0 \mod 2$.

3: It's possible to write $3$ as $1+2$, or just like $3$. But, the problem requires also that $[n-9, n-3, n+3, n+9]$ are all prime numbers. So $n \not\equiv 0 \mod 3$.

4: It can't be written as other way than just $4$, so $n \equiv 0 \mod 4$.

We can continue with it, but we have enough conditions for now.

$$n \equiv 4 \mod 12 \Vert n \equiv 8 \mod 12$$

Additionally, you can see that by the condition of primality that $n \equiv 0 \mod 5$. Then:

$$n \equiv 20 \mod 60 \Vert n \equiv 40 \mod 60$$

Practical numbers can be checked with the following equation:

$$p_{i} \leq 1 + \sigma(p_{1}^{\alpha_{1}}p_{2}^{\alpha_{2}}\ldots p_{i-1}^{\alpha_{i-1}}) = 1 + \prod_{j=1}^{i-1} \frac{p_{j}^{\alpha_{j}+1} - 1}{p_{j}-1}$$

I made a code for this approach in C++, but the numbers are really big so I’m having memory issues. So I had to switch to using Mathematica instead, in both cases, the bottleneck is the generation of primes.

For example, taking 10000000 primes and applying our modulus condition, we reduce the size to 1250278.

Applying the primality test for the next numbers, reduce the size to 7324.

Applying the condition of successive primes, reduce the size to 5412.

However, you need to use even more primes than just $10^{7}$ to reach the first answer.

Here is the code in C++, modifying the sieve to a segmented sieve should be enough, but I don’t want to waste more time on it.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
typedef pair<ll,ll> pll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef vector<pii> vii;
typedef vector<pll> vll;
#define N 230000005

int lpf[N];
vi primes;
int pi[N];
vi possible_ans;

ll fastexp(ll x, ll y){
    ll ans = 1;
    while(y > 0){
        if(y&1) ans = (ans*x);
        y = y>>1;
        x = (x*x);
    }
    return ans;
}

void sieve(){
    for(int i=2; i<N; i++){
        if(!lpf[i]){
            primes.push_back(i);
            lpf[i] = i;
            pi[i] = pi[i-1]+1;
            for(ll j=1LL*i*i; j<N; j+=i){
                if(lpf[j] == 0) lpf[j] = i;
            }
        }
        else pi[i] = pi[i-1];
    }
    for(auto p: primes){
        if(p%60 == 11 or p%60 == 31){
            if(lpf[p+6] == p+6 and lpf[p+12] == p+12 and lpf[p+18] == p+18){
                if(pi[p+6] - pi[p] == 1 and pi[p+12] - pi[p+6] == 1 and pi[p+18] - pi[p+12] == 1){
                    possible_ans.push_back(p+9);
                }
            }
        }
    }
}

vll factorize(ll n){
    vll ans;
    while(n>1){
        int d = lpf[n];
        int c = 0;
        while(n%d == 0){
            c++;
            n /= d;
        }
        ans.push_back(make_pair(d, c));
    }
    return ans;
}

bool is_practical(ll n){
    vll factors = factorize(n);
    ll current_sum = (fastexp(factors[0].first, factors[0].second+1)-1)/(factors[0].first-1);
    for(int i=1; i<factors.size(); i++){
        int p = factors[i].first;
        int a = factors[i].second;
        if(p > 1 + current_sum){
            return false;
            break;
        }
        current_sum *= (fastexp(p, a+1) - 1)/(p-1);
    }
    return true;
}

int main(){
    sieve();
    cout << possible_ans.size() << '\n';
    for(auto pr: possible_ans){
        ll value = pr;
        bool take = true;
        for(int i=-8; i<=8; i+=4){
            if(!is_practical(value+i)){
                take = false;
                break;
            }
        }
        if(take){
            cout << value << '\n';
        }
    }
    return 0;
}

The code can find the first answer on my computer: 219869980.