PARI/GP и Теория Чисел

[wrest]

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

Кажется Вы перепутали глубину рекурсии и общее количество вызовов функции.

Ну наверное. Просто вызов функции всегда, кроме первого раза, делается из самой функции. Соответственно, сохраняется состояние вызвавшей функци, создаются новые экземпляры всей внутренности функции (всё что в my())и т.п.
Ну может я просто не понимаю...

-- 29.07.2025, 19:05 --

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

Этот метод (проверка по модулю, младших цифр) вообще довольно часто применим.

Конечно, более того и в issquare() эти проверки делаются по нескольким модулям (16 и 64 емнип)
Не, ну задним-то умом оно конечно очевидно, когда уже сделано и работает.
Но надевать кольца снаружи это прям огонь! Меня вот не осенило.
Саме мысль мне приходила, более того я был уверен что как-то отрезать минимум на порядок объем можно. Я пробовал и проверять по модулю, но уже почти составленную пирамиду (т.е. тратил время на строительство), и получил жалкие проценты ускорения. А тут - радикально, на два порядка!
Урок мне и ТСу (надеюсь) хороший.

-- 29.07.2025, 19:16 --

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

За 16ч посчитались все пирамиды по 29 включительно (по 14 колец, k_max=14). Новых решений не найдено!

Ну вот на этом месте пора перейти к поиску доказательства что решений больше и нет :)

[Dmitriy40]

wrest в сообщении #1695784 писал(а):

Конечно, более того и в issquare() эти проверки делаются по нескольким модулям (16 и 64 емнип)

Забавно что на асме проверка по модулю 64 делается буквально одной командой - bt <маска_в_x64_регистре>,<младшие_64_бита_проверяемого_числа_в_x64_регистре>, потом сразу j(n)c next (jump) и привет. Под x32/x86 тоже работает, с маской 32 бита. Если нужен модуль больше 64 (но степень двойки чтобы команда деления с остатком не испортила всю малину), то можно маску хранить и в памяти, тогда она и модуль вплоть до $2^{64}$ битов могут быть (если памяти хватит!), но это лишний произвольный (что плохо, они отвратительно кэшируются) доступ к памяти.

[gris]

