не уверен насчет оптимизатора js, но после оптимизатора с++ не будет ни того, ни другого: деление превратится в умножение на reciprocal (throughput = 1 инструкция/такт), а if - в conditional move (2 инструкции/такт). Мне кажется, ни то, ни то не будет узким местом, скорее уж обращение к памяти

я пытался строить все варианты с половиной кубиков (18 для n=8, 22 или 23 для n=9).
В вашем случае границы, как на рис.4, сложность в том, что на любой позиции над средней линией можно поставить квадратик размера 1 (допустим, он бы еще не был использован), и тоже получится правильная граница. То есть, "жадное" заполнение сверху вниз и слева направо не работает, надо перебирать больше вариантов, по крайней мере, близко к границе. Или я чего-то не учитываю?

А сколько у вас времени работало для n=8?

Я немного поэкспериментировал с задачей, но комментарии хабра слишком коротки, чтобы подробно описать результат, поэтому вот очень примерное описание решения задачи для n=8 (n=9 я пока даже не трогал):

1. Неоптимимизированное решение: стандартный DFS (depth first search with backtracking), состояние представленно как три 16-байтных массива (текущие высоты, длины, и количество оставшихся квадратов). После заполнения текущего нулевого уровня, все заполненные строки схловываются, как в тетрисе, так что заполняется всегда нулевой уровень, что упрощает код.

Время работы: 24 мин, количество шагов DFS: 64G, код

2. Симметрия: рассматриваются только такие состояния, в которых самый первый квадрат не меньше второго и не меньше последнего в первом ряду, остальные состояния достраиваются с помощью симметрии.

Время: 14 мин, шагов DFS: 35G, код

3. Эвристика: после добавления квадрата, рассматриваются два предыдущих: если средний квадрат в этой тройке ниже, чем крайние, то проверяется, что у нас еще остались достаточно маленьуие квадраты, чтобы в будущем заполнить этот промежуток.

Время: 12 мин, шагов DFS: 23G, код

4. Низкоуровневые оптимизации перехода состояний иcпользуют 128-битные __uint128_t и регистры SSE, а также ускоренное нахождение размеров квадратов, которые можно использовать на текущем шаге, через битовые маски

Время: 9.5 мин, шагов DFS: 23G, код

Еще пробовал, но не получилось:
5. Еще более низкоуровневые оптимизации: все операции в 128-бит или SSE, уменьшение размера состояния до 32 байт, уменьшение дорогостоящих операций обновления состояния: ускорение получается меньше 10%, а код очень сильно усложняется.
6. Я пробовал делать DFS в половину глубины, а потом сшивать попарно. К сожалению, идея не работает: обратный DFS не может сгененировать некоторые варианты (например, когда в первой половине есть "яма" длиной 2 и высотой 4)

Спасибо, интересная задача, крутое исследование! Вы, наверно, можете в oeis.org завести соответствующую последовательность (будет что-то типа 1,0,0,0,0,0,0,2332,216285)
Я тоже хочу попробовать ускорить ее на CPU, как будет больше времени)) Тут возможны низкоуровневые оптимизации для представления текущего заполнения, интересно, насколько они будут полезны.
Еще можно попробовать такую оптимизацию: при добавлении нового квадрата, если он получается выше предыдущего, проверять, что существуют свободные квадраты, которые можно поставить на предыдущий квадрат (если пред-предыдущий квадрат выше) - но не факт, что это сильно сократит количество вариантов для перебора.

Я написал простую реализацию на c++ (однопоточную и без учета симметрий), она решила n=8 за 25 минут.

то добавляем к задаче ещё 1 квадратик (всевозможные его варианты) и тем самым подразбиваем задачу примерно на 9 подзадач

Если я правильно понял, то вы перешли от dfs (depth-first search), не очень подходящего для GPU, к bfs (breadth-first search)?

216'285 вариантов заполнения ..., или в 8 раз больше, если включить в ответ все повороты и симметрии

