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

[gris]

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

\\ алгоритм заключается в том, что на основе уже сформированных пирамид, то есть вектора всех чисел нечётной длины, являющихся квадратами вместе со всеми "квадратными" кольцами, формируется новый вектор из всех чисел длиной на 2 большей. Попутно очередная пирамидка проверяется на квадратность, но никак не отмечается, так как квадратность может случиться после надевания нескольких колец. Самый первый вектор [1,4,9]. Процесс продолжается, пока у компьютера будет хватать памяти. Необходимо хранить предыдущий вектор длиной  3*6^kl и новый вектор длиной в 6 раз большей.

vf=[1,2,3,4,6,8]; \\Первая цифра кольца
vt=[6,5,6,9,4,1]; \\Последняя цифра кольца
kp=3; \\ длина вычисляемой пирамидки
vp=[10,40,90]; \\ очередной вектор уже вычисленных пирамид. Удесятерённых. В самом начале — удесятерённые середины

forstep(kl=3,17,2, \\ длина вычисляемых новых пирамидок
  print("lenght=",kl); 
  kn=kp*6; vn=vector(kn); \\длина вычисляемого вектора и выделение места для него
  k=0;  \\ счётчик в этом векторе.
  p10=10^(kl-1); \\коэффициент умножения первой цифры внешнего кольца
  for(i=1,6, \\ цикл по добавляемым кольцам
    m=vf[i]*p10+vt[i]; \\"масса" кольца
    for(j=1,kp, \\цикл по уже вычисленному вектору пирамидок
      k++; \\счётчик
      vn[k]=m+vp[j]; \\добавление кольца. На каждую пирамидку надевается 6  колец
      if(issquare(vn[k],&s), print(vn[k]," = ", s,"^2") ); \\проверка на квадратность очередной новой пирамидки.
    ); \\ конец цикла по старому вектору
  ); \\ конец цикла по шести надеваемым кольцам
  vp=vn*10; \\ новый вектор пирамидок удесятеряется и становится старым.
  kp=kn; \\ и длина его
);\\ конец цикла по длинам пирамидок

Надеюсь, что вам еще не совсем надоело
Я пробовал перед проверкой issquare учитывать некоторые признаки неквадратности, но не заносить такое числа в вектор нельзя, разве что на последнем шаге. Мне бы хотелось всё-таки посмотреть на пирамидке длиной (высотой?) ещё не встреченной. Хотя основное это тренировки в улучшении алгоритмы и применении увы неосвоенного инструментария.

[wrest]

gris в сообщении #1695640 писал(а):

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

