Пентадекатлон мечты

[Huz]

Yadryara в сообщении #1556552 писал(а):

Эта тема как раз и посвящена нахождению T(6,15). Просто найдя три таких числа, все кроме меня поиск прекратили. И начали обсуждать, искать и успешно находить другие цепочки чисел.

Proving that T(6,15) exists by finding an example chain was a great achievement, and proved that A119479(12) = 15. But T(6,15) itself is defined as the minimal value, so maybe we should use U(n,k) to refer to the best known upper bound, to avoid confusion.

I find it interesting to extend A292580 as far as I can. Maybe T(6,15) will remain out of reach, but there is still much room for improvement in my code (including rewriting all the Perl code to C); and maybe tomorrow someone will find a new U(6,15) that reduces the difficulty by a factor of $10^{10}$.

-- 06.06.2022, 14:20 --

EUgeneUS в сообщении #1556607 писал(а):

Нашлась 9-ка на 84 делителя (M84n9, или верхняя граница для

T(42,9)

в нотации а-файла)

Позиция $n+10$ была проверена отдельно, там 84 делителя. Но прошу перепроверить.

Confirmed, congratulations. :)

[EUgeneUS]

Huz в сообщении #1556612 писал(а):

Confirmed, congratulations. :)

Thanks! It was quickly. I checked it about 15-20 min

[VAL]

T(36,12) <= 71961074668242582598821715699989485453741795326938184455293694715
Сейчас поищу оценки сверху еще и для цепочек длины, 9-11.
Полагаю, можно будет найти и 13. Далее без ускорителей тяжело.

[Huz]

VAL в сообщении #1556628 писал(а):

T(36,12) <= 71961074668242582598821715699989485453741795326938184455293694715

Congratulations. :)

[VAL]

Huz в сообщении #1556634 писал(а):

Congratulations. :)

Thanks!

T(36,9) <= 2696228936872882173807209448676125713560660220
T(36,10) <= 29402851804008858033637044945858759579501933787851516

[Huz]

VAL в сообщении #1556643 писал(а):

T(36,8) <= 2696228936872882173807209448676125713560660220

I'll take that as a typo for T(36,9), thanks.

[Dmitriy40]

Huz
Hugo, похоже я так и не ответил Вам об устройстве ускорителей и возможности их применения для других паттернов (patterns). Отвечу здесь, возможно это будет интересно и другим.

Предварительные замечания.
Задача в нахождении цепочки чисел с постоянным количеством делителей. Поделим числа в цепочке на два класса: проверяемые (checked) и непроверяемые (nonchecked), в первых осталось найти только одно большое простое в первой степени, во вторых простых может быть два и более. Ускорять внешней программой будем лишь поиск чисел первого класса, остальную часть работы оставим за PARI/GP, который заодно и перепроверит и найденные числа первого класса.
В соответствии с КТО первое число в цепочке не может быть произвольным, а должно соответствовать форме $n=n_0+km$, где $n_0$ это некое начальное значение в зависимости от конкретного расположения малых простых в паттерне цепочки, а $m$ это модуль/шаг. Таким образом достаточно перебирать числа $n=n_0+km$ для всех $k=1,2,3...$ и проверять количество делителей чисел $n+0,n+1,n+2,...$.
Для проверяемых (checked) чисел количество делителей однозначно определяется простотой большого числа, равного $x_i=(n+i)/v_i$, где $v_i$ это коэффициент с ним, а $i$ это номер проверяемого числа в цепочке (с нуля). Таким образом поиск кандидата в цепочку свёлся к поиску некоторого количества больших простых чисел $x_i=(n_0+km+i)/v_i$ в строгом соотношении друг с другом. Найденный кандидат в цепочку потом будет допроверен в PARI/GP и только там станет (или не станет) найденной цепочкой.

