Некий китаец посчитал количество решений до n=59 (oeis.org/A071983), и у меня нет идей, как он в принципе мог это сделать. Там похоже, что количество решений растет как 4^n, поэтому любое перечисление решений будет неэффективно. Тем не менее, простой DFS с оптимизациями (представление ребер графа и подмножеств вершин битовыми масками) на практике оказался почти самым быстрым: 4 часа для n=46 (во всех алгоритмах я считаю однопоточное исполнение). До этого его скорость растет примерно как 3^n, но понятно, что на больших n она выростет как минимум до 4^n.
DP оказалось неэффективно: для n=39 у меня уже не хватает памяти (растет экспоненциально), а решение для n=38 занимает 4.5 минуты.
Самое эффективное решение, которое я нашел, - комбинация DFS и DP: для нечетного n ищем с помощью DFS количество цепочек F(a, S), начинающихся в a, для всех подмножнств S мощности (n+1)/2. Эти позволяет рассчитать количество гамильтоновых путей, имеющих a ровно посередине. Ищет n=47 за 6 часов, до этого растет примерно как (1.7)^n. К сожалению, память растет тоже экспоненциально (но медленнее, чем стандартный DP).
Самый лучший теоретический алгоритм основан на принципе inclusion/exclusion, он имеет сложность (2^n)*(n^3) при полиномиальной памяти, но он требует возведения матриц nxn в степень n, поэтому он на практике слишком медленный (реализация на питоне не доходит даже до n=30). Рано или поздно он обгонит по скорости все прочие алгоритмы, но, скорее всего, время работы будет превышать время существование вселенной. Для этого алгоритма, правда, есть огромное поле оптимизаций (полное распараллеливание на 2^n умножений матриц, которые можно делать на GPU) - но по моим подсчетам, все равно будет слишком медленно, если только не найти эвристик, которые позволят не считать большую часть подмножеств.
А почему у Вас получилась такая оценка? всего вариантодля еще дают , и для построения каждого надо примерно операций (примерное количество ребер из ). То есть, получается И кажется, можно обойтись без - если взять - количество гамильтоновых путей в , начинающихся в .
Попытки сделать криптовалюту на простых числах были, например, gapcoin. Они ищут наибольшую разницу между ближайшими простыми числами (точнее, prime gap merit), и утверждают, что их proof of work - не просто сжигание электричества, но приносит пользу науке. Но, похоже, не прижилось
Томаш Рокицки уже провел некоторое исследование и OEIS-тификацию
Тут хочется добавить, что это тот самый Tomas Rokicki, главный гуру в recreational computer science, автор программы Golly (самая известная программа по симуляции игры "жизнь" Конвея); также, он доказал, что кубик Рубика собирается за не больше, чем 20 ходов из любой позиции (они использовали мощности Гугла, и в сумме затратили 35 cpu-лет). Так что, если такой человек заинтересовался этой задачей, то уже можно считать, что она имеет ненулевое значение в научной и околонаучной тусовке.
Я запустил ваш код на моей RTX 4050 на ноутбуке: 12 -> 126.820. Если экстраполировать, то N=16 он обработает за 14 дней. А на какой видеокарте вы запускали?
Я придумал еще несколько оптимизаций, на моем ноутбуке в 4 раза быстрее, чем ваша версия.
Чтобы a было в нужном диапазоне, надо rev(b) < b, и можно рассматривать только те b, у которых последняя цифра не больше первой (а также не 0).
Как вы доказали, для каждого из остатков b%99 есть 0 или 1 остаток a%99. Их можно предрассчитать и сохранить в static storage в видеокарте. У трети значений остатка нет, их можно пропускать.
я использую int128 для вычисления точного диапазона, где может лежать a. Если нужный остаток a%99 не попадает в этот диапазон, то вычисления можно пропустить, и переходить к следующему b. Оказывается, только примерно каждое 75-е число попадает в этот диапазон.
Каждый поток проверяет числа в диапазоне, кратном 9900. Таким образом, и последняя цифра во всех потоках блока одинаковая, и остаток от деления на 99 тоже, поэтому и conditional flow у них одинаковый (кроме предыдущего пункта, когда надо проверять само условие ab=10^N*rev(b)+rev(a), что случается редко). И доступ к массиву остатков происходит по одному адресу для всех потоков (multicast), что эффективно.
rev(b), b%99 считаются инкрементально, что быстрее, чем их вычислять для каждого b
rev считается в два приема по предвычисленным кешам для чисел длины N/2, которые хранятся в global memory. Вроде бы он не так часто выполняется, и пока потоки ждут ответа от памяти, другие потоки могут выполнять вычисления. Но я не профилировал, поэтому не уверен.
Если интересно, вот код, для простоты поддержал только четные N
Hidden text
#include <vector>
#include <iostream>
#include <cuda.h>
#include <cuda_runtime.h>
using namespace std;
void ERR() {
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
cerr << cudaGetErrorString(err) << endl;
throw runtime_error(cudaGetErrorString(err));
}
}
static constexpr int N = 12;
static_assert(N % 2 == 0);
static constexpr int64_t Exp10(int n) {
return n == 0 ? 1 : 10 * Exp10(n - 1);
}
static_assert(Exp10(6) == 1000000);
static_assert(Exp10(12) == 1000000000000ll);
static constexpr int PART_SIZE = 99000;
static_assert(PART_SIZE % 9900 == 0);
static constexpr int rangeMax(int b0) {
return (10 + b0 - 1) / b0;
}
vector<int> vect_remainders_map {0, -1, 47, 24, -1, 89, 84, -1, 59, 9, -1, 11, 60, -1, 35, 30, -1, 95, 72, -1, 20, -1, -1, 71, 3, -1, 41, 45, -1, 83, 15, -1, -1, 66, -1, 14, 90, -1, 56, 51, -1, 26, 75, -1, 77, 27, -1, 2, 96, -1, 62, 39, -1, 86, -1, -1, 38, 69, -1, 8, 12, -1, 50, 81, -1, -1, 33, -1, 80, 57, -1, 23, 18, -1, 92, 42, -1, 44, 93, -1, 68, 63, -1, 29, 6, -1, 53, -1, -1, 5, 36, -1, 74, 78, -1, 17, 48, -1, -1,};
__constant__ int remainders_map[99];
static constexpr int64_t HALF = Exp10(N / 2);
static int REVERSE[HALF];
int* reverse_mem;
int64_t make_reverse(int64_t z, int n) {
int64_t rev = 0;
for (int i = 0; i < n; i++) {
int d = z % 10;
rev = rev * 10 + d;
z /= 10;
}
return rev;
}
__device__ static constexpr int64_t MIN = Exp10(N - 1);
__device__ static constexpr int64_t MAX = Exp10(N);
__device__ int64_t reverse(int64_t n) {
int64_t rev = 0;
for (int i = 0; i < N; i++) {
int64_t next = n / 10;
int d = n - next * 10;
rev = rev * 10 + d;
n = next;
}
return rev;
}
__device__ int64_t reverse2(int64_t n, int* mem) {
return
mem[n % HALF] * HALF +
mem[n / HALF];
}
__device__ void found(int64_t a, int64_t b) {
printf("%lld %lld\n", a, b);
}
__device__ void process_100(int b0, int range_max, int64_t b_base, int* mem) {
constexpr int range_min = 1;
int64_t inv_b_base = reverse2(b_base, mem);
int64_t inv_b = inv_b_base;
int64_t b_rem99_base = b_base % 99;
int64_t b_rem99 = b_rem99_base;
for (int b_second = 0; b_second < 100; b_second += 10) {
for (int b_first = 1; b_first <= b0; b_first++) {
int64_t b = b_base + b_second + b_first;
b_rem99++;
if (b_rem99 == 99) b_rem99 = 0;
inv_b += MIN;
if (inv_b > b) continue;
auto a_rem99 = remainders_map[b_rem99];
if (a_rem99 == -1) continue;
auto R = (int64_t)((__int128_t)MAX * inv_b / b);
int64_t R_rem99 = R % 99;
int64_t rem_diff = a_rem99 < R_rem99 ? 99 + a_rem99 - R_rem99 : a_rem99 - R_rem99;
if (rem_diff < range_min || rem_diff > range_max) continue;
int64_t a = R + rem_diff;
auto mul = (__int128_t)a * b;
int64_t inv_a = reverse2(a, mem);
if (mul == (__int128_t)MAX * inv_b + inv_a) {
found(a, b);
}
}
inv_b_base += (MIN / 10);
inv_b = inv_b_base;
b_rem99_base += 10;
if (b_rem99_base >= 99) b_rem99_base -= 99;
b_rem99 = b_rem99_base;
}
}
__device__ void process_9900(int b0, int range_max, int64_t b_base, int64_t b_max, int* mem) {
for (int i = 0; i < PART_SIZE; i += 100) {
if (b_base + i >= b_max) break;
process_100(b0, range_max, b_base + i, mem);
}
}
__global__ void process(int b0, int range_max, int* mem) {
int64_t b_min = b0 * MIN;
int64_t b_max = b_min + MIN;
auto tid = (int64_t)blockIdx.x * blockDim.x + threadIdx.x;
int64_t b_base = b_min + tid * PART_SIZE;
if (b_base >= b_max) return;
process_9900(b0, range_max, b_base, b_max, mem);
}
void main_process(int b0) {
cerr << "processing first digit: " << b0 << endl;
int64_t count = (MIN + PART_SIZE - 1) / PART_SIZE;
int threadsPerBlock = 256;
int64_t blocksPerGrid = (count + threadsPerBlock - 1) / threadsPerBlock;
process<<<blocksPerGrid, threadsPerBlock>>>(b0, rangeMax(b0), reverse_mem);
ERR();
cudaDeviceSynchronize();
ERR();
}
int main(int argc, char *argv[]) {
for (int64_t i = 0; i < HALF; i++) {
REVERSE[i] = make_reverse(i, N / 2);
}
cerr << "N = " << N << endl;
cudaMemcpyToSymbol(remainders_map, &vect_remainders_map[0], vect_remainders_map.size() * sizeof(int));
ERR();
cudaMalloc(&reverse_mem, HALF * sizeof(int));
ERR();
cudaMemcpy(reverse_mem, REVERSE, HALF * sizeof(int), ::cudaMemcpyHostToDevice);
ERR();
for (int b0 = 1; b0 <= 9; b0++) {
main_process(b0);
}
return 0;
}
А можете рассказать, как Вы ускорили? можно код посмотреть? Я пытался написать код на c++ CUDA, но у меня медленнее получилось (но я не double использовал, а 128-битные числа)
я придумал усовершенствованный алгоритм, который позволил получить значение для n=6: 4505706 за 62 минуты
Можно сделать быстрее, если перебирать все пары, но более эффективно считать количество цифр в z=x*y.
Количество цифр в числе представить в виде одного 40-битного числа, где каждые 4 бита отведены под количество нулей, единиц итд.
Заранее рассчитать эти значения для чисел от 0 до N, и проверять counts[x] + counts[y] == counts[z mod 10^N] + counts[z / 10^N]
Для n=6 на моем ноутбуке задача считается за 11 минут (в одном потоке).
Можно сделать еще быстрее, если вместо z = x*y считать z_hi = x*y/10^N, z_lo=x*y mod 10^N, и воспользоваться тем, что x*(y+1) = x*y + x, чтобы обновлять значения в цикле, а не пересчитывать заново. Это избавляет от модуля и деления, и снижает время до 7 минут.
Доступ к массиву counts получается cache friendly для x, y, z_hi, но не для z_lo. Его можно разделить на z_lo_hi, z_lo_lo, каждый из которых будет в диапазоне до 10^((N+1)/2) - и их тоже обновлять в цикле. Так как диапазон гораздо меньше, он поместится в кеш процессора.
А я вот хочу автору сказать спасибо за статью, и позволю не согласиться, что статья ни о чем. Что я из нее понял:
Автор рассматривает обобщение гипотезы Коллатца на произвольный нечетный множитель и рассмаривает, наверное, самую простую задачу (не считая тривиальной задачи нахождения неподвижных точек) нахождения циклов длины 2. Я быстрым гуглопоиском не нашел, чтобы кто-то эту задачу до него сформулировал. И даже эта задача, по-видимому, не решается аналитически. Ну и наличие решений для и само по себе неожиданно.
Автор применил интересный трюк, показав, что . Тем самым он ускорил алгоритм по нахождению циклов с началом, меньшим , с у алгоритма "в лоб", до . Это позволило на питоне проверить все числа до , что позволяет высказать гипотезу, что циклов длины 2, скорее всего, больше нет.
Еще мне было интересно ускорить сам алгоритм, избавившись от извлечения корня и умножений с делениями (мне не удалось избавиться только от одного ). Мой вариант на с++ с первой попавшейся из интернета имплементацией bigint работает примерно в 10 раз быстрее версии автора на питоне.
Как к задаче подступиться аналитически, совершенно непонятно. Кажется, что (фигурные скобки - это дробная часть), и вероятность найти стремится к нулю, если в двоичном представлении корня из двух нет очень больших последовательностей из нулей, и . Это предположение верно, если корень из двух - это нормальное число в двоичной системе счисления, но это не доказано.
В общем, если то, что сделал автор, аккуратно сформулировать и красиво оформить, может, могла бы получиться интересная научная статья.
А Вы не сравнивали с их реализацией? они дают на нее в статье ссылку. Правда, их статья десятилетней давности, а код последний раз обновлялся 5 лет назад. Ну и всегда можно напрямую автору написать, у Линдена есть профиль на линкедине
Про SIMD можно найти прямо вашу задачу на StackOverflow. Хотя, кажется, можно проще, чем там предлагают, если использовать __mm_shuffle_epi8. Но не уверен, что это сильно поможет - возможно, вы уже приблизились к границе пропускной способности памяти
А сколько времени в итоге занимает обработка тестового файла?
Распаковка и чтение с диска не будет занимать больше времени, чем, собственно, сам подсчет?
Спасибо за статью, я раньше думал, что memory-mapped файлы не очень хорошо работают на Windows.
База на диске — хороший компромисс между скоростью работы и простотой реализации, но для некоторых задач она не нужна:
— Если у вас в основном последовательный доступ, то есть, вам надо перебрать простые числа на интервале [X… Y], то имеет смысл их каждый раз заново просеивать: достаточно заранее в памяти просчитать и держать список простых чисел до sqrt(N), в вашем случае — до миллиона.
— Если у вас полностью рандомный доступ, то выгоднее воспользоваться критерием простоты чисел: en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases. Аккуратно написанный, он дает скорость 20 микросекунд на число (в отличие от 2 миллисекунд на seek на HDD): ceur-ws.org/Vol-1326/020-Forisek.pdf
Альтернативно, базу можно сократить в два раза, если хранить биты только на каждое нечетное число, или даже в три раза, если хранить только числа, дающие при делении на 6 остатки 1 или 5 (остальные не могут быть простыми, кроме чисел 2 и 3). А можно вообще упороться, и хранить только числа, дающие один из восьми остатков от деления на 30, что как раз занимает один байт. Я это все описывал в habr.com/ru/post/526924, там есть код на C#, который можно использовать для этого.
вот тут я вижу как достигают скорости 2.5 миллиона evaluated cells per second на CPU
все-таки там 2.5 миллиарда клеток в секунду, а не миллиона.
А с использованием SSE-команд можно еще быстрее, до 15 миллиардов в секунду: habr.com/ru/post/505606
А в этой статье показывают, что на GPU можно достичь скорости, если я не ошибся в подсчетах, полтора триллиона клеток в секунду!
примерно после девяти миллиардов итераций и двенадцати часов работы алгоритм нашел score=8 ݻටຈݹ௫ටຯѧ௨پଈݶಏقງݹலටະऄ௨پໃݢ௨ගɑݸҰقະࡨ௨sໃɥ௨ڏІݟ௨ගɱݸతටࢩݹࢴටแɕ௧ڄໃףۑඤÐݸలقลڝ௨ටໃݹ௨ටࠋ
Отличная задача на поиск оптимальной стратегии! Основная сложность — для каждой позиции определение следующей фигурки требует перебора всех фигурок по всем смещениям и вращениям, а количество позиций растет экспоненциально с увеличением глубины колодца. Используя очень эффеткивный алгоритм (перебирающий около двухсот тысяч позиций в секунду), мне удалось решить задачу для глубины 5 (score=4, всего пришлось перебрать около четырёхсот миллионов позиций за 30 минут): ޛටדݹ௨ඞТஈ௧ƣເঅ௩ɦܘݶహइໂђ௨ටໞݹஜɦܘݚذචʅݹਐටࢺݹસقໃݹ௨ටໃݹਆ
И я в процессе решения задачи для глубины 6. Пока что он нашел score=6 (поле перебора трех миллиардов позиций) ܮටลຍ௨ಀඍݸహටາɝɷචŒݹமϺາŦ௨ಀІݟௐقГঅஶѼໃӌ௨ගɱݸƓΟເח௧ڈໃӻ௨ගȣݸಏටງђஶɈժށѻටໃݹ௨ටເ0
но алгоритм еще не закончил работу, пока перебрав восемь миллиардов позиций за 11 часов.
Я попробовал еще больше соптимизировать однопоточную версию, используя вместо буфера на 15 чисел три буфера на 100 чисел, пользуясь тем, что остаток от деления на 15 у чисел n и n + 300 совпадают. Поэтому, если есть буфер чисел, например, [10000..10099], то от него можно перейти к буферу с числами [10300..10399] только поменяв цифру 0 на 3 в позициях, соответствующих числам, не делящимся на 3 и 5. Потом от нее можно перейти к [10600...10699] и т.д. При переполнении младшего разряда предыдущий увеличивается на единичку (получился алгоритм сложения с числом 300 в столбик). Например, переход от [19900..19999] к [20200..20299] получается заменой 9 -> 2 на третей позиции, 9 -> 0 на второй, и 1 -> 2 на первой для всех чисел.
На моем компьютере (у меня Windows/cygwin, поэтому мне сложно сравнить с вашими результатами напрямую) этот алгоритм работает 1.7 секунд. Я Вам послал пул реквест, было бы интересно узнать, будет ли он обгонять reusebuffer в ваших тестах. Правда, сейчас он работает только для LIMIT = 10^k, но, думаю, его можно доделать для любого диапазона чисел без ухудшения производительности.
decimal.c
#include <stdio.h>
#include <string.h>
#define LIMIT 1000000000
char *FIZZ = "Fizz";
char *BUZZ = "Buzz";
struct Segment {
char buffer[4096]; // 100 numbers
int buffer_length;
int positions[100]; // offsets in buffer1 of third digit of all numbers
int positions_length; // number of items in positions
};
int get_number_length(int number) {
int number_length = 1;
while (number >= 10) {
number /= 10;
number_length++;
}
return number_length;
}
char* print_number(int number, int numberLength, char* buf) {
for(int i = 0; i < numberLength; i++) {
int digit = number % 10;
buf[numberLength - i - 1] = (digit + '0');
number /= 10;
}
buf[numberLength] = '\n';
return buf + numberLength + 1;
}
char* print_fizz(char* buf) {
*(int*)buf = *(int*)FIZZ;
buf[4] = '\n';
return buf + 5;
}
char *print_buzz(char *buf) {
*(int *)buf = *(int *)BUZZ;
buf[4] = '\n';
return buf + 5;
}
char *print_fizz_buzz(char *buf) {
*(int *)buf = *(int *)FIZZ;
*(int *)(buf + 4) = *(int *)BUZZ;
buf[8] = '\n';
return buf + 9;
}
char *print(int num, int num_length, char *ptr) {
if (num_length == -1)
num_length = get_number_length(num);
if (num % 15 == 0)
ptr = print_fizz_buzz(ptr);
else if (num % 3 == 0)
ptr = print_fizz(ptr);
else if (num % 5 == 0)
ptr = print_buzz(ptr);
else
ptr = print_number(num, num_length, ptr);
return ptr;
}
char *print_first_segment(char *ptr) {
for(int num = 1; num < 100; num++)
ptr = print(num, -1, ptr);
return ptr;
}
char *print_last_number(char *ptr) {
return print(LIMIT, -1, ptr);
}
void print_initial_segment(int first_number, struct Segment* segment) {
int digits = get_number_length(first_number);
segment->positions_length = 0;
char *ptr = segment->buffer;
for (int i = first_number; i < first_number + 100; i++) {
ptr = print(i, digits, ptr);
if (i % 3 != 0 && i % 5 != 0)
segment->positions[segment->positions_length++] = ptr - segment->buffer - 4;
}
segment->buffer_length = ptr - segment->buffer;
}
void print_initial_segments(int segment, struct Segment* s1, struct Segment* s2, struct Segment* s3) {
print_initial_segment(segment, s1);
print_initial_segment(segment + 100, s2);
print_initial_segment(segment + 200, s3);
}
void next(struct Segment* s) {
char c = s->buffer[s->positions[0]];
char nextC = c + 3;
if(c < '7') {
for (int i = 0; i < s->positions_length; i++)
s->buffer[s->positions[i]] = nextC;
} else if (s->buffer[s->positions[0] - 1] < '9') {
nextC -= 10;
int nextD = s->buffer[s->positions[0] - 1] + 1;
for (int i = 0; i < s->positions_length; i++) {
s->buffer[s->positions[i]] = nextC;
s->buffer[s->positions[i] - 1] = nextD;
}
} else {
nextC -= 10;
for (int i = 0; i < s->positions_length; i++) {
int pos = s->positions[i];
s->buffer[pos--] = nextC;
while(s->buffer[pos] == '9')
s->buffer[pos--] = '0';
s->buffer[pos]++;
}
}
}
int main(void) {
char* ptr;
struct Segment s1, s2, s3;
ptr = print_first_segment(s1.buffer);
fwrite(s1.buffer, 1, ptr - s1.buffer, stdout);
for (int segment = 100; segment < LIMIT; segment *= 10) {
print_initial_segments(segment, &s1, &s2, &s3);
fwrite(s1.buffer, 1, s1.buffer_length, stdout);
fwrite(s2.buffer, 1, s2.buffer_length, stdout);
fwrite(s3.buffer, 1, s3.buffer_length, stdout);
int repeat_count = 3 * (segment / 100);
for(int i = 1; i < repeat_count; i++) {
next(&s1);
next(&s2);
next(&s3);
fwrite(s1.buffer, 1, s1.buffer_length, stdout);
fwrite(s2.buffer, 1, s2.buffer_length, stdout);
fwrite(s3.buffer, 1, s3.buffer_length, stdout);
}
}
ptr = print_last_number(s1.buffer);
fwrite(s1.buffer, 1, ptr - s1.buffer, stdout);
return 0;
}
Некий китаец посчитал количество решений до n=59 (oeis.org/A071983), и у меня нет идей, как он в принципе мог это сделать.
Там похоже, что количество решений растет как 4^n, поэтому любое перечисление решений будет неэффективно.
Тем не менее, простой DFS с оптимизациями (представление ребер графа и подмножеств вершин битовыми масками) на практике оказался почти самым быстрым: 4 часа для n=46 (во всех алгоритмах я считаю однопоточное исполнение). До этого его скорость растет примерно как 3^n, но понятно, что на больших n она выростет как минимум до 4^n.
DP оказалось неэффективно: для n=39 у меня уже не хватает памяти (растет экспоненциально), а решение для n=38 занимает 4.5 минуты.
Самое эффективное решение, которое я нашел, - комбинация DFS и DP: для нечетного n ищем с помощью DFS количество цепочек F(a, S), начинающихся в a, для всех подмножнств S мощности (n+1)/2. Эти позволяет рассчитать количество гамильтоновых путей, имеющих a ровно посередине. Ищет n=47 за 6 часов, до этого растет примерно как (1.7)^n. К сожалению, память растет тоже экспоненциально (но медленнее, чем стандартный DP).
Самый лучший теоретический алгоритм основан на принципе inclusion/exclusion, он имеет сложность (2^n)*(n^3) при полиномиальной памяти, но он требует возведения матриц nxn в степень n, поэтому он на практике слишком медленный (реализация на питоне не доходит даже до n=30). Рано или поздно он обгонит по скорости все прочие алгоритмы, но, скорее всего, время работы будет превышать время существование вселенной.
Для этого алгоритма, правда, есть огромное поле оптимизаций (полное распараллеливание на 2^n умножений матриц, которые можно делать на GPU) - но по моим подсчетам, все равно будет слишком медленно, если только не найти эвристик, которые позволят не считать большую часть подмножеств.
А почему у Вас получилась такая оценка? всего
вариантодля 
еще 
дают 
, и для построения каждого 
надо примерно 
операций (примерное количество ребер из 
). То есть, получается 
- если взять 
- количество гамильтоновых путей в 
, начинающихся в 
.
И кажется, можно обойтись без
Попытки сделать криптовалюту на простых числах были, например, gapcoin.
Они ищут наибольшую разницу между ближайшими простыми числами (точнее, prime gap merit), и утверждают, что их proof of work - не просто сжигание электричества, но приносит пользу науке.
Но, похоже, не прижилось
Тут хочется добавить, что это тот самый Tomas Rokicki, главный гуру в recreational computer science, автор программы Golly (самая известная программа по симуляции игры "жизнь" Конвея);
также, он доказал, что кубик Рубика собирается за не больше, чем 20 ходов из любой позиции (они использовали мощности Гугла, и в сумме затратили 35 cpu-лет).
Так что, если такой человек заинтересовался этой задачей, то уже можно считать, что она имеет ненулевое значение в научной и околонаучной тусовке.
Я запустил ваш код на моей RTX 4050 на ноутбуке:
12 -> 126.820. Если экстраполировать, то N=16 он обработает за 14 дней. А на какой видеокарте вы запускали?Я придумал еще несколько оптимизаций, на моем ноутбуке в 4 раза быстрее, чем ваша версия.
Чтобы
aбыло в нужном диапазоне, надоrev(b) < b, и можно рассматривать только теb, у которых последняя цифра не больше первой (а также не 0).Как вы доказали, для каждого из остатков
b%99есть 0 или 1 остатокa%99. Их можно предрассчитать и сохранить в static storage в видеокарте. У трети значений остатка нет, их можно пропускать.я использую int128 для вычисления точного диапазона, где может лежать
a. Если нужный остатокa%99не попадает в этот диапазон, то вычисления можно пропустить, и переходить к следующемуb. Оказывается, только примерно каждое 75-е число попадает в этот диапазон.Каждый поток проверяет числа в диапазоне, кратном 9900. Таким образом, и последняя цифра во всех потоках блока одинаковая, и остаток от деления на 99 тоже, поэтому и conditional flow у них одинаковый (кроме предыдущего пункта, когда надо проверять само условие
ab=10^N*rev(b)+rev(a), что случается редко). И доступ к массиву остатков происходит по одному адресу для всех потоков (multicast), что эффективно.rev(b),b%99считаются инкрементально, что быстрее, чем их вычислять для каждогоbrev считается в два приема по предвычисленным кешам для чисел длины
N/2, которые хранятся в global memory. Вроде бы он не так часто выполняется, и пока потоки ждут ответа от памяти, другие потоки могут выполнять вычисления. Но я не профилировал, поэтому не уверен.Если интересно, вот код, для простоты поддержал только четные
NHidden text
А можете рассказать, как Вы ускорили? можно код посмотреть? Я пытался написать код на c++ CUDA, но у меня медленнее получилось (но я не double использовал, а 128-битные числа)
у меня не только при посте но и, например, при нажатии на "нравится". Рефреш страницы временно спасает
Спасибо, интересные задачи!
Можно сделать быстрее, если перебирать все пары, но более эффективно считать количество цифр в
z=x*y.Количество цифр в числе представить в виде одного 40-битного числа, где каждые 4 бита отведены под количество нулей, единиц итд.
Заранее рассчитать эти значения для чисел от 0 до
N, и проверятьcounts[x] + counts[y] == counts[z mod 10^N] + counts[z / 10^N]Для
n=6на моем ноутбуке задача считается за 11 минут (в одном потоке).Можно сделать еще быстрее, если вместо
z = x*yсчитатьz_hi = x*y/10^N, z_lo=x*y mod 10^N, и воспользоваться тем, чтоx*(y+1) = x*y + x, чтобы обновлять значения в цикле, а не пересчитывать заново. Это избавляет от модуля и деления, и снижает время до 7 минут.Доступ к массиву
countsполучается cache friendly дляx, y, z_hi, но не дляz_lo. Его можно разделить наz_lo_hi, z_lo_lo, каждый из которых будет в диапазоне до10^((N+1)/2)- и их тоже обновлять в цикле. Так как диапазон гораздо меньше, он поместится в кеш процессора.У меня этот вариант работает за 5 минут.
А я вот хочу автору сказать спасибо за статью, и позволю не согласиться, что статья ни о чем.
Что я из нее понял:
Автор рассматривает обобщение гипотезы Коллатца на произвольный нечетный множитель и рассмаривает, наверное, самую простую задачу (не считая тривиальной задачи нахождения неподвижных точек) нахождения циклов длины 2. Я быстрым гуглопоиском не нашел, чтобы кто-то эту задачу до него сформулировал. И даже эта задача, по-видимому, не решается аналитически. Ну и наличие решений для
и 
само по себе неожиданно.
Автор применил интересный трюк, показав, что
. Тем самым он ускорил алгоритм по нахождению циклов с началом, меньшим 
, с 
у алгоритма "в лоб", до 
. Это позволило на питоне проверить все числа до 
, что позволяет высказать гипотезу, что циклов длины 2, скорее всего, больше нет.
Еще мне было интересно ускорить сам алгоритм, избавившись от извлечения корня и умножений с делениями (мне не удалось избавиться только от одного
). Мой вариант на с++ с первой попавшейся из интернета имплементацией bigint работает примерно в 10 раз быстрее версии автора на питоне.
Как к задаче подступиться аналитически, совершенно непонятно. Кажется, что
(фигурные скобки - это дробная часть), и вероятность найти 
стремится к нулю, если в двоичном представлении корня из двух нет очень больших последовательностей из нулей, и 
. Это предположение верно, если корень из двух - это нормальное число в двоичной системе счисления, но это не доказано.
В общем, если то, что сделал автор, аккуратно сформулировать и красиво оформить, может, могла бы получиться интересная научная статья.
А Вы не сравнивали с их реализацией? они дают на нее в статье ссылку.
Правда, их статья десятилетней давности, а код последний раз обновлялся 5 лет назад. Ну и всегда можно напрямую автору написать, у Линдена есть профиль на линкедине
Про SIMD можно найти прямо вашу задачу на StackOverflow. Хотя, кажется, можно проще, чем там предлагают, если использовать __mm_shuffle_epi8.
Но не уверен, что это сильно поможет - возможно, вы уже приблизились к границе пропускной способности памяти
Распаковка и чтение с диска не будет занимать больше времени, чем, собственно, сам подсчет?
Доказано, что эта задача — NP-полная, en.wikipedia.org/wiki/Addition-chain_exponentiation:
База на диске — хороший компромисс между скоростью работы и простотой реализации, но для некоторых задач она не нужна:
— Если у вас в основном последовательный доступ, то есть, вам надо перебрать простые числа на интервале [X… Y], то имеет смысл их каждый раз заново просеивать: достаточно заранее в памяти просчитать и держать список простых чисел до sqrt(N), в вашем случае — до миллиона.
— Если у вас полностью рандомный доступ, то выгоднее воспользоваться критерием простоты чисел: en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases. Аккуратно написанный, он дает скорость 20 микросекунд на число (в отличие от 2 миллисекунд на seek на HDD): ceur-ws.org/Vol-1326/020-Forisek.pdf
Альтернативно, базу можно сократить в два раза, если хранить биты только на каждое нечетное число, или даже в три раза, если хранить только числа, дающие при делении на 6 остатки 1 или 5 (остальные не могут быть простыми, кроме чисел 2 и 3). А можно вообще упороться, и хранить только числа, дающие один из восьми остатков от деления на 30, что как раз занимает один байт. Я это все описывал в habr.com/ru/post/526924, там есть код на C#, который можно использовать для этого.
все-таки там 2.5 миллиарда клеток в секунду, а не миллиона.
А с использованием SSE-команд можно еще быстрее, до 15 миллиардов в секунду: habr.com/ru/post/505606
А в этой статье показывают, что на GPU можно достичь скорости, если я не ошибся в подсчетах, полтора триллиона клеток в секунду!
Не совсем: проверка числа на простоту, в отличие от факторизации, решается за полиномиальное время
ݻටຈݹ௫ටຯѧ௨پଈݶಏقງݹலටະऄ௨پໃݢ௨ගɑݸҰقະࡨ௨sໃɥ௨ڏІݟ௨ගɱݸతටࢩݹࢴටแɕ௧ڄໃףۑඤÐݸలقลڝ௨ටໃݹ௨ටࠋޛටדݹ௨ඞТஈ௧ƣເঅ௩ɦܘݶహइໂђ௨ටໞݹஜɦܘݚذචʅݹਐටࢺݹસقໃݹ௨ටໃݹਆИ я в процессе решения задачи для глубины 6. Пока что он нашел score=6 (поле перебора трех миллиардов позиций)
ܮටลຍ௨ಀඍݸహටາɝɷචŒݹமϺາŦ௨ಀІݟௐقГঅஶѼໃӌ௨ගɱݸƓΟເח௧ڈໃӻ௨ගȣݸಏටງђஶɈժށѻටໃݹ௨ටເ0но алгоритм еще не закончил работу, пока перебрав восемь миллиардов позиций за 11 часов.
На моем компьютере (у меня Windows/cygwin, поэтому мне сложно сравнить с вашими результатами напрямую) этот алгоритм работает 1.7 секунд. Я Вам послал пул реквест, было бы интересно узнать, будет ли он обгонять reusebuffer в ваших тестах. Правда, сейчас он работает только для LIMIT = 10^k, но, думаю, его можно доделать для любого диапазона чисел без ухудшения производительности.