Comments 52
Всё выше изложенное жёстко связано с системой счисления, поэтому не представляет настолько большого интереса, как глобальные задачи теории чисел
Всё равно можно сформулировать гипотезу и назвать её своим именем (или ником). Например, так.
Малая гипотеза Света Тьмы:
Пар чисел перемножающихся образом abcd*ef = fedcba (во множителях разное количество знаков) бесконечно много вне зависимости от того, какая позиционная система счисления выбрана
Большая гипотеза Света Тьмы:
Пар чисел перемножающихся образом abc*def = fedcba (во множителях одинаковое количество знаков) бесконечно много вне зависимости от того, какая позиционная система счисления выбрана.
Интересно, чат GPT может эти теории как то истолковать, объяснить или новые найти ?
Думаю, что не может. На math.stackexchange.com мне вообще рекомендовали к нему не обращаться)
Вам и тут так посоветуют. Эти GPT лишь подбирают слоги так, чтобы текст был более менее статистически похож на то, что они видели в обучающей выборке. Все. Никакого понимания, знаний там даже близко нет. Если спрашивать что-то очень часто встречающееся в выборке - оно вам может даже точно воспроизведет. Иначе - сгаллюцинирует красивый бред.
Оно годится только там, где важна форма, а не содержание. Стихи какие-нибудь оно отлично генерирует, потому что вообще не важно, что там конкретно. А вот математические какие-то толкования тут вообще не получаются.
Стихи оно вообще не генерирует
Просто подбирают слоги на основании обучения цепи Маркова. Те, кто видел их типичный выхлоп, понимают, каких космических масштабов разница между ними и языковыми моделями.
Ну обратим внимание на то, что вы при написании вашего комментария всего лишь последовательно генерировали нажатия на клавиши, никакого понимания там близко нет...
Цепи маркова слишком локальны. Языковые модели лишь позволяют более гибко подгонять статистические свойства для более широкого контекста.
Мне обещали слоги. А вы говорите, контекст. Как только появляется контекст, то модель перестаёт "лишь подбирать слоги". Контекст и есть понимание.
Хорошо: она подбирает по одному слогу так, чтобы итоговый текст статистически походил на обучающую выборку.
Цепь Маркова тоже это делает. И то, что пишите вы, тоже обладает этим свойством.
статистически походил на обучающую выборку.
Вот только проблема в том, что статистическими методами нельзя получить длинный осмысленный текст, который отвечает на заданный вопрос.
А тот же chatgpt разве не длинные осмысленные тексты дает? Статистически. Да, часто выдает с ошибками.... Но я вот точно знаю, что люди частенько тоже не могут выдать связный осмысленный текст. Особенно в темах, где они "плавают", не разбираются, там где не хватает "обучения". Думаю вы и сами это неоднократно наблюдали, в школе, в вузе, в быту и т.д.
Просто наш биологический механизм значительно сложнее, чем текущие языковые модели и ёмкость побольше. Я как-то читал одну статью по исследованию эмоций. Как они появляются, как развиваются, трансформируются и т.д. Оказывается, даже эмоции людей - это результат обучения. Кроме некоторых врождённых реакций. И это обучение эмоциям такой сложный путь проходит, столько слоев накручивается внутри нас в процессе взросления, что к зрелому состоянию там в мозгу образуется мощнейший клубок, в котором множество слоев, связей между ними. И все эмоции (слои обученных реакций) как-то связаны ещё и с жизненным опытом, памятью, и даже биохимией (может пульс ускориться, гормоны выплеснуться, кожа покраснеть и т.д.). Причем, это все работает с обратной связью, т.е. , например, недостаток гормона какого-то может вызвать психологическую и эмоциональную реакцию (отсюда и появилось лечение антидепрессантами), голод может вызвать раздражение и т.д. Воспоминания тем же свойством обладают, поэтому неприятное воспоминание может провоцировать физическую реакцию тела, и негативно влиять на реальное состояние тела.
Так что, видится мне, что от текущих нейросетей мы отличаемся лишь невероятной биологической сложностью и вычислительной мощью для такого компактного размера. И то что текущие нейросети проще человеческой, это наоборот хорошо. И я надеюсь, что ученные будут тщательно контролировать эту сложность (структурную, а не количественную). Не думаю, что кому-то бы понравилось, если бы его автономный автомобиль испытывал раздражение при снижении уровня заряда))) Хотя, на ум при этом приходят некоторые сюжеты из фантастики, про живые космические корабли....
Вам шашечки или ехать? Знаний нет, сознания нет, понимания нет. Какая разница? Она решает задачу? Решает. С вот такой вот точностью и вероятностью. И быстрее чем мясной человек. Не везде. Вот это ещё не может. И вот это. А это может. Вот тут и используем. Может быть, что с увеличением количества параметров она догонит и перегонит мозг. А может быть и нет.
Не обязательно быть мозгом, чтобы действовать как мозг.
Она решает задачу? Решает.
Вот в том-то и проблема, что нет! Точность и вероятность там очень маленькие. А главная проблема - даже в случае ошибки оно выдает правдоподобный результат (в смысле: похожий на нужный ответ по форме).
Поэтому нельзя пользоваться им, если вы не достаточно разбираетесь в теме, чтобы уметь отличить галлюцинацию от случайного попадания в ответ. Но в этом случае оно вам и не надо было особо.
Так я и говорю о том, что вот инструмент, его удобно применять тут и не удобно применять тут. Стиральная машина стирает сама, но загружать/выгружать шмот (пока) не у умеет.
3D принтер может напечатать крючок для полотенца, но не может напечатать вентилятор в сборе вместе с подшипниками и двигателями, потому что принтеров, которые печатают и проводником, и изолятором, пока ещё нету и делать почему‑то никто особо не торопится (хотя, имхо, когда это появится массово, это перевернёт мир похлеще появления паровых двигателей — много чего электромеханического можно будет печатать по тупо по слоям, поменяются сами принципы электромеханики). Но это не значит, что 3D принтеры — плохие, и «точность и вероятность» решения проблемы там маленькая. Есть круг того, что они могут напечатать. Который постепенно растёт.
Мне надо было на C++ написать поиск процессов через WinAPI — я пнул ChatGPT, он написал код, код работает. Я бы его сам написал, но не за минуту, а подольше, потому что надо курить доки WinAPI. Мне надо было по документации на 50 страниц ответить на вопрос как конфигурировать то‑то — он дал ответ, причем почти мгновенно. Инструмент работает? Работает.
А есть масса ситуаций, где он не работает. Это вполне нормально, таким свойством обладает любой инструмент созданный человечеством. Она ошибается, она косячит, она несёт лютую дичь и её надо перепроверять. А люди разве так не делают? Вот‑вот.
Если одни начинают продавать лопаты и орать что лопата может всё, и варить двери, и суп, и пускать спутники, и рыть ямы, и мяукать, и ей можно бриться и красить дома, и вы должны немедленно её купить, и они, при этом, не правы — это вовсе не означает, что лопата не может делать вообще ничего.
Пока одни говорят, что по формуле Циолковского можно доказать, что лопата не может запускать спутники, и добавляют, что она ни разу ещё не мяукнула, связывая это с отсутствием усов, кто‑то уже копает котлован под третий дом.
Причина, по которой этот инструмент называют «интеллектом», заключается в особенностях мира, в котором мы живём: сделал что‑то → надо это продать. А слово «интеллект» очень хорошо помогает донести до большинства полезность инструмента. Да, это преувеличение. Но нет, тот факт, что это преувеличение, не является доказательством того, что этот инструмент не решает задачи.
Он решает определённый круг задач. Который постепенно растёт.
Интересны детали алгоритма. Почему для n > 6 требуется много памяти.
Если вкратце, в произведении i*j=m перебираются всевозможные i и все возможные первые t=4 цифры числа m (обозначим за x). По этим значениям определяется, каким может быть число j – оно лежит в некотором диапазоне и при этом должно обязательно содержать цифры, которые входят в x и не входят в i. Списки чисел, которые содержат определённый набор цифр предрассчитываются, поэтому требуется много памяти при больших n.
Позже повнимательнее посмотрел - на самом деле требуется не так уже много памяти для n=7 - чуть более 2 Гб. Поэтому запустил код на 11 ядрах и где-то за 12 часов воспроизвёл результат @gchebanov213002162.
Кстати, довольно просто показать, что решение первоначальной задачи существует для любого n >=2. Если у нас есть решение a * b с n цифрами, то 10a * 10b будет решением для случая n+1 цифр. И отсюда же следует, что количество решений с ростом n будет, как минимум, не убывать.
Посмотрел в базе числовых последовательностей и совпадений не нашёл. Советую обратиться к математикам-блоггерам с определением формулы. Кстати, вот знатные любители таких задачек. Ну и в целом публикуй для потомков в базе. Круто, но странно что это не придумали раньше.
Нашёл последовательность для переворачивающихся чисел не обязательно одинаковой длины: A009944 - OEIS
182738252812 туда тоже входит, но его официально не подтвердили, поэтому по нему не ищется.
вычисления длились 5 часов на 11-ядрах процессора
Чет медленно. Сейчас сделаю побыстрее.
@cuda.jit
def cuda_run_core(A, oA, oB, B):
B0 = B // BASE
pos = cuda.grid(1)
if pos + 1 < A.size:
b_lo = A[pos]
b_hi = A[pos + 1]
for b in range(b_lo, b_hi):
if b % BASE == 0:
continue
rb = rev(b)
a0 = uint64(math.ceil(rb * 1.0 * B / b))
for a in range(a0, a0 + BASE):
if not (B0 <= a < B):
continue
ra = rev(a)
err = uint64(a * b - rb * B)
if err == ra:
oA[pos] = a
oB[pos] = b
if err > B:
breakПроверяет все до m=12 за минуту на видюшке. Посмотрим что будет с m=14 через пару часов.
Hidden text
@cuda.jit(device=True)
def cuda_rev(x):
r = uint64(0)
while x > 0:
r = r * BASE + x % BASE
x = x // BASE
return uint64(r)
def run_cuda(d):
# 10 -> 0.638
# 11 -> 6.601
# 12 -> 68.787
# 13 -> 775.907
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)))
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(8, 16):
t0 = time.perf_counter()
run_cuda(d)
t1 = time.perf_counter()
print(f'{d} -> {t1-t0:.3f}')Ничего не нашлось, а пятнадцать нужно 16 часов ждать. Расширил пространство на несколько порядков, вдруг пригодится потомкам.
Другие основания, из не указанного выше
M=2
(40) 1100111001101110111101000000000100000011 * 1010110100001000010111011110000111010001 = 1000101110000111101110100001000010110101 1100000010000000001011110111011001110011
(42) 111010101010110001000001000101111001101011 * 101100101100110011100110011110011111000101 = 101000111110011110011001110011001101001101 110101100111101000100000100011010101010111
M=3
(18-26) None
M=4
(11) 20221122201 * 33313301202 = 20210331333 10222112202
(12) 220032112122 * 220231112121 = 121211132022 221211230022
(13) 3033122211023 * 2021000003321 = 1233000001202 3201122213303
(19) 3123010113323033123 * 1223000330230123211 = 1123210320330003221 3213303233110103213
(20) 21123300031031323121 * 33203011213013222012 = 21022231031211030233 12132313013000332112
(21) 231122200023011022212 * 200211130021013310311 = 113013310120031112002 212220110320002221132
(21) 131012330133313121331 * 220013120002032330201 = 102033230200021310022 133121313331033210131
M=5
(10) 4402131124 * 1403420431 = 1340243041 4211312044
(10) 4424311302 * 2243042322 = 2232403422 2031134244
(13) 1322024130211 * 4242031404321 = 1234041302424 1120314202231
(15) 204304333140412 * 314110410304121 = 121403014011413 214041333403402
(16-18) None
M=6
(10) 2344420112 * 2433414111 = 1114143342 2110244432
(10) 2230332102 * 3120354411 = 1144530213 2012330322
(10) 1440535111 * 5043054031 = 1304503405 1115350441
(14) 21240011132212 * 42005450351431 = 13415305450024 21223111004212
M=7
(13) 5054030524163 * 5540446553414 = 4143556440455 3614250304505
(14) 26623462035461 * 52046210055512 = 21555001264025 16453026432662
M=8
(10) 6774770256 * 1540253631 = 1363520451 6520774776
(10) 6003775002 * 4000002003 = 3002000004 2005773006
(10) 5470206021 * 7431456225 = 5226541347 1206020745
(13) 6000377750002 * 4000000020003 = 3000200000004 2000577730006
M=9
(7-13) None
M=10
(11-14) None
M=11
(12) 534338726a98 * 598050086292 = 292680050895 89a627833435
M=12
(11) a0a40435495 * 3067a602862 = 268206a7603 59453404a0a
(11) 96124b45043 * 4a8b05766a3 = 3a66750b8a4 34054b42169
M=13
(11) 628ba56c626 * 22a0bc87901 = 10978cb0a22 626c65ab826
M=14
(8-11) None
M=15
(10) 59bcc30787 * e5a4923465 = 5643294a5e 78703ccb95
M=16
(3) 269 * e12 = 21e 962
(6) ab5af2 * 848b85 = 58b848 2fa5ba
(9) 6ebade3a2 * 82b559883 = 388955b28 2a3edabe6
M=17
(2) 6a * c4 = 4c a6
(5) 473d7 * ceg53 = 35gec 7d374
(8) 88b918f2 * 8g554484 = 484455g8 2f819b88
(9) gad3cc3g9 * f8701ac2f = f2ca1078f 9g3cc3dag
(10) None
M=18
(2) 4b * a2 = 2a b4
(5) c060c * 19001 = 10091 c060c
(9) ahb305ahe * 35c8a6702 = 2076a8c53 eha503bha
(9) ad95ea864 * c531fe167 = 761ef135c 468ae59da
(9) a3h94a05h * eh2g90f88 = 88f09g2he h50a49h3a
M=19
(2) 62 * b3 = 3b 26
(3-10) None
M=20
(3) ac5 * j2a = a2j 5ca
(4) ce24 * 4j23 = 32j4 42ec
(5) g040g * 15001 = 10051 g040g
(8) 44ch32i2 * ac2jgh42 = 24hgj2ca 2i23hc44
(9) 4f55610hc * 8bb1h2h02 = 20h2h1bb8 ch01655f4
M=21
(3) 7g7 * ce4 = 4ec 7g7
(7) h41830j * 306fh92 = 29hf603 j03814h
M=22
(6) g5ldhf * 5d5334 = 4335d5 fhdl5g
M=23
(2-8) None
M=24
(5) i060i * 18001 = 10081 i060i
(5) g080g * 1c001 = 100c1 g080g
(5) m823h * ff4de = ed4ff h328m
(5) l7ch3 * h305f = f503h 3hc7l
(5) 6in21 * m4h66 = 66h4m 12ni6
(7) 68kcfj3 * 9jbkbe2 = 2ebkbj9 3jfck86
(8) 62c9lnci * c6f1gm23 = 32mg1f6c icnl9c26
(9) NoneВаше умение обращаться с GPU впечатляет!
Я ускорил ваш алгоритм примерно в 4 раза и расширил его до n=16 (в умножении типа double уже не хватает точности в 15 знаков). Таким образом удалось проверить n=15 и n=16. На это ушло примерно 3,5 дня.
Для n=15 ничего не нашлось, для n=16 нашлась одна пара:
8830950118748589 * 3376117662341892 = 29814326671167339858478110590388
Надеюсь, ничего не пропустил.
А можете рассказать, как Вы ускорили? можно код посмотреть? Я пытался написать код на c++ CUDA, но у меня медленнее получилось (но я не double использовал, а 128-битные числа)
Вот код:
Hidden text
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 раза меньше витков цикла.
Я запустил ваш код на моей 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. Вроде бы он не так часто выполняется, и пока потоки ждут ответа от памяти, другие потоки могут выполнять вычисления. Но я не профилировал, поэтому не уверен.
Если интересно, вот код, для простоты поддержал только четные 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;
}
Я использовал nVidia A100 40Gb :) Мой код отрабатывает за 12 -> 23.936
Спасибо за интересные оптимизации.
Чтобы
aбыло в нужном диапазоне, надоrev(b) < b, и можно рассматривать только теb, у которых последняя цифра не больше первой (а также не 0).
Я тоже такое пробовал, но почему-то это не давало прироста в скорости, хотя я и выравнивал интервалы по значениям, кратным 30. Наверное, нужно было как у вас явно выделить b_first.
Переписав на С++, на вычисление ушло 25 минут
Лобовой подсчет без каких-либо алгоритмов на видяшке работает 4 секунды! (задача с парами). Если не лень ждать 16 часов сможете узнать ответ для n=7.
7
2 -> 0.675
156
3 -> 0.001
3399
4 -> 0.044
112025
5 -> 4.365
4505706
6 -> 537.973Hidden text
@cuda.jit(device=True)
def cuda_pair_check(i, j, k, d, base):
a = cuda.local.array(base, np.uint8)
for z in range(base):
a[z] = 0
for z in range(d):
a[int(i % base)] += 1
i = i // base
a[int(j % base)] += 1
j = j // base
for z in range(2*d):
a[int(k % base)] -= 1
k = k // base
nz = 0
for z in range(base):
if a[z] == 0:
nz += 1
return nz == base
@cuda.jit
def cuda_pairs_core(A, B, d, base):
assert base == 10
p_i, p_j = cuda.grid(2)
t_i, t_j = cuda.threadIdx.x, cuda.threadIdx.y
low_lim = 1
for k in range(2*d - 1):
low_lim *= base
if p_i + 1 < A.size and p_j + 1 < A.size:
i0, i1 = A[p_i], A[p_i + 1]
j0, j1 = A[p_j], A[p_j + 1]
cnt = uint64(0)
for i in range(i0, i1):
i = uint64(i)
for j in range(j0, j1):
j = uint64(j)
if i <= j:
k = uint64(i * j)
if k < low_lim: continue
if cuda_pair_check(i, j, k, d, 10):
cnt += 1
cuda.atomic.add(B, (t_i, t_j), cnt)
def run_cuda_pairs(d):
num = 2**10
A = np.linspace(pow(BASE, d - 1), pow(BASE, d), num).astype(np.uint64)
threads = (2**5, 2**4)
B = np.zeros(threads, dtype=np.int32)
blocks = (
(A.size - 1 + (threads[0] - 1)) // threads[0],
(A.size - 1 + (threads[1] - 1)) // threads[1],
)
A, B = cuda.to_device(A), cuda.to_device(B)
cuda_pairs_core[blocks, threads](A, B, d, BASE)
B = B.copy_to_host()
return np.sum(B)
def main_pairs():
print(f'start {datetime.datetime.now()}')
for d in range(2, 8):
t0 = time.perf_counter()
cnt = run_cuda_pairs(d)
t1 = time.perf_counter()
print(cnt)
print(f'{d} -> {t1-t0:.3f}')
Можно код для Python?
Исходная задача
Hidden text
for n in range(1, 4+1):
p = 10**(n-1)
count = 0
for i in range(p, p*10):
for j in range(i, p*10):
if sorted(str(i)+str(j)) == sorted(str(i*j)):
count += 1
print(n, count)Переворачивающиеся при умножении числа
Hidden text
MIN = 1
MAX = 8
divisors_mod10 = [[[] for j in range(10)] for i in range(10)]
for i in range(10):
for j in range(10):
divisors_mod10[(i*j)%10][i].append(j)
def f(i, si, n, k, sj, accum):
if k == n:
if sj[0] >= sj[n-1]:
j = int("".join(map(str, sj)))
sm = "".join(map(str, (si + sj)[::-1]))
if int(sm) == i*j:
print(f"{i} * {j} = {sm}")
return
for k1 in range(k):
accum += si[n-1-k+k1]*sj[n-1-k1]
vd = divisors_mod10[(si[k] - (accum % 10) + 10) % 10][si[n-1]]
for divisor in vd:
sj[n-1-k] = divisor
if k == 0 and si[0] < sj[n-1]:
continue
f(i, si, n, k+1, sj, (accum + si[n-1]*sj[n-1-k])//10)
for n in range(MIN, MAX+1):
print(f"n={n}")
p = 10**(n-1)
for i in range(p+1, p*10):
if i % 3 == 1:
continue
if i % p == 0:
print(i)
si = list(map(int, str(i)))
sj = [0]*len(si)
f(i, si, n, 0, sj, 0);
Есть довольно простое решение, основанное на том, что если a*b=rev(a)rev(b), то b = rev(a)*10^n/a + rev(b)/a. Второе слагаемое от 0 до 9. А первое можно округлить делением нацело. Т.е. можно для каждого a перебрать до 10 значений b:
N = 8
low = pow(10, N-1)
high = pow(10, N)
for a in range(low, high):
ra = int(str(a)[::-1])
b = ra*high // a
for t in range(0, 9):
rb = int(str(b)[::-1]);
if a*b == high*ra+rb:
print (a, b)
b += 1Работает примерно так же как и ваше решение (для N=7 ваше отработало за 20 секунд, мое за 24). Но на порядки проще. Можно чуть чуть ускорить, введя более логичные границы на сдвиг t (смотря на последнюю цифру a, например). Если там брать 9//a[0] а не 9, то работает уже за 8 секунд. Плюс, перевороты через строки и обратно - это очень тупо. Но я не питонист.
Интересная идея!
Понимаю, что вы не питонист, но по-моему, range должен идти от 0 до 10+1, поскольку последнее значение не включается и теоретически может так оказаться, что дробная часть первого и второго слагаемого в сумме дают ещё один. В усовершенствованном варианте получается range(0, 10//a[0] + 1).
Ну и если добавить проверки на a % 3 != 1 и b % 3 != 1, то получится ещё почти в 2 раза ускорение.
А, да, про это не заметил. Но до N=8 это на ответ не повлияло нигде.
И еще, могли бы вы объяснить логику этого вашего решения? Что считает функция f, что храниться в si, sj и т.д?
Я рассматривал равенство i*j = rev(j)rev(i) и перебирал всевозможные варианты числа i, то есть в отличие от вас, я перебирал младшие разряды произведения. Допустим для n = 4: abcd*xyzt = tzyxdcba, где xyzt нам неизвестно. Сначала рассматриваем разряд единиц: a = (d*t)%10. Чтобы найти t, я предрассчитал все варианты t для всех пар a и d. В цикле перебираем все варианты t и запускаем рекурсию, переходя к следующему разряду.
Разряд десятков: b = (c*t + d*z + (d*t)//10)%10. Нам уже известно значение t, поэтому можем в модульной арифметике записать (d*z) % 10 = (b - c*t + (d*t)//10) % 10. И снова находим z пользуясь предрассчитанными значениями, в цикле перебираем все варианты и запускаем рекурсию, переходя к следующему разряду.
Разряд сотен: c = (b*t + c*z + d*y + (c*t + d*z + (d*t)//10)//10)%10. Снова находим d*y в модульной арифметике и так далее.
Значение в скобках я обозначил accum - его не надо пересчитывать полностью, а только добавлять новые слагаемые и делить на 10.
В некоторых случаях оказывается, что вариантов очередной цифры нет и рекурсия просто завершается (например (2*t)%10 = 3 не имеет решений). Если дошли до конца, то проверяем, выполняется ли требуемое равенство для найденного значения xyzt.
P.S. si и sj хранят цифры чисел i и j.
А вы не против, если я ваш код приведу как пример на oeis.org?
Я вот эту последовательность хочу зарегистрировать https://oeis.org/draft/A370723
Либо можете сами зарегистрироваться и вписать в раздел PROG - там любой человек может черновик редактировать.
Спасибо, интересные задачи!
я придумал усовершенствованный алгоритм, который позволил получить значение для 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 минут.
Попробовал ваш алгоритм на видяшке, работает 4.5 секунды для n=6.
Заодно и осилил n=7
7
2 -> 0.871
156
3 -> 0.002
3399
4 -> 0.002
112025
5 -> 0.018
4505706
6 -> 4.435
213002162
7 -> 974.439Hidden text
@cuda.jit
def cuda_pairs_core(A, B, C, C1, C2, d, base):
# ...
if p_i + 1 < A.size and p_j + 1 < A.size:
i0, i1 = A[p_i], A[p_i + 1]
j0, j1 = A[p_j], A[p_j + 1]
cnt = uint64(0)
for i in range(i0, i1):
k = uint64(i * max(j0, i))
k_hi = uint64(k // b_lim)
k_lo = uint64(k % b_lim)
for j in range(max(j0, i), j1):
if k >= low_lim:
if C[k_hi] + C1[k_lo // c2_lim] + C2[k_lo % c2_lim] == C[i] + C[j]:
cnt += uint64(1)
k += i
k_lo += i
if k_lo >= b_lim:
k_lo -= b_lim
k_hi += uint64(1)
cuda.atomic.add(B, (t_i, t_j), cnt)
@njit
def init_c(C, d, BASE):
p = 1
for i in range(d):
for j in range(p - 1, -1, -1):
for z in range(BASE - 1, -1, -1):
C[j * BASE + z] = C[j] + (uint64(1) << (z * 4))
p *= BASEВсе 4 оптимизации очень достойные!
У меня интересная ошибка хабра происходит. Если запостить любой комментарий, то весь код на странице сжимается в одну строку. У кого-нибудь еще это просиходит?
Было бы любопытно узнать, как график ведёт себя на бесконечности
Прогнал Монте-Карло на равномерных неупорядоченных парах, (это не точно учтёт пары где x=y, но так как их немного, погрешность должна быть небольшой).
Просто взял p:=cnt_good/cnt_iters; estimate:=81*10^{d-2}*p/2;
Скорость симуляции получилась около 5 миллиардов итераций в секунду, что не впечатляет. Пришлось немного попотеть что бы правильно работать с 36-значными числами.
Monte-Carlo
# [x y]
# d=d elapsed_time p_good estimate_pairs
# [269 581]
# d=3 129.385 3.851820e-04 1.559987e+02
# [4395 8691]
# d=4 154.701 8.392064e-05 3.398786e+03
# [87365 79139]
# d=5 155.051 2.766043e-05 1.120248e+05
# [775100 162389]
# d=6 154.746 1.113085e-05 4.507995e+06
# [2530020 4761510]
# d=7 207.304 5.262449e-06 2.131292e+08
# [85693413 36003975]
# d=8 233.018 2.778503e-06 1.125294e+10
# [345731723 932322371]
# d=9 234.242 1.595132e-06 6.460286e+11
# [8388016499 5721446594]
# d=10 237.009 9.736462e-07 3.943267e+13
# [40264350366 98521508913]
# d=11 237.098 6.248657e-07 2.530706e+15
# [919177426937 734068423163]
# d=12 237.305 4.179230e-07 1.692588e+17
# [8616586083842 7797719745830]
# d=13 267.747 2.884882e-07 1.168377e+19
# [45126923716086 53984469270216]
# d=14 270.713 2.058496e-07 8.336909e+20
# [204357829294211 514926408573260]
# d=15 269.294 1.502094e-07 6.083479e+22
# [6896406641595237 8139155118610023]
# d=16 319.639 1.116254e-07 4.520827e+24
# [62459283066899325 52619943312077259]
# d=17 343.897 8.466919e-08 3.429102e+26
# [409783614423487398 781679240519529663]
# d=18 344.870 6.535767e-08 2.646985e+28Относительная точность падает при уменьшении p, поэтому оценку для d=18 получилось посчитать с точностью ±0.3%
Еще я сомнительным способом генерировал случайное целое, uniform_uint64() % (HI-LO) + LO В районе 10^18 там уже значительный перекос в пользу меньших чисел.
Спасибо! При d=18 коэффициент роста последовательности оценки количества пар равняется примерно 77,2. Попробую построить прогноз, к чему придём на бесконечности.
Я взял оценку числа неупорядоченных пар для n=3..18, затем для n=4..18 посчитал коэффициент роста этой последовательности (относительно предыдущего значения) и построил график. Поскольку на бесконечности коэффициент роста не ожидается более 100, здесь явно не логарифмическая функция. Я предположил, что это экспоненциальная функция со смещением, и приближение получилось довольно точным. Приближал методом градиентного спуска с клипированием. Таким образом на бесконечности коэффициент роста ожидается приблизительно 81,6.
Переворачивающиеся при умножении числа