Что делают мои ускорители.
В PARI/GP, из-за того что это интерпретатор, все операции довольно медленны, потому основной идеей ускорителей было замена внутреннего цикла перебора $isprime((n_0+km+i)/v_i)$ на внешнюю программу, которая возвращает вероятно подходящие $k$, относительно быстро допроверяемые потом в PARI/GP.
Для этого достаточно передать внешней программе величины $n_0$ и $m$, а она найдёт вероятно подходящие $k$. Для деления работы на кусочки можно передавать и диапазон $[k_1;k_2]$ в котором искать вероятно подходящие $k$. В текущей реализации ускорителей величины $n_0$ и $m$ жёстко зашиты в исходный код программы и передаются только $[k_1;k_2]$.
Так как все $x_i$ являются гарантированно целыми (это гарантировала КТО), то при переходе от $k$ к $k+1$ величины $x_i$ изменятся всегда на одно и то же целое значение (каждая на своё). Потому вместо вычисления $(n_0+km+i)/v_i$ на каждом шаге и для каждого проверяемого числа можно инициализировать вектор проверяемых чисел начальными значениями $x_0_i=(n_0+k_1 m+i)/v_i$ и потом лишь прибавлять к нему величину $m/v_i$. Это эффективно реализуется в SSE/AVX одной командой векторного сложения.
Далее, для проверки чисел $x_i$ на простоту используется крайне простой алгоритм проверки делимости на малые простые. Но вместо делений $x_i/p_j$ (где $p_j$ это малые простые до некоего порога) используется арифметика по модулую $p_j$, т.е. вместо $x_i$ хранятся $x_{i,j}\pmod{p_j}$ и при переходе от $k$ к $k+1$ они все пересчитываются по формуле $x_{i+1,j}=x_{i,j}+m/v_{i,j}\pmod{p_j}$. Как только ни одно из чисел $x_{i,j}$ не стало равным $0$ (равенство означает что число делится на соответствующее $p_j$) — значит нашли кандидата в цепочку. Операция суммирования по модулю реализуется тремя векторными командами AVX (в SSE шестью).
Итак, поиск кандидата в цепочку свёлся к правильной инициализации нескольких таблиц и потом циклу суммирования чисел по модулю малых простых с проверкой обнуления каждого из чисел. Так как простые малы, то достаточно отвести под каждое одно слово (word, 16 битов), что позволит проверять делимость кандидатов в цепочки на простые до $2^{15}=32768$ (почему не $2^{16}$ поясню ниже). До этого порога 3512 простых, или меньше 220 слов AVX по 256 битов каждое. Но это только на одно проверяемое число в цепочке, для 11-ти проверяемых чисел слов AVX будет уже $220\times11=2420$ или всего 75КБ.
Основной цикл перехода от $k$ к $k+1$ (порты выполнения и времена указаны для архитектуры Haswell):

Используется синтаксис ASM

;                                                               Ports   Time
.add0:          mov             EDI,0                           ;0156   1/0.25
.add:           vmovdqu         ymm15,[save_ymm+EDI]            ;23     3/0.5
                vpaddw          ymm15,ymm15,[steps+EDI]         ;15     1/0.5
                vpsubw          ymm14,ymm15,[primes+EDI]        ;15     1/0.5
                vpminuw         ymm15,ymm14,ymm15               ;15     1/0.5   текущее значение остатков по модулю простых
                vmovdqu         [save_ymm+EDI],ymm15            ;4      1/1
                add             EDI,32                          ;0156   1/0.25
                cmp             EDI,stop1                       ;0156   1/0.25
                jc              .add                            ;06     1/0.5

