Моё решение задачи 146
Необходимо найти сумму всех натуральных \(n\), что \(n^2+1\), \(n^2+3\), \(n^2+7\), \(n^2+9\), \(n^2+13\), и \(n^2+27\) будут последовательными простыми числами.
Полное условие можно найти тут
Хочется отметить, что сложность у задачи 50%, а на текущий момент её решило меньше 4000 человек. Тем не менее, мне она показалось простой. Простейшее решение отработало очень быстро.
Для начала, можно отметить, что в лоб проверять условие очень долго. Проверять на простоту числа порядка \(10^{15}\) достаточно сложно, поэтому их нужно как-то отсеять.
Самое простое — не рассматривать те \(n\), что хотя бы одно из \(n^2+1\), \(n^2+3\), \(n^2+7\), \(n^2+9\), \(n^2+13\), и \(n^2+27\) будет заведомо делиться на какое-то маленькое простое число. Это даёт достаточно хорошие результаты: из 150 миллионов чисел, после отсеивания по простым числам \(< 3000\) (этот параметр я подбирал уже после решения задач: если он слишком маленький, то будет слишком много проверок на простоту, если же слишком большой, то мы делаем слишком много работы, чтобы отсеять несколько чисел), останется меньше \(2000\) чисел. Их уже можно проверить непосредственно.
Тогда алгоритм может быть таким:
- Находим простые числа меньше \(3000\).
- Для каждого из них находим допустимые остатки.
- Для каждого из чисел от \(1\) до \(n\) проверяем, что остатки по всем простым хорошие.
- Непосредственно проверяем условие. Важно не забыть проверить непростоту оставшихся нечётных чисел из диапазона \(n^2 + 1 \ldots n^2 + 27\) там могут быть (и будут!) другие простые числа.
Непосредственно сам поиск такой клики можно реализовать тривиально. Ниже мой код на C++11 с использованием библиотек Flint и primesieve. Распараллеливание хоть и просится, но смысла не имеет, т.к. я получил ответ менее, чем за 5 секунд.
/*
* Problem 146 on Project Euler
* Aleksey Lobanov (c) 2016
*/
#include <iostream>
#include <vector>
#include <cstdint>
#include <set>
#include <iomanip>
#include <algorithm>
#include "fmpzxx.h"
#include "arithxx.h"
#include "primesieve.hpp"
using namespace std;
using namespace flint;
bool is_prime(int64_t num) {
fmpz_factorxx fact;
fact.set_factor(num);
if ( fact.size() != 1 )
return false;
if ( fact.exp(0) != 1)
return false;
return true;
}
bool is_possible(int64_t num, const vector<int64_t> &to_add) {
for (auto &&add: to_add)
if ( !is_prime(num*num + add) )
return false;
// primes must be consecutive
// so we need check, that other numbers like n^2 + i is not primes
vector<int64_t> other_adds;
for (size_t i = to_add[0] + 2; i < to_add[to_add.size() - 1]; i += 2)
if ( !binary_search(to_add.begin(), to_add.end(), i) )
other_adds.push_back(i);
for (auto &&add: other_adds)
if ( is_prime(num*num + add) )
return false;
return true;
}
int main() {
const vector<int64_t> to_add = {1, 3, 7, 9, 13, 27};
const int64_t MAX_N = 1l*150*1000*1000;
const int64_t MAX_PRIME = 3000;
vector<int64_t> sieve_primes;
primesieve::generate_primes(MAX_PRIME, &sieve_primes);
vector< vector <int64_t> > good_remainders;
for (auto &&prime: sieve_primes) {
set<int64_t> remainders;
for(int64_t rem = 0; rem < prime; ++rem)
remainders.insert(rem);
set<int64_t> base_remainders;
for (auto &&base: to_add)
base_remainders.insert(base % prime);
for (int64_t rem = 0; rem < prime; ++rem)
for (auto &&base_rem: base_remainders)
if ( (rem*rem + base_rem) % prime == 0 )
remainders.erase(rem);
good_remainders.push_back(vector<int64_t>(remainders.begin(), remainders.end()));
}
size_t cnt = 0;
size_t sum = 0;
/* WARNING
* for small n can be that
* n^2 + [1 or 3 or .. ] is prime in sieve primes
* but there is only one n < 315410 is 10,
* so we need add 10 to n
*/
sum += 10;
for (int64_t i = 1; i < MAX_N; ++i) {
bool is_good = true;
for (size_t prime_ind = 0; prime_ind < sieve_primes.size(); ++prime_ind) {
if ( !binary_search(good_remainders[prime_ind].begin(), good_remainders[prime_ind].end(), i % sieve_primes[prime_ind]) ) {
is_good = false;
break;
}
}
if ( is_good ) {
cnt++;
if ( is_possible(i, to_add) ) {
sum += i;
}
}
}
cout << "count = " << cnt << endl;
cout << "Result is: " << sum << endl;
return 0;
}
Ответ: 676333270