Решил поискать, что нового есть по этой задаче в Интернете, и обнаружил объяснение, почему задача известна как "Partridge puzzle". Я знал, что слово "partridge" переводится как "куропатка", и предположил, что это связано с пестротой рисунка. А оказалось, что это слово из песни "Двенадцать дней Рождества" ("12 Days of Christmas") [слова], [википедия], которая строится по принципу игры "Снежный ком": сначала любимый дарит 1 подарок на Рождество (куропатку на грушевом дереве), потом 2 подарка, потом 3 подарка, и так далее до 12 подарков. Не сложно заметить схожесть с формулировкой задачи.
Во-первых, результат для n=9 подтвердился. Я взял код №2 от @lightln2, распараллелил его с помощью обвязки на Питоне и запустил на 10 ядрах. Результат сошёлся с моим.
Да, вы правы, вариантов с 1, которые продлить нельзя, будет много, но, наверное, их можно отфильтровать дополнительным условием. Может тогда и для 9 бы сработало, если бы я об этом знал.
Не совсем понял про "яму" длиной 2 и высотой 4. Если это про метод, проиллюстрированный на Рис. 4, то это решается рассмотрением вариантов квадратиков, у которых основание лежит на нужной половине, а центр уже за границами половины, и исключением их из финальной фигуры.
Вы, наверно, можете в oeis.org завести соответствующую последовательность (будет что-то типа 1,0,0,0,0,0,0,2332,216285)
Да, я как раз подал на рассмотрение в oeis.org такую последовательность.
Я написал простую реализацию на c++ (однопоточную и без учета симметрий), она решила n=8 за 25 минут.
Это очень впечатляет!
Если я правильно понял, то вы перешли от dfs (depth-first search), не очень подходящего для GPU, к bfs (breadth-first search)?
По сути, да. Делаю 1 шаг поиска в ширину, а потом каждую подзадачу пробую решить поиском в глубину.
А всех вариантов действительно будет в 8 раз больше? Это верно, только каждый вариант разбиения несимметричен ни в какой симметрии, так как иначе ему будет соответствовать меньше, чем 8 вариантов. Я проверил для n=8 и действительно получил 18656 = 2332*8 вариантов. Интересно, можно ли доказать, что не существует симметричных разбиений?
Поясню этот момент. Это утверждение мне казалось очевидным, но в реальности оно чуть более сложное.
Во-первых, поворот на 90 и 180 градусов не может привести к тому же заполнению, т.к. у нас более одного квадратика в нечётном количестве, т.е. нечётное количество квадратиков одного размера стоит вне центра, и поэтому хотя бы один квадратик не найдёт соответствия, т.к. все циклы "кто в кого переходит при повороте" чётной длины.
Во-вторых, при n=8 (и других n, когда сторона квадрата чётная) симметрии быть не может, т.к. квадратики нечётного размера не могут стоять в центре, а их нечётное количество.
В-третьих, при n=9 (и других n, когда сторона квадрата нечётная) симметрия потенциально может быть двух видов: относительно горизонтали/вертикали и относительно диагонали. В первом случае на центральной линии, а во втором на диагонали должны стоять только квадратики размера 1, 3, 5, 7, 9. Проблема возникает с квадратиком размера 1. Если симметрия относительно горизонтали/вертикали, то он будет стоять между двумя квадратиками большего размера, что приведёт к "пазу" толщины 1, или же он будет у границы. В случае диагонали, несложно понять, что квадратик размера 1 также не может стоять в таком положении, т.к. это приведёт к "пазу" толщины 1.
Я же правильно понимаю, что в некоторых случаях хеджинг может не помочь? Например, когда нагрузка небольшая, каждый запрос выполняется на одном ядре и есть существенное различие в количестве вычислений разных запросов (скажем, в 2-3 раза).
Я использовал nVidia A100 40Gb :) Мой код отрабатывает за 12 -> 23.936
Спасибо за интересные оптимизации.
Чтобы a было в нужном диапазоне, надо rev(b) < b, и можно рассматривать только те b, у которых последняя цифра не больше первой (а также не 0).
Я тоже такое пробовал, но почему-то это не давало прироста в скорости, хотя я и выравнивал интервалы по значениям, кратным 30. Наверное, нужно было как у вас явно выделить b_first.
import time
import numpy as np
import math
from numba import cuda
from numpy import uint64, int64, int32, uint8
BASE = uint64(10)
@cuda.jit
def cuda_run_core(A, oA, oB, B, d):
mod9 = cuda.local.array(9, uint64)
mod9[0] = 0
mod9[1] = 0
mod9[2] = 2
mod9[3] = 6
mod9[4] = 0
mod9[5] = 8
mod9[6] = 3
mod9[7] = 0
mod9[8] = 5
mod11 = cuda.local.array(11, uint64)
if d % 2 == 0:
mod11[0] = 0
mod11[10] = 0 # stub
mod11[9] = 9
mod11[8] = 4
mod11[7] = 6
mod11[6] = 7
mod11[5] = 1
mod11[4] = 8
mod11[3] = 2
mod11[2] = 3
mod11[1] = 5
else:
mod11[0] = 0
mod11[1] = 0 # stub
mod11[2] = 9
mod11[3] = 4
mod11[4] = 6
mod11[5] = 7
mod11[6] = 1
mod11[7] = 8
mod11[8] = 2
mod11[9] = 3
mod11[10] = 5
div2_mod9_x11 = cuda.local.array(18, uint64) # x/2 mod 9
for i in range(2):
div2_mod9_x11[i*9+0] = 0*11
div2_mod9_x11[i*9+1] = 5*11
div2_mod9_x11[i*9+2] = 1*11
div2_mod9_x11[i*9+3] = 6*11
div2_mod9_x11[i*9+4] = 2*11
div2_mod9_x11[i*9+5] = 7*11
div2_mod9_x11[i*9+6] = 3*11
div2_mod9_x11[i*9+7] = 8*11
div2_mod9_x11[i*9+8] = 4*11
shift = cuda.shared.array(99, uint8)
b = cuda.threadIdx.x
if b < 99:
a = mod11[b % 11]
a += div2_mod9_x11[9 + mod9[b % 9] - a % 9]
shift[b] = a
reverse = cuda.shared.array(100, uint64)
i = cuda.threadIdx.x
if i < 100:
reverse[i] = cuda_rev2(i)
BASE2 = uint64(pow(BASE, 2))
cuda.syncthreads()
B0 = B // BASE
pos = cuda.grid(1)
if pos + 1 < A.size:
b_lo = A[pos]
b_hi = A[pos + 1]
b = b_lo - 1
while b < b_hi - 1:
b += 1
if b % 3 == 1:
b += 1
if b % BASE == 0:
b += 1
if b >= b_hi:
break
rb = cuda_rev(b, reverse, BASE2)
a0 = uint64(math.ceil(rb * 1.0 * B / b))
a_lo = a0 - 20
a_hi = a0 + BASE + 20
a = a_lo - a_lo % 99 + shift[b % 99]
if a < a_lo:
a += 99
if a < a_hi and B0 <= a < B:
ra = cuda_rev(a, reverse, BASE2)
err = uint64(a * b - rb * B)
if err == ra:
oA[pos] = a
oB[pos] = b
@cuda.jit(device=True)
def cuda_rev2(x): # precomputation for 2-digit values
r = uint64(0)
for i in range(2):
r = r * BASE + x % BASE
x = x // BASE
return uint64(r)
@cuda.jit(device=True)
def cuda_rev(x, reverse, BASE2): # 33% speedup for 100-valued reverse
r = uint64(0)
while x > 0:
r = r * BASE2 + reverse[x % BASE2]
x = x // BASE2
while r % BASE == 0:
r = r // BASE
return uint64(r)
def rev(x): # CPU version
r = 0
while x > 0:
r = r * BASE + x % BASE
x = x // BASE
return r
def run_cuda(d):
num = 2**18
A = np.linspace(pow(BASE, d - 1), pow(BASE, d), num).astype(np.uint64)
oA, oB = np.zeros_like(A), np.zeros_like(A)
threadsperblock = 2**9
blockspergrid = (A.size - 1 + (threadsperblock - 1)) // threadsperblock
A, oA, oB = cuda.to_device(A), cuda.to_device(oA), cuda.to_device(oB)
cuda_run_core[blockspergrid, threadsperblock](A, oA, oB, uint64(pow(BASE, d)), uint64(d))
A, oA, oB = A.copy_to_host(), oA.copy_to_host(), oB.copy_to_host()
for x, y in zip(oA, oB):
if x:
print(x, '*', y, '=', rev(y), rev(x))
def main():
for d in range(6, 17):
t0 = time.perf_counter()
run_cuda(d)
t1 = time.perf_counter()
print(f'{d} -> {t1-t0:.3f}')
if __name__ == '__main__':
main()
Используются свойства остатков от деления на 3, 9, 11. Остаток отделения на 3 и 9 равен остатку от деления на 3 и 9 суммы цифр числа. Остаток от деления на 11 равен остатку от деления на 11 знакопеременной суммы цифр начиная с конца. Пусть a, b: a*b = rev(b)rev(a). Тогда для каждого b возможно не более одного остатка от деления a на 9 и не более одного остатка от деления a на 11. Это легко проверить, перебрав все пары остатков от деления и сравнив их произведения с суммой или разностью (в случае нечётных чисел и деления на 11). Если эти свойства объединить, то получаем для каждого b не более 1 остатка от деления на 99. Его можно получить из уравнения a = 99t + 9k + mod9[b] = 99t + 11l + mod11[b]. Далее, поскольку 0 <= l < 9, можем перейти к уравнению по модулю 9 и выразить l через деление на 2 по модулю 9.
Остаток от деления a и b на 3 не может быть равен 1. Однако если мы в цикле сделаем continue, то код не будет быстрее работать, поскольку в соседних нитях остаток от деления может быть не равен 1, а, если я правильно помню из курса по CUDA, все нити одного блока всегда выполняют одинаковую строчку кода. Поэтому вместо continue используется += 1 и цикл while вместо for.
Для каждого b в оригинальной версии кода рассматривался диапазон от a0 до a0+BASE. Но при n=16 у нас последняя цифра может быть неправильной, поскольку точность типа double в общем случае только 15 знаков. Поэтому я взял диапазон от a0-20 до a0+BASE+20, чтобы заведомо покрыть все нужные значения (предпоследняя цифра вроде тоже может быть округлена в +-1). В этот диапазон попадает не более 1 остатка от деления на 99.
Я также ускорил функцию rev, предрассчитав обратные для пар цифр, таким образом делается в 2 раза меньше витков цикла.
Я ускорил ваш алгоритм примерно в 4 раза и расширил его до n=16 (в умножении типа double уже не хватает точности в 15 знаков). Таким образом удалось проверить n=15 и n=16. На это ушло примерно 3,5 дня.
Для n=15 ничего не нашлось, для n=16 нашлась одна пара:
Позже повнимательнее посмотрел - на самом деле требуется не так уже много памяти для n=7 - чуть более 2 Гб. Поэтому запустил код на 11 ядрах и где-то за 12 часов воспроизвёл результат @gchebanov213002162.
Я взял оценку числа неупорядоченных пар для n=3..18, затем для n=4..18 посчитал коэффициент роста этой последовательности (относительно предыдущего значения) и построил график. Поскольку на бесконечности коэффициент роста не ожидается более 100, здесь явно не логарифмическая функция. Я предположил, что это экспоненциальная функция со смещением, и приближение получилось довольно точным. Приближал методом градиентного спуска с клипированием. Таким образом на бесконечности коэффициент роста ожидается приблизительно 81,6.
Спасибо! При d=18 коэффициент роста последовательности оценки количества пар равняется примерно 77,2. Попробую построить прогноз, к чему придём на бесконечности.
Решил поискать, что нового есть по этой задаче в Интернете, и обнаружил объяснение, почему задача известна как "Partridge puzzle". Я знал, что слово "partridge" переводится как "куропатка", и предположил, что это связано с пестротой рисунка. А оказалось, что это слово из песни "Двенадцать дней Рождества" ("12 Days of Christmas") [слова], [википедия], которая строится по принципу игры "Снежный ком": сначала любимый дарит 1 подарок на Рождество (куропатку на грушевом дереве), потом 2 подарка, потом 3 подарка, и так далее до 12 подарков. Не сложно заметить схожесть с формулировкой задачи.
Есть 2 хорошие новости:
Во-первых, результат для n=9 подтвердился. Я взял код №2 от @lightln2, распараллелил его с помощью обвязки на Питоне и запустил на 10 ядрах. Результат сошёлся с моим.
Во-вторых, последовательность приняли в OEIS: https://oeis.org/A381976
Да, вы правы, вариантов с 1, которые продлить нельзя, будет много, но, наверное, их можно отфильтровать дополнительным условием. Может тогда и для 9 бы сработало, если бы я об этом знал.
Не совсем понял про "яму" длиной 2 и высотой 4. Если это про метод, проиллюстрированный на Рис. 4, то это решается рассмотрением вариантов квадратиков, у которых основание лежит на нужной половине, а центр уже за границами половины, и исключением их из финальной фигуры.
Уточнил, n=8 занимает 8 минут на GPU.
При этом у вас, как я понимаю, работает на одном ядре. Поэтому 9,5 мин - это очень крутой результат!
15 минут было до того, как я применил тест на каноничность, поэтому, предполагаю, около 7 минут. Завтра смогу написать точнее.
Спасибо. Очень интересно! По времени уже близко к тому, как мой код работает на GPU.
Спасибо за интерес к статье!
Да, я как раз подал на рассмотрение в oeis.org такую последовательность.
Это очень впечатляет!
По сути, да. Делаю 1 шаг поиска в ширину, а потом каждую подзадачу пробую решить поиском в глубину.
Поясню этот момент. Это утверждение мне казалось очевидным, но в реальности оно чуть более сложное.
Во-первых, поворот на 90 и 180 градусов не может привести к тому же заполнению, т.к. у нас более одного квадратика в нечётном количестве, т.е. нечётное количество квадратиков одного размера стоит вне центра, и поэтому хотя бы один квадратик не найдёт соответствия, т.к. все циклы "кто в кого переходит при повороте" чётной длины.
Во-вторых, при n=8 (и других n, когда сторона квадрата чётная) симметрии быть не может, т.к. квадратики нечётного размера не могут стоять в центре, а их нечётное количество.
В-третьих, при n=9 (и других n, когда сторона квадрата нечётная) симметрия потенциально может быть двух видов: относительно горизонтали/вертикали и относительно диагонали. В первом случае на центральной линии, а во втором на диагонали должны стоять только квадратики размера 1, 3, 5, 7, 9. Проблема возникает с квадратиком размера 1. Если симметрия относительно горизонтали/вертикали, то он будет стоять между двумя квадратиками большего размера, что приведёт к "пазу" толщины 1, или же он будет у границы. В случае диагонали, несложно понять, что квадратик размера 1 также не может стоять в таком положении, т.к. это приведёт к "пазу" толщины 1.
Вот пример, когда все четвёрки отдельно
Это вы точно заметили - в 99.99% случаев всегда присутствует пара четвёрок с общей стороной. Для остальных размеров квадратиков статистика такая:
2 - 93.72%
3 - 99.31%
4 - 99.99%
5 - 99.95%
6 - 99.92%
7 - 100%
8 - 100%
9 - 100%
Спасибо за высокую оценку!
Спасибо, интересная статья!
Я же правильно понимаю, что в некоторых случаях хеджинг может не помочь? Например, когда нагрузка небольшая, каждый запрос выполняется на одном ядре и есть существенное различие в количестве вычислений разных запросов (скажем, в 2-3 раза).
Я использовал nVidia A100 40Gb :) Мой код отрабатывает за
12 -> 23.936Спасибо за интересные оптимизации.
Я тоже такое пробовал, но почему-то это не давало прироста в скорости, хотя я и выравнивал интервалы по значениям, кратным 30. Наверное, нужно было как у вас явно выделить b_first.
Вот код:
Hidden text
Используются свойства остатков от деления на 3, 9, 11. Остаток отделения на 3 и 9 равен остатку от деления на 3 и 9 суммы цифр числа. Остаток от деления на 11 равен остатку от деления на 11 знакопеременной суммы цифр начиная с конца. Пусть a, b: a*b = rev(b)rev(a). Тогда для каждого b возможно не более одного остатка от деления a на 9 и не более одного остатка от деления a на 11. Это легко проверить, перебрав все пары остатков от деления и сравнив их произведения с суммой или разностью (в случае нечётных чисел и деления на 11). Если эти свойства объединить, то получаем для каждого b не более 1 остатка от деления на 99. Его можно получить из уравнения a = 99t + 9k + mod9[b] = 99t + 11l + mod11[b]. Далее, поскольку 0 <= l < 9, можем перейти к уравнению по модулю 9 и выразить l через деление на 2 по модулю 9.
Остаток от деления a и b на 3 не может быть равен 1. Однако если мы в цикле сделаем continue, то код не будет быстрее работать, поскольку в соседних нитях остаток от деления может быть не равен 1, а, если я правильно помню из курса по CUDA, все нити одного блока всегда выполняют одинаковую строчку кода. Поэтому вместо continue используется += 1 и цикл while вместо for.
Для каждого b в оригинальной версии кода рассматривался диапазон от a0 до a0+BASE. Но при n=16 у нас последняя цифра может быть неправильной, поскольку точность типа double в общем случае только 15 знаков. Поэтому я взял диапазон от a0-20 до a0+BASE+20, чтобы заведомо покрыть все нужные значения (предпоследняя цифра вроде тоже может быть округлена в +-1). В этот диапазон попадает не более 1 остатка от деления на 99.
Я также ускорил функцию rev, предрассчитав обратные для пар цифр, таким образом делается в 2 раза меньше витков цикла.
Я ускорил ваш алгоритм примерно в 4 раза и расширил его до n=16 (в умножении типа double уже не хватает точности в 15 знаков). Таким образом удалось проверить n=15 и n=16. На это ушло примерно 3,5 дня.
Для n=15 ничего не нашлось, для n=16 нашлась одна пара:
8830950118748589 * 3376117662341892 = 29814326671167339858478110590388
Надеюсь, ничего не пропустил.
Позже повнимательнее посмотрел - на самом деле требуется не так уже много памяти для n=7 - чуть более 2 Гб. Поэтому запустил код на 11 ядрах и где-то за 12 часов воспроизвёл результат @gchebanov213002162.
Я взял оценку числа неупорядоченных пар для n=3..18, затем для n=4..18 посчитал коэффициент роста этой последовательности (относительно предыдущего значения) и построил график. Поскольку на бесконечности коэффициент роста не ожидается более 100, здесь явно не логарифмическая функция. Я предположил, что это экспоненциальная функция со смещением, и приближение получилось довольно точным. Приближал методом градиентного спуска с клипированием. Таким образом на бесконечности коэффициент роста ожидается приблизительно 81,6.
Спасибо! При d=18 коэффициент роста последовательности оценки количества пар равняется примерно 77,2. Попробую построить прогноз, к чему придём на бесконечности.
А вы не против, если я ваш код приведу как пример на oeis.org?
Я вот эту последовательность хочу зарегистрировать https://oeis.org/draft/A370723
Либо можете сами зарегистрироваться и вписать в раздел PROG - там любой человек может черновик редактировать.