В save_ymm сидят текущие значения $x_{i,j}$, без разделения к какому из проверяемых чисел и какому $p_j$ оно относится, это деление выполняется внешней программой подготовки таблиц, для ускорителя это деление по индексам совершенно неважно, он проверяет факт обнуления любого из этих чисел.
В steps сидят величины $m/v_{i,j}$, т.е. приращения чисел $x_{i,j}$ при переходе от $k$ к $k+1$.
В primes сидят значения малых простых, по модулю которых и производятся вычисления.
После сложения с приращением вычитаем из результата величину модуля (малого простого) и берём минимум из них, при этом если сумма превышает простое, то будет взято уменьшенное значение (т.е. как раз по модулю простого), если же сумма не превышает простое, то результат вычитания станет больше суммы и беззнаковый (unsigned) минимум возьмёт сумму. Именно в этом месте появляется ограничение на величину малых простых, чтобы сумма не переполняла отведённое поле и правильно работала команда беззнакового (unsigned) минимума.
Можно было бы сразу здесь же проверить результат на ноль в любой позиции и если он там есть, то переходить к следующему $k$, но из-за принятого решения заменить умножения $km$ на сложения $+m$ необходимо всегда подсчитывать все остатки $x_{i,j}$, даже если уже понятно что текущее $k$ не подходит.
Это обходится частичным разворачиванием цикла (loop unrolling) по $k$ по длине $l$ и дублированием таблиц в памяти $l$ раз. При этом можно сразу выкинуть из таблиц те $k$, которые точно не подойдут, например по модулю некоторых простых (если взять $l$ таким простым), что позволит ещё немного ускорить вычисления. При этом самый внутренний цикл проверяет лишь неравенство нулю чисел $x_{i,j}$ и при обнаружении нуля в любой позиции сразу же переходит к следующему $k$, не выполняя увеличения $x_{i,j}$, главное правильно инициализировать все таблицы, что выполняет внешняя программа, ускоритель же пользуется уже готовыми таблицами. В результате полное вычисление всех $x{i,j}$ выполняется в $l$ раз реже чем проверка их на нули. В таблицах при этом сидят значения для $k+0,k+1,k+2,...,k+l-1$ (если вдруг они все допустимы).
Итак, самый внутренний цикл проверки $x_{i,j}$ на нули:

Используется синтаксис ASM

.next:          inc             R15                             ;0156   1/0.25  k++
                add             RSI,const1                      ;0156   1/0.25
                cmp             RSI,stop2                       ;0156   1/0.25
                jnc             .add0                           ;06     1/0.5
                mov             EDI,0                           ;0156   1/0.25
.check:         vmovdqu         ymm15,[save_ymm+EDI]            ;23     3/0.5
                vpcmpeqw        ymm15,ymm15,[RSI+EDI]           ;15     1/0.5
                vptest          ymm15,ymm15                     ;0+5    2/1
                jnz             .next                           ;06     1/0.5   если нули есть, то что-то где-то разделилось и переходим к следующему
                add             EDI,32                          ;0156   1/0.25
                cmp             EDI,stop1                       ;0156   1/0.25
                jc              .check                          ;06     1/0.5

На самом деле значения $x_{i,j}$ проверяются не на нули, а на значения дающие ноль после суммирования с $km/v_{i,j}$, так как после разворачивания цикла (по $k$ на длину $l$, в RSI начало таблицы для соответствующего $k$) при каждом значении RSI числа в $x_{i,j}$ должны быть разными чтобы при суммировании давать ноль.

Полный кусок цикла поиска подходящего $k$:

Используется синтаксис ASM

;                                                               Ports   Time
.add0:          mov             EDI,0                           ;0156   1/0.25
.add:           vmovdqu         ymm15,[save_ymm+EDI]            ;23     3/0.5
                vpaddw          ymm15,ymm15,[steps+EDI]         ;15     1/0.5
                vpsubw          ymm14,ymm15,[primes+EDI]        ;15     1/0.5
                vpminuw         ymm15,ymm14,ymm15               ;15     1/0.5   текущее значение остатков по модулю простых
                vmovdqu         [save_ymm+EDI],ymm15            ;4      1/1
                add             EDI,32                          ;0156   1/0.25
                cmp             EDI,stop1                       ;0156   1/0.25
                jc              .add                            ;06     1/0.5
                mov             RSI,table1                      ;0156   1/0.25
                jmp             .check                          ;06     1/0.5