Проверки на неквадратность в стороне от issquare() не помогут в оптимизации, т.к. встроенная функция issquare() их и так делает, причём быстрее чем вы сможете реализовать в интерпретируемом коде. Там нет такого, что сразу вычисляется квадратный корень с плавающей точкой и потом проверяется не получилась ли равной нулю дробная часть, или извлекается целочисленный корень и обратным возведением в квадрат сравнивается с исходным. Там всё по-умному, сперва проверяются остатки от деления на несколько небольших модулей и только потом начинается извлечение корня. Можете просто закомментировать проверку issquare() или заменить функцию на просто ноль (всегда ложно) т.е. вместо if(issquare(),... написать if(0,...)и посмотреть, насколько ускорился перебор (который теперь ничего не найдёт). Такой метод позволяет понять куда на самом деле расходуется время (т.к. каких-то развитых профайлеров или нет или мы о них не знаем).

[wrest]

gris в сообщении #1695640 писал(а):

kn=kp*6; vn=vector(kn); \\длина вычисляемого вектора и выделение места для него

В итоге у вас должен получиться вектор из $3\cdot 6^9$ 19-значных чисел, если я правильно понял.
Да, тут есть заметное ускорение за счет переиспользования ранее сгенерированных кандидатов, но ценой невозможности вычислений для кандидатов длины больше 17 знаков...

[gris]

А вот пример.
Число 6314218_1_1659664=25128108^2 получилось после семи шагов
_1_
8_1_1
18_1_16
218_1_165
4218_1_1659
14218_1_16596
314218_1_165966
добавлением колец. Но ни одно из предшественников не является квадратом. Если бы я выкидывал эти числа, то не получил бы квадратную 19-пирамидку. Или я что-то не то делаю?
Хранить числа с предпредыдущего шага не нужно. Они сохранены.
У меня уже четвёртый час считается без сохранения сплошным перебором длина 23 для безнулевых и 19 для нулесодержащих Надоело, но жалко останавливать. Потом проверю время проверки на квадратность.

[wrest]

gris в сообщении #1695669 писал(а):

Надоело, но жалко останавливать.

Каждый шаг минимум в 5 раз дольше предыдущего. Уверены, что дождётесь?
В общем, ваша программа норм, но ест слишком много памяти, что препятствует росту высоты "пирамиды".
То есть, надо сокращать использование памяти, возможно в ущерб производительности.

Касательно генерации. Квадрат не может оканчивать на 5, т.к. это требует предпоследней цифры 2, а у нас таких в арсенале нет. Вот вам первый шаг к "математической" оптимизации: не генерируйте числа, у которых самое внешнее "кольцо" 25. Это минимум 1/6 всех генерируемых чисел.
Последним кольцом после 81 может быть 16, 36, 64 и не может быть 25, 49, 81
Последним кольцом после 25, 49, 81 может быть только 16 и не могут остальные
Вот в эту сторону надо бы смотреть, на мой взгляд.

[gris]

Расчётное время ещё часа два. Вот тут ещё одна проблема. Забыл, как внутри программы время выводить. Есть некий момент, скажем переход на следующий виток и надо "напечатать" время, которое затрачено начиная сначала. Не могу залезть в FM.

Увы, с нулями длину 19 считал 6 часов, без нулей 21 и23 тоже долго. Нет результатов :(

[wrest]

Новая версия
Ускоренная примерно в три раза.

1. find_special_squares(n_circles) ={ 
2.     /* Основные параметры */ 
3.     my(axis_digits = [1, 4, 9]);  /* Цифры оси*/ 
4.     my(circle_digits = [16, 25, 36, 49, 64, 81]);  /* Цифры колечек */ 
5.     /* Разделяем кольца на начало и конец */ 
6.     my(two_digit_pairs = vector(#circle_digits, i, [circle_digits[i]\10, circle_digits[i]%10])); 
7.     my(num=0); 
8.     my(count=0); 
9.      
10.     /* Проверяем только для заданной длины */ 
11.         print("Начинаю поиск для  ", n_circles, " колец"); 
12.          
13.         /* Готовим оси для соответствующей длины кандидатов */ 
14.         my(axis=vector(#axis_digits,i,axis_digits[i]*10^n_circles)); 
15.  
16.         /* Готовим кольца соответствующих диаметров*/ 
17.         my(circles=vector(n_circles,i,vector(#two_digit_pairs,j,two_digit_pairs[j][1]*10^(n_circles+i)+two_digit_pairs[j][2]*10^(n_circles-i)))); 
18.          
19.         /* Генерируем все комбинации колец */ 
20.         forvec(pairs = vector(n_circles, i, [1, #circle_digits]), 
21.                 /* Готовим первое кольцо */ 
22.                 num=circles[1][pairs[1]]; 
23.                 /* Нанизываем остальные кольца */ 
24.                 for(i = 2, n_circles, num+=circles[i][pairs[i]]); 
25.                 /* Каркас пирамиды готов */ 
26.                 for(c = 1, #axis, 
27.                            /* Проверяем получился ли квадрат c каждой из осей*/ 
28.                            if(issquare(num+axis[c]), print(num+axis[c], " = ", sqrtint(num+axis[c]), "^2");count++);  
29.                         ) 
30.                   ); 
31.        print("Найдено: ", count); 
32. } 
33. /* Запуск поиска */ 
34. find_special_squares(10); 

Ускорение за счёт замены вложенных циклов последовательными.
Предыдущая версия сначала вычисляла "ось" например 400 а потом "нанизывала" внешние "кольца", получая например 1960 и затем 31966. Заоем следующая ось и опять "нанизывание".
Теперь же сперва вычисляется 31066 а потом прибавляются "оси" 100, 400, 900.
Для например 6 "колец" по первому варианту получается 3х6=18 сложений а по второму 3+6=9 сложений.
Кроме того, я решил не делать цикл по каждой длине, а запускать только для одной конкретной длины.
21-значные числа "без нулей" программа на моём андроид-планшете проходит за 8 с четвертью минут:

? find_special_squares(10);
Начинаю поиск для  10 колец
Найдено: 0
time = 8min, 13,430 ms.
?

Предоставляю вам запуск для 11 (ожидаемое время на моем планшете 48 минут), 12 колец (ожидаемое время 5 часов).

[Dmitriy40]

Имеет смысл нанизывать кольца снаружи, а не изнутри: так можно проверять что младшие цифры могут быть получены возведением некоего числа в квадрат.
Так можно уменьшить количество кандидатов для полной проверки и не перебирать внутренние кольца если внешние уже недопустимы.
Выигрыш:

d=[16,25,36,49,64,81];\\Кольца
dl=[t%10|t<-d];\\Младшие цифры колец
for(k=2,9,\\Количество внешних колец
   n=0;
   forvec(x=vector(k,i,[1,#d]),
      if(!issquare(Mod(fromdigits([dl[t]|t<-x]),10^#x)), next);
      n++;
   );
   print("6^",k,"=",6^k," / ",n," = ",6.0^k/n);
);
print(strtime(gettime()));

Запуск:
6^2=36 / 15 = 2.4000000000000000000000000000000000000
6^3=216 / 72 = 3.0000000000000000000000000000000000000
6^4=1296 / 270 = 4.8000000000000000000000000000000000000
6^5=7776 / 1377 = 5.6470588235294117647058823529411764706
6^6=46656 / 6804 = 6.8571428571428571428571428571428571428
6^7=279936 / 38637 = 7.2452830188679245283018867924528301887
6^8=1679616 / 218700 = 7.6800000000000000000000000000000000000
6^9=10077696 / 1292517 = 7.7969543147208121827411167512690355330
20,950 ms

Разумеется можно и интерактивно накапливать кольца от внешних к внутренним с оставлением только допустимых вариантов (1.3млн вместо 10млн для 19-ки). Тогда при нахождении допустимого варианта подставить три оси и проверить не решение ли и если да, то вывести и продолжить перебор.

[wrest]

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

Имеет смысл нанизывать кольца снаружи, а не изнутри: так можно проверять что младшие цифры могут быть получены возведением некоего числа в квадрат.

Хорошая мысль! Будет зависеть от скорости fromdigits(), надо попробовать.
Я вообще вот не уверен что строковые манипуляции с последующим fromdigits() медленнее арифметики pari/gp, заточенной под произвольную точность.

[Dmitriy40]

Вот реализация идеи, выводит все возможные пирамиды вплоть до k_max колец:

k_max=10;\\Сколько максимум колец проверять

c=[1,4,9];\\Ось
d=[16,25,36,49,64,81];\\Кольца
dh=[(t\10)%10|t<-d];\\Старшие цифры колец
dl=[t%10|t<-d];\\Младшие цифры колец

{Piramida(x)=my(y,z,t,k=#x);
   for(i=1,#d,
      y=10^k*dl[i]+fromdigits([dl[t]|t<-x]);
      if(!issquare(Mod(y,10^(k+1))), next);
      for(j=1,#c,\\Проверка не получилась ли пирамида с одной из осей
         z=fromdigits([dh[t]|t<-Vecrev(x)])*10^(k+3)+dh[i]*10^(k+2)+10^(k+1)*c[j]+y;
         t=0; if(issquare(z,&t), print(z," = ",t,"^2"); );
      );
      if(k+1<k_max, Piramida(concat([i],x)); );\\Продолжим нанизывать кольца
   );
}

Piramida([]);\\Начинаем с пустой пирамиды
print(strtime(gettime()));

Запуск для 10 колец:
196 = 14^2
38416 = 196^2
631421811659664 = 25128108^2
841 = 29^2
1min, 25,940 ms

Запуск для 11 колец:
--/--/--
9min, 6,455 ms

Думаю можно и ещё ускорить заменив оба fromdigits на передачу параметрами снаружи, ведь там они и так нужны для проверки, т.е. лишнего времени на вычисления не нужно.

[Dmitriy40]

Заменил, ускорилось:

c=[1,4,9];\\Ось
d=[16,25,36,49,64,81];\\Кольца
dh=[(t\10)%10|t<-d];\\Старшие цифры колец
dl=[t%10|t<-d];\\Младшие цифры колец

{Piramida(k_max,hi,lo,k)=my(y,z,t,q,kk=10^(k+1));
   for(i=1,#d,
      y=10^k*dl[i]+lo;
      if(!issquare(Mod(y,kk)), next);
      q=hi*100+dh[i]*10*kk;
      for(j=1,#c,
         z=q+c[j]*kk+y;
         t=0; if(issquare(z,&t), print(z," = ",t,"^2"); );
      );
      if(k+1<k_max, Piramida(k_max,q,y,k+1); );
   );
}

Piramida(11,0,0,0);\\Начинаем с пустой пирамиды, до 11 колец
print(strtime(gettime()));

Запуск для 10 колец:
196 = 14^2
38416 = 196^2
631421811659664 = 25128108^2
841 = 29^2
42,791 ms

Запуск для 11 колец:
--/--/--
4min, 16,949 ms

[wrest]

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

Заменил, ускорилось:

Прекрасно, я считаю! Добавил счётчики

? {Piramida(k_max,hi,lo,k)=
   my(y,z,t,q,kk=10^(k+1));
   for(i=1,#d,
      y=10^k*dl[i]+lo;
      count3++; /* отработало 9 350 394 */
      if(!issquare(Mod(y,kk)), count1++;next); /* отказ 154 992 */
      q=hi*100+dh[i]*10*kk;
      for(j=1,#c,
         z=q+c[j]*kk+y;
         t=0; if(issquare(z,&t), print(z," = ",t,"^2"); );
      ); count2++; /* отработало 9 195 402 */
      if(k+1<k_max, Piramida(k_max,q,y,k+1); );
   );
}
? Piramida(10,0,0,0);
196 = 14^2
38416 = 196^2
631421811659664 = 25128108^2
841 = 29^2
time = 1min, 42,391 ms.
? count1
%36 = 154992
? count2
%37 = 9195402
? count3
%38 = 9350394
? 3*6^10
%39 = 181398528
?

Без отсева, для 10 колец надо проверить

? total=0;for(i=1,10,total+=3*6^i);total
%40 = 217678230
?

Так что математиака тут сыграла бОльшую роль чем прогаммирование.
По программированию. Ну лично мне не нравятся функции которые не содержат в себе всё что им надо.
Тут рекурсия срабатывает (для 10 колец) 1,5 млн. раз. Это довольно сильно потребляет память на стек. У меня на планшете, например, становится плохо браузеру. Хотя казалось бы - в андроиде изоляция, управление памятью неплохое, да и памяти в планшете 12 гигов при 6 свободных (parisizemax установлен в 2 гигабайта)...

[wrest]

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

Это довольно сильно потребляет память на стек. У меня на планшете, например, становится плохо браузеру.

Хотя... может это не браузеру плохеет а глюки dxdy. Так что этот момент снимаю с повестки.

Но лично я бы предпочёл нерекурсивные (в смысле явной рекурсии - вызова функцией самой себя) функции там, где глубина рекурсии может быть очень большая.

[wrest]

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

Имеет смысл нанизывать кольца снаружи, а не изнутри:

Гениальная всё-таки мысль. Браво!

[Dmitriy40]

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

-- 29.07.2025, 17:49 --

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

где глубина рекурсии может быть очень большая.

Здесь глубина рекурсии ограничена количеством нанизываемых колец, два десятка уровней это смешно.

-- 29.07.2025, 17:54 --

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

Тут рекурсия срабатывает (для 10 колец) 1,5 млн. раз. Это довольно сильно потребляет память на стек.

Памяти слёзы (сотни байтов на вызов), а вот трафик памяти - да. Но вряд ли тормозит именно передача параметров ... всё же внутри достаточно сложных вычислений, основное время должны занимать они.
Кажется Вы перепутали глубину рекурсии и общее количество вызовов функции.

-- 29.07.2025, 18:17 --

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

Гениальная всё-таки мысль.

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