Несмотря на большую загруженность посторонними делами в конце месяца, я читаю и стараюсь проникать в подробности. И тщательно разбираться в текстах программ. Даже по мелочам, типа gettime(). Спасибо. Кстати, закомментил проверку на квадратность и занимет она 8% времени (для размеров 9-17. Идею внешних колец начал было пробовать. Хотел даже предвычислять все толстые (до 8 слоёв) кольца анализируя посторонние цифры в концах квадратов. Например, для двойного получается [16, 41, 44, 49, 56, 61, 64, 69, 96], а начало определяется, хоть и неоднозначно из-за 6. А потом впихивать середину. Запутался. У меня для безнулёвок получилось только для 25 и 27. А нулёвки на 19 забросил. А ведь ещё и пирамидки с чётным количеством цифр. Очень это похоже на кортежи! Отличие в конечности перебора .

[maxal]

Вот код, который фиксирует суффикс и префикс длины $s$ искомого квадрата и использут функцию sqrtintmod2(N,q) и параллелизацию для эффективного нахождения таких квадратов. Если ищем квадраты длины $L$, то $s$ нужно оптимально брать равным примерно $L/4$.

export(isNULL);
export(asqrtintmod2);

{ piramid_len(L) = my(rings, fd, ld, s);
    rings = vector(10, i, (i-1)^2);
    fd = rings \ 10; \\ first digits
    ld = rings % 10; \\ last digits
    s = max(1,L\4);
    parforvec(v=vector(s,i,[if(i<s,1,5),#ld]),
        my(suff, pref);
        suff = fromdigits(apply(x->ld[x],v));
        pref = fromdigits(Vecrev(apply(x->fd[x],v)));
        foreach(asqrtintmod2(10^s,suff), r,
            forstep(t=sqrtint(pref*10^(L-s)-1)+1, sqrtint((pref+1)*10^(L-s)-1), r,
                my( d = digits(t^2) );
                if(!issquare(d[(L+1)\2]),next);
                for(i=s+1,(L-1)\2, if( !issquare(d[i]*10+d[L+1-i]), next(2) ); );
                print([pref,suff]," ",r," ",t," ",t^2);
            );
        );
    );
}

Результаты такие:

? piramid_len(3)
[1, 6] Mod(4, 10) 14 196
[8, 1] Mod(9, 10) 29 841
? piramid_len(5)
[38, 16] Mod(46, 50) 196 38416
? piramid_len(7)
[60, 44] Mod(12, 50) 2462 6061444
[83, 61] Mod(31, 50) 2881 8300161
? piramid_len(9)
[380, 16] Mod(496, 500) 19496 380094016
? piramid_len(11)
[668, 144] Mod(488, 500) 258488 66816046144
? piramid_len(13)
[3003, 6116] Mod(454, 2500) 1732954 3003129566116
? piramid_len(15)
[6314, 9664] Mod(608, 2500) 25128108 631421811659664
? piramid_len(17)
[41260, 94569] Mod(25587, 50000) 203125587 41260004094094569
? piramid_len(19)
? piramid_len(21)
? piramid_len(23)
? piramid_len(25)
[3003260, 4456196] Mod(1635486, 2500000) 1732991635486 3003260008664441094456196
? piramid_len(27)
[1208060, 401156] Mod(1237534, 2500000) 10991178737534 120806010040419494060401156
? piramid_len(29)
[3001800, 4116996] Mod(2018114, 2500000) 173257034518114 30018000010010946096194116996
[6620468, 1494544] Mod(5612, 625000) 257302705005612 66204682003204990560951494544
? piramid_len(31)
? ##
  ***   last result: cpu time 50min, 15,050 ms, real time 7min, 45,920 ms.

Как видим, поиск квадратов длины 31 потребовал около 8 минут. Так что, можно поискать и дальше, если есть желание.

[wrest]

maxal в сообщении #1698588 писал(а):

Вот код, который

У меня на планшете не завелось :(

? export(asqrtintmod2);
? piramid_len(11)
***   at top-level: piramid_len(11)
***                 ^---------------
***   in function piramid_len: ...ly(x->fd[x],v)));foreach(asqrtintmod2(10^s,suf
***                                                        ^---------------------
***   not a function in function call
? \v
                                                    GP/PARI CALCULATOR Version 2.17.2 (released)
                                     aarch64-linux running android (portable C/GMP-6.3.0 kernel) 64-bit version
compiled: Mar  6 2025, Android (12470979, +pgo, +bolt, +lto, +mlgo, based on r522817c) clang version 18.0.3 (https://android.googlesource.com/toolchain/llvm-project d8003a456d14a3deb8054cdaa529ffbf02d9b262)
                                                              threading engine: single
                                                   (readline v8.1 enabled, extended help enabled)
?

[maxal]

wrest в сообщении #1698608 писал(а):

У меня на планшете не завелось :(

Так нужно сначала функцию asqrtintmod2 определить - я дал ссылку выше.

[wrest]

maxal в сообщении #1698613 писал(а):

Так нужно сначала функцию asqrtintmod2 определить - я дал ссылку выше.

Получилось.
Считаю "без нулей", т.е.
rings = vector(6, i, (i+3)^2);
И в один поток (forvec вместо parforvec).
? piramid_len(31)
time = 15,945 ms.
? piramid_len(33)
time = 14,226 ms.
? piramid_len(35)
time = 1min, 44,231 ms.
? piramid_len(37)
time = 1min, 30,985 ms.
?
Удивительно...

[maxal]

wrest в сообщении #1698617 писал(а):

Считаю "без нулей", т.е.
rings = vector(6, i, (i+3)^2);

Если без нулей, то еще if(i<s,1,5) нужно заменить на просто 1.

-- Mon Aug 18, 2025 04:34:48 --

А также после my( d = digits(t^2) ); добавить проверку на наличие нулевых цифр:
if(!vecmin(d),next);.

[maxal]

maxal в сообщении #1698588 писал(а):

Как видим, поиск квадратов длины 31 потребовал около 8 минут. Так что, можно поискать и дальше, если есть желание.

Немного соптимизировал и дошел до длины 35:

? piramid_len2(31)
? ##
  ***   last result: cpu time 20min, 18,319 ms, real time 2min, 38,044 ms.
? piramid_len2(33)
[34010803, 64196996] Mod(9315614, 25000000) 18442018059315614 340108030100123245656446064196996
? ##
  ***   last result: cpu time 20min, 13,936 ms, real time 2min, 40,059 ms.
? piramid_len2(35)
? ##
  ***   last result: cpu time 3h, 34min, 23,487 ms, real time 27min, 54,997 ms.

[wrest]

maxal
Но вот что интересно: время растёт скачками. 15 и 16 "колец" считаются одинаковое время - меньше трех минут, 17 колец считаются 27 минут. У меня (в версии без нулей и без предложенных вами поправок) получилось так же по соотнолениям. И 18 колец посчитались даже быстрее чем 17.

[maxal]

wrest в сообщении #1698791 писал(а):

Но вот что интересно: время растёт скачками. 15 и 16 "колец" считаются одинаковое время - меньше трех минут, 17 колец считаются 27 минут. У меня (в версии без нулей и без предложенных вами поправок) получилось так же по соотнолениям. И 18 колец посчитались даже быстрее чем 17.

15 и 16 колец соответствуют длинам 29 и 31, которые дают одно и тоже $s = 7$. Следующая длина дает $s=8$, откуда и скачок. Можно поиграться со значениеми $s$, так как выбор $s=[L/4]$ основан на грубой прикидке.

Замена if(i<s,1,5) на 1 необходима, иначе вы теряете решения, и замеры времени не имеют смысла.