.next:          inc             R15                             ;0156   1/0.25  k++
                add             RSI,size_table1                 ;0156   1/0.25
                cmp             RSI,stop2                       ;0156   1/0.25  а не достигли ли конца разворота внутреннего цикла
                jnc             .add0                           ;06     1/0.5   если достигли, уходим на пересчёт всех x_i,j
                mov             EDI,0                           ;0156   1/0.25
.check:         vmovdqu         ymm15,[save_ymm+EDI]            ;23     3/0.5
                vpcmpeqw        ymm15,ymm15,[RSI+EDI]           ;15     1/0.5   проверка на (-km/v_i,j)%p_j - т.е. даст ли ноль после суммирования для текущего k
                vptest          ymm15,ymm15                     ;0+5    2/1
                jnz             .next                           ;06     1/0.5   если нули есть, то что-то где-то разделилось и переходим к следующему k
                add             EDI,32                          ;0156   1/0.25
                cmp             EDI,stop1                       ;0156   1/0.25
                jc              .check                          ;06     1/0.5   цикл по j проверки делимости на p_j
.ok:            ;здесь должен быть вывод найденного k
                jmp             .next

В table1 сидят значения $-km/v_{i,j}\pmod{p_j}$, т.е. те, которые дадут ноль при суммировании для текущего $k$.
Порядок расположения кусков кода специально подобран для лучшего предсказания условных переходов при наиболее вероятном пути выполнения.

Собственно это всё, остальное вынесено в генератор таблиц для асм кода.
Ну и пару мелких оптимизаций (кэширование в регистрах ymm первых значений save_ymm, разворачивание первого десятка итераций цикла .check, объединение нескольких vpcmpeqw в группу с одним vptest+jnz после неё, использование 8-ми битных простых вместо 16-ти битных) не показал так как они существенно усложняют код не давая принципиального выигрыша.

PS. Асм код правильный, а вот мои пояснения к нему могут быть не везде точными, особенно про индексы, кажется я в них запутался.
PPS. По скорости: наиболее часто выполняются несколько (1-2-3-5) итераций цикла .check, потом уход на .next, и лишь один из $l$ раз происходит уход на вычисление всех $x_{i,j}$ в .add0. В итоге, для 11-ти проверяемых чисел и порога малых простых до 3584 вместо 32768 средняя скорость составила 1.5 варианта $k$ за такт или примерно 5 млрд $k$ в секунду в одном потоке (thread).

[mathematician123]

Нашёл обсуждение уравнения $x^n - 2y^n = \pm 1$: https://mathoverflow.net/questions/257124/are-there-any-solutions-to-the-diophantine-equation-xn-2yn-1-with-x1-and. Там утверждается, что это уравнение не имеет решений при $x > 1$.

[Huz]

Dmitriy40 в сообщении #1556654 писал(а):

Huz
Hugo, похоже я так и не ответил Вам об устройстве ускорителей и возможности их применения для других паттернов (patterns). Отвечу здесь, возможно это будет интересно и другим.

Thank you very much. Google garbles the translation in parts, but I can understand much of it - I will need to study it further. I use quite a few of these techniques already (much slower, of course, in Perl), but I can see there is an opportunity to go further.

I clearly have some catching up to do if I want to take advantage of modern CPU instruction sets - the only book I have on 80x86 instructions was published in 1991.

[VAL]

Dmitriy40 в сообщении #1556654 писал(а):

Асм код правильный, а вот мои пояснения к нему могут быть не везде точными

Надо их аккуратно проверить. И Ваша часть статьи готова (с точностью до перевода). Завидую!

PS: T(36,11) <= 8807537694598756405759232050726160223249927395349997820
Как я уже писал, для 72 делителей без ускорителей реально построить цепочку еще, разве что, из 13 чисел.