А всех вариантов действительно будет в 8 раз больше? Это верно, только каждый вариант разбиения несимметричен ни в какой симметрии, так как иначе ему будет соответствовать меньше, чем 8 вариантов. Я проверил для n=8 и действительно получил 18656 = 2332*8 вариантов. Интересно, можно ли доказать, что не существует симметричных разбиений?

Продолжаю свой крестовый поход по NP-полным задачам

Есть хорошая книга про NP-полные задачи и экспоненциальные алгоритмы: Fedor V. Fomin and Dieter Kratsch, Exact Exponential Algorithms
там есть все, от классических решений и стандартных оптимизаций до последних достижений науки (на 2010 год)

Кажется, что это не верно для 2S, а ваша оценка точная. Например, для набора (6, 22, 33, 66): S=66, 2S=132, sum(S-a)+1 = 138. Но 106 + 222 + 33 = 137 больше 2S, но его нельзя представить с номиналом 66.

Можно получить явную формулу для решения.
Как уже отмечали в комментариях, есть четыре класса A={0}, B={4,6}, C={1,3,7,9}, D={2,8}, связанных рекуррентным соотношением

Если взять квадрат преобразования, то соответствующий граф разбивается на две части:

Это, кстати, следует из того, что граф переходов получается двудольный, то есть, рёбра есть только между двумя множествами вершин {0,1,3,7,9} и {2,4,6,8}
Получившиеся матрицы 2x2 легко диагонализуются (с собственными числами ), и можно получить формулу, аналогичную Фибоначчи, для каждого класса. Например, количество путей из цифры 1:

Так как первое слагаемое меньше единицы, то можно формулу представить в таком виде:

Естесственно, как и для Фибоначчи, эти формулы быстро теряют точность даже для небольших n, так что быстрое возведение соответствующих матриц в степень n по-прежнему самый быстрый способ.

Некий китаец посчитал количество решений до 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

#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) - и их тоже обновлять в цикле. Так как диапазон гораздо меньше, он поместится в кеш процессора.

У меня этот вариант работает за 5 минут.

А я вот хочу автору сказать спасибо за статью, и позволю не согласиться, что статья ни о чем.
Что я из нее понял:

1. Автор рассматривает обобщение гипотезы Коллатца на произвольный нечетный множитель и рассмаривает, наверное, самую простую задачу (не считая тривиальной задачи нахождения неподвижных точек) нахождения циклов длины 2. Я быстрым гуглопоиском не нашел, чтобы кто-то эту задачу до него сформулировал. И даже эта задача, по-видимому, не решается аналитически. Ну и наличие решений для и само по себе неожиданно.

2. Автор применил интересный трюк, показав, что . Тем самым он ускорил алгоритм по нахождению циклов с началом, меньшим , с у алгоритма "в лоб", до . Это позволило на питоне проверить все числа до , что позволяет высказать гипотезу, что циклов длины 2, скорее всего, больше нет.

3. Еще мне было интересно ускорить сам алгоритм, избавившись от извлечения корня и умножений с делениями (мне не удалось избавиться только от одного ). Мой вариант на с++ с первой попавшейся из интернета имплементацией bigint работает примерно в 10 раз быстрее версии автора на питоне.

4. Как к задаче подступиться аналитически, совершенно непонятно. Кажется, что (фигурные скобки - это дробная часть), и вероятность найти стремится к нулю, если в двоичном представлении корня из двух нет очень больших последовательностей из нулей, и . Это предположение верно, если корень из двух - это нормальное число в двоичной системе счисления, но это не доказано.

В общем, если то, что сделал автор, аккуратно сформулировать и красиво оформить, может, могла бы получиться интересная научная статья.

А Вы не сравнивали с их реализацией? они дают на нее в статье ссылку.
Правда, их статья десятилетней давности, а код последний раз обновлялся 5 лет назад. Ну и всегда можно напрямую автору написать, у Линдена есть профиль на линкедине

Про SIMD можно найти прямо вашу задачу на StackOverflow. Хотя, кажется, можно проще, чем там предлагают, если использовать __mm_shuffle_epi8.
Но не уверен, что это сильно поможет - возможно, вы уже приблизились к границе пропускной способности памяти