[EUgeneUS]

mathematician123 в сообщении #1556664 писал(а):

Там утверждается, что это уравнение не имеет решений при $x > 1$.

Прекрасно! Там не просто утверждается, а даётся ссылка на работу Беннетта (2001), из которой это следует.

2 All

по цепочкам на 84 делителя. Вчера нашлась 8-ка - кандидатная в 11-ку. Но недостающие три числа раскладывались очень долго. Только сегодня утром разложились два и стало ясно, что там очередная 9-ка А третье число и за примерно 12 часов не разложилось.
К сегодняшнему вечеру досчитается до 400e105. Думаю, до конца недели дойду до 1000e105. А дальше посмотрим.
Похоже, с ростом чисел будут всё бОльшие проблемы с факторизацией.

[EUgeneUS]

mathematician123
Несколько комментариев:
1. В работе Беннетта (2001) есть ссылка на более раннюю работу (его с соавтором), от 1999 года. Там доказывается чуть менее общий результат, но он нам тоже подходит. (если не найдёте - скажите, я приведу со всеми ссылками, как до норамальной клавиатуры доберусь).

2. Применительно к теореме тут, закрываются все случаи $k \ge 3$, насколько понимаю. Но $k=2$ всё также нужно рассматривать отдельно.

[Dmitriy40]

VAL в сообщении #1556682 писал(а):

Dmitriy40 в сообщении #1556654 писал(а):

Асм код правильный, а вот мои пояснения к нему могут быть не везде точными

Надо их аккуратно проверить. И Ваша часть статьи готова (с точностью до перевода). Завидую!

Не, там более-менее понятно только до первого куска кода и о причине разворота цикла, а вот детали его разворота дальше очень сумбурно, я прямо сам читаю и не понимаю, хоть и в курсе что там и как. Впрочем для математической статьи дальше можно и не писать. Но и чуть выше переход от двухмерной матрицы проверяемых чисел и малых простых к одномерному суммированию векторов плохо объяснён. Так что прилизать ещё много где надо, всё же это писалось прямо из головы.

Huz в сообщении #1556674 писал(а):

I clearly have some catching up to do if I want to take advantage of modern CPU instruction sets - the only book I have on 80x86 instructions was published in 1991.

Hugo, now I checked all my text, Google quite clearly translated it into English, in some places the terms are incorrect (the replacement of "prime numbers" with something else is especially incomprehensible), but the context is clear.
About books, I use Intel official documentation and materials of Agner Fog, they are of course too complicated to get started, but there are a lot of books like "Introduction to AVX2" on the net, unfortunately I can’t recommend a specific one, but I think there are no frankly bad ones among them, for a start any will be enough, the main thing is to understand the principles of vector computing. It is not necessary to learn assembly language, Intel provides an "C/C++ Compiler Intrinsic Equivalents" for using AVX2 commands directly from C/C++.

EUgeneUS, mathematician123
Подскажите, а это доказательство получится расширить на случай произвольного количества сомножителей в $k$ или оно заточено только на $k=2pq$? А то в разложении $k$ на 10 случаев явно используется что там всего два простых, а потом появляются формулы типа $\gcd(p_i-1)>2$, в которых количество сомножителей уже никак не оговаривается ... Одно дело если это просто сокращение текста, другое если в этих местах привязки к двум простым не нужно.

-- 07.06.2022, 13:14 --

EUgeneUS в сообщении #1556683 писал(а):

Вчера нашлась 8-ка - кандидатная в 11-ку. Но недостающие три числа раскладывались очень долго. Только сегодня утром разложились два и стало ясно, что там очередная 9-ка А третье число и за примерно 12 часов не разложилось.

Да, вот поэтому и сделал факторизацию ограниченной, ведь не слишком редко встречаются и легко/быстро разложимые варианты, и так остаётся выбор ждать ли их или пытаться самому разложить трудные варианты до конца.

-- 07.06.2022, 13:37 --

Dmitriy40 в сообщении #1556337 писал(а):

VAL в сообщении #1556333 писал(а):

А как обстоят дела со 164-я, 172-я и 188-ю делителями?

164 считается, 5-ка и 6-ка найдены (за десяток минут и пару часов), а вот 7-ки так и нет (как и второй 6-ки кстати), уже больше суток. Похоже всё же идея с простыми в первой степени плоха. Для 172-х делителей от неё откажусь. А 164 пусть уж досчитывается как может, даже интересно когда же найдётся 7-ка.

За 90ч семёрки с 164 делителями так и не нашлось, вчера переделал паттерны обратно на 5 проверяемых чисел, запустил поиск семёрок, но что-то за 18ч только два не прошедших кандидата и всё, маловато. Правда 5-ка и 6-ка нашлись примерно вовремя.
Параллельно запустил и 172 делителя, но по паттернам с 7-ю проверяемыми числами, обратно переделывать пока не стал, за 4 суток нашлось под две сотни 5-ок и 7шт 6-ок, семёрки также нет.

[mathematician123]

EUgeneUS в сообщении #1556694 писал(а):

1. В работе Беннетта (2001) есть ссылка на более раннюю работу (его с соавтором), от 1999 года. Там доказывается чуть менее общий результат, но он нам тоже подходит. (если не найдёте - скажите, я приведу со всеми ссылками, как до норамальной клавиатуры доберусь).

Нашёл работу M. Bennett and B. M. M. de Weger, On the Diophantine equation $|ax^n-by^n| = 1$, Math. Comp. 67
(1998). Это она?

EUgeneUS в сообщении #1556694 писал(а):

2. Применительно к теореме тут
, закрываются все случаи $k \ge 3$, насколько понимаю. Но $k=2$ всё также нужно рассматривать отдельно.

В нашём проекте в Papeeria я написал доказательство. Там этот случай разобран. Ещё есть тема https://dxdy.ru/topic149661.html.

-- 07.06.2022, 13:53 --

Dmitriy40
Есть два доказательства. Первое: $M(2pq) \le 3$ для любых простых $p, q > 5$. Второе: если $d \equiv \pm 2 \pmod{12}$ и $gcd(p_i-1) > 2$, где $p_i$ - нечётные простые делители $d$, то $M(d) \le 3$.

-- 07.06.2022, 13:59 --

EUgeneUS
Есть ещё один вопрос. Как правильно: Михайлеску или Михэйлеску? Вы пишете через "а", в Википедии через "э". Спрашиваю, так как в моём доказательстве тоже есть ссылка на эту теорему.

[EUgeneUS]

mathematician123 в сообщении #1556702 писал(а):

Нашёл работу M. Bennett and B. M. M. de Weger, On the Diophantine equation $|ax^n-by^n| = 1$, Math. Comp. 67
(1998). Это она?

Да, верно. В работе Беннета (2001) ссылка на неё указана так (со страницами):

Цитата:

[BdW] M. Bennett and B. M. M. de Weger, On the Diophantine equation $|ax^n-by^n| = 1$, , Math. Comp. 67
(1998), 413-438.

Сорри, я ранее ошибся с годом.

-- 07.06.2022, 14:27 --

mathematician123 в сообщении #1556702 писал(а):

Есть ещё один вопрос. Как правильно: Михайлеску или Михэйлеску?

А я не знаю, как правильно и принято транслитирировать эту фамилию на русский.
Знаю, что в Румынии был король-комсомолец Михай I, поэтому пишу через "а", иногда вспоминаю, как написано на википедии - пишу через "э". Тут скорее уважаемый VAL может подсказать (насколько помню, он писАл, как в википедии).

-- 07.06.2022, 14:30 --

mathematician123 в сообщении #1556702 писал(а):

для любых простых $p, q > 5$

$p, q \ge 5$ (или так: $p, q > 3$)