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

[gris]

Я начал изучение ТЧ по учебнику Александра Бухштаба и хочу через любимую PARI/GP прорабатывать материал практически. Уже получил интересные примеры и тщательно их изучаю. Надеюсь, что и в дальнейшем уважаемые известно кто не откажутся дать мне полезный совет и укажут на ошибки.
Пока немного забавного. В первой главе "Простые Числа" задумался над фразой в доказательстве Теоремы 22:
"Таким образом, каждое из этих чисел имеет положительного делителя, отличного от 1 и самого себя". Сразу вспомнилось про подвалов и шкафов . А делители ведь одушевлённые конструкты! Я это всегда чувствовал.

[wrest]

gris
потому что имеет "что", а не имеет - "чего".
"имеет делитель" но "не имеет делителя"

[gris]

wrest, согласен! Вася имеет (в хор. см.) кота, а Серж не имеет кота. Даже, если кот мёртвый. Но Вася имеет компьютер, а Серж не имеет компьютера.
Но там речь о существовании вектора любой длины из последовательных составных чисел. То есть речь идёт о составных числах, которые-таки имеют собственного делителя. И поделом ему Не подумайте, что я выискиваю какие-то особенности. Просто внимательно читаю и обдумываю.
А вот отвлекусь на минутку и сооружу понятно какую последовательность.

{\\Первое появление к составных чисел подряд
d=vector(1000); 
d[1]=4;
p=5;
forprime(pn=7,1000,
  dp=pn-p-1;
  if(d[dp],,d[dp-1]=p+1;d[dp]=p+1);
  p=pn;
);
print(d[1..13]);
}
[4, 8, 8, 24, 24, 90, 90, 140, 140, 200, 200, 114, 114]

Всё-таки $24=4!$ Факториал используется в доказательстве.
Конечно, можно заметить, что чётных провалов между простыми и не бывает, ну кроме нуля между 2 и 3.

+++ Yadryara, я рад, что вы заглянули. У вас я научился слову гэп и даже хотел пересчитать всё на первое появление гэпа длины к. В моём учебнике пока не нашёл про Платта. Но есть же поговорка, что простых больше, чем звёзд на небе! Как же можно их количество посчитать? Разве что число простых хоккейных голов

[Yadryara]

gris, очень рад, что Вы эту тему внимательно изучаете. Может поймёте как например у Платта посчитано количество простых чисел.

[wrest]

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

Может поймёте как например у Платта посчитано количество простых чисел.

Там же лютый матан с интегральными логарифмами, остаточными членами во всевозможных формах и вот этим всем, а не стройная чистая ТЧ

[Yadryara]

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

Там же лютый матан с интегральными логарифмами, остаточными членами во всевозможных формах и вот этим всем, а не стройная чистая ТЧ

Ну так Вы разобрались в этом лютом матане? Хотя бы частично?

Если правильно помню, vicvolf говорил что там доходчиво написано и надо как раз читать Бухштаба.

[Yadryara]

Только сейчас заметил:

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

У вас я научился слову гэп

У меня? А как же знаменитая Prime gaps — A001223.

Почему аномально большой гэп 14 появляется после рекордного 8 ? Сумма по дзета-нулям увеличивается ещё и в точке 121 (квадрат простого), и в точке 125 (куб простого), и в точке 128 (7-я степень простого). Интересно как это работает?

[gris]

Ах, если бы это было правильное Захарово, но, увы, только симметричное. Подсказать было некому... Я далеко... И вместо памятника ВИЛ памятник юному АСП. Ну да ладно.
В общем, поиск квадратных пирамидок. Они пришли из сопредельных пространств.
Что такое квадратная пирамидка? Это ось: 0, 1, 4, 9. И квадратные надетые кольца: 16, 25, 36, 49, 64, 81. Например: 196. Легко приходит на ум. Дальше 841 и 38416. Ведущие нули запрещены и лучше вообще без нулей.
631421811659664 = 25128108^2. Каждое колечко — квадрат. И само число квадрат. Я программку написал, но уж больно долго работает.

\\квадратные пирамидки
{mdg=15;
vm=[0,1,4,9];
vf=[0,0,0,0,1,2,3,4,6,8];
vt=[0,1,4,9,6,5,6,9,4,1];
for( m=1,mdg, print("l=",2*m+1);
  ww=vector(m,i,[1,10]); ww[m]=[5,10];
  foreach( vm,md, 
    nm=md*10^m;
    forvec( w=ww,
      s=nm;
      for( k=1,m,
        s+=vf[w[k]]*10^(m+k)+vt[w[k]]*10^(m-k)
      ); 
      if(issquare(s),print(s," ", sqrtint(s),"^2"));
    );
  );
);
}

l=3
841 29^2
196 14^2
l=5
38416 196^2
l=7
8300161 2881^2
6061444 2462^2
l=9
380094016 19496^2
l=11
66816046144 258488^2
l=13
3003129566116 1732954^2
l=15
631421811659664 25128108^2
l=17
41260004094094569 203125587^2

Народу понравилось: много различных циклов. Хотел вернумшись улучшить, но пока другие заботы
Дайте совет!

[wrest]

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

Каждое колечко — квадрат.

Моих телепатических способностей не хватило для расшифровки вашего желания. Колечки, оси, пирамидки...

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

Дайте совет!

Совет: пишите яснее.

Обратился к ИИ. Его вердикт такой
Код ищет **симметричные числа, являющиеся квадратами**, но делает это полным перебором. Его можно ускорить в **10–100 раз**, если:
1. Сузить перебор, используя свойства квадратов.
2. Заменить генерацию+проверку на целевой подбор `k²` с нужной структурой цифр.
3. Добавить модульные фильтры и параллельные вычисления.

[gris]

У нас есть полные натуральные квадраты из одной десятиричной цифры: 0, 1, 4, 9. И квадраты, состоящие из двух таких цифр: 00, 01, 04, 09, 16, 25, 36, 49, 64, 81.
Берём одноцифренный квадрат в качестве центра. И на него начинаем последовательно нанизывать двухцифренные квадраты: первая цифра впереди, вторая сзади. Например:
$ 4 \to 040, 041, 044, 049, 146, 245, 346, 449, 644, 841 \to 00400, 00401 ... 88411$
Каждое получающееся число проверяем нет ли ведущего нуля (они могут потом пригодиться) и является ли оно полным квадратом. Таких очень мало. Но надо найти побольше. Вот первые пирамидки (не я так назвал, честно!)
841, 196, 38416, 8300161, 6061444, 380094016...
само число квадрат и каждая симметричная относительно центра пара цифр составляет квадрат. И сам центр квадрат. Иногда и другие квадраты просматриваются. Проблема в организации перебора.
ИИ немножко не про то. Я для контроля пробовал и полый перебор, но это после Е12 очень долго. Но результаты совпадают. Наверное, надо выращивать дерево, чтобы не применять постоянно fromdigits.
Для меня сами числа не важны, а важен алгоритм и его реализация. Можно брать числа из чётного количества цифр или искать кубы.
Вот куб с тремя кубичными кольцами: 73 6 754 6 2 8 7 4 480 4 59 = 419219^3: 8, 27, 64, 64. Увы, из девяти.

[Yadryara]

(gris)

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

Ах, если бы это было правильное Захарово, но, увы, только симметричное. Подсказать было некому... Я далеко... И вместо памятника ВИЛ памятник юному АСП.

Это просто блестяще

Давненько я не писывал текста, который может понять лишь один форумчанин

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

Обратился к ИИ.

Лучше спросите у него, что значит АСП. До Пушкина он может и догадается, но скажу по секрету, у меня это было сокращение от АнтиСпецПотрох.

[wrest]

gris
Вот без нулей и немного быстрее

find_special_squares(max_length) ={
    /* Основные параметры */
    my(axis_digits = [1, 4, 9]);  /* Цифры оси*/
    my(circle_digits = [16, 25, 36, 49, 64, 81]);  /* Цифры колечек */
    my(two_digit_pairs = vector(#circle_digits, i, [circle_digits[i]\10, circle_digits[i]%10]));
    

    /* Основной цикл */
    for(n_circles = 1, max_length,
        print("\nПоиск для длины ", 2*n_circles+1);
        
        /* Генерируем все комбинации пар */
        forvec(pairs = vector(n_circles, i, [1, #two_digit_pairs]),
            /* Перебираем центральные цифры */
            for(c = 1, #axis_digits,
                /* Кладём ось */
                my(num = axis_digits[c]);
                /* Нанизываем пары цифр */
                for(i = 1, n_circles,
                    my(pair = two_digit_pairs[pairs[i]]);
                    num=num*10+pair[1]*10^(2*i) + pair[2]
                 );
                                /* Проверяем получился ли квадрат*/
                if(issquare(num), print(num, " = ", sqrtint(num), "^2"))
            )
        )
    )
;
}
/* Запуск поиска */
find_special_squares(8);

Работает так:

? find_special_squares(8);

Поиск для длины 3
196 = 14^2
841 = 29^2

Поиск для длины 5
38416 = 196^2

Поиск для длины 7

Поиск для длины 9

Поиск для длины 11

Поиск для длины 13

Поиск для длины 15
631421811659664 = 25128108^2

Поиск для длины 17
time = 1min, 14,877 ms.
?

[wrest]

В общем пока не забыл, есть у меня такие мысли по оптимизации.
1. Несмотря на то, что кажется, что проверка на квадрат затратная, это не так. Встроенная issquare() работает чрезвычайно быстро. На неё уходит не более 5% времени.
2. Основное время уходит на генерацию чисел для проверки. Но квадратов, в диапазоне $10^{17}..10^{18}$ всё равно больше чем чисел, которые надо проверить, поэтому перебор квадратов вместо перебора кандидатов в квадраты, скорее всего, не поможет.
3. При генерации чисел происходит умножение на степени десятки. Вместо этого, лучше предвычислить то, что нужно прибавить для "нанизывания" кольца. Например, чтобы "нанизать" 16 на 146, нужно 146 умножить на 10 и прибавить 10006 -- вот эти 10006, 20005 и т.п. надо предвычислить и просто прибавлять.
4. Надо проверить, что возможно работа со строками будет быстрее арифметики. Но тут я не уверен в целесообразности.
5. Для 3 "осевых" цифр [1, 4, 9] и 6 "колец" без нулей [16, 25, 36, 49, 64, 81], количество комбинаций для 8 надетых колец (т.е. 17-значных чисел) составляет $3\cdot6^8=5,038,848$. А для 9 колец уже $3\cdot 6^9=30,233,088$ Тут надо думать как это количество вариантов радикально уменьшить, чтобы даже не генерировать числа, которые точно не могут быть квадратами. Это не вопрос программирования, а вопрос математический.

[gris]

wrest, спасибо. Я слегка уже думал, как это дело ускорить. Но другие напасти навалились. За две недели моего отсутствия комп обленился и стал отключаться внезапно и постоянно. Возможно от жары. Мышь капризничает: у неё хвост из трёх частей 6 метров. Клавиатура работает хорошо, но она сама хорошая, а мыши дешёвые и падают. Моник разбитый. В общем, надо всё радикально менять. И поставить наконец 64-битное(-ую, -ый) PARI/GP. В общем, такие смешные проблемы.
Вчера так и не дождался окончания обсчёта 19-значных . Но я пока всерьёз и не ковырялся в алгоритме. С наскока получил пятнашку (какие, однако, параллели ) и успокоился. Я ещё подумаю над всем этим.

+++++++++++
wrest, учёл ваши советы и слепил снегурочку
Но она растаяла на длине 19

(12:29) gp > {
vf=[1,2,3,4,6,8];
vt=[6,5,6,9,4,1];
kp=3;
vp=[10,40,90];
forstep(kl=3,21,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];
      if(issquare(vn[k],&s), print(vn[k]," = ", s,"^2") );
    );
  );
  vp=vn*10;
  kp=kn;
);
}
lenght=3
196 = 14^2
841 = 29^2
lenght=5
38416 = 196^2
lenght=7
lenght=9
lenght=11
lenght=13
lenght=15
631421811659664 = 25128108^2
  *** _*_: Warning: increasing stack size to 32000000.
lenght=17
  *** _*_: Warning: increasing stack size to 256000000.
lenght=19
  ***   at top-level: ..."lenght=",kl);kn=kp*6;vn=vector(kn);k=0;p10=10
  ***                                             ^---------------------
  *** vector: overflow in lg().
(12:30) gp >

Никак нет научусь использовать процедуры.
Памяти не хватает? Очистка предыдущего не поможет, наверное, да я и не умею. Может быть попробовать map?
Продолжил с предварительно рассчитанного 17 до 19 и 21 на лету без записи вектора. За 8 минут ничего не нашлось. Попробую 23...

[wrest]

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

Может быть попробовать map?

map с точки зрения контроля использования памяти в pari/gp недоработанная немного, насколько я помню.
При переполнении pari/gp просто падает с segmentation fault.
Вот вам реализация пункта

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

3. При генерации чисел происходит умножение на степени десятки. Вместо этого, лучше предвычислить то, что нужно прибавить для "нанизывания" кольца. Например, чтобы "нанизать" 16 на 146, нужно 146 умножить на 10 и прибавить 10006 -- вот эти 10006, 20005 и т.п. надо предвычислить и просто прибавлять.

1. find_special_squares(max_length) ={ 
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.      
8.     /* Основной цикл */ 
9.     for(n_circles = 1, max_length, 
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.             for(c = 1, #axis, 
23.                 /* Кладём подготовленную ось */ 
24.                 my(num=axis[c]); 
25.                 /* Нанизываем подготовленные кольца */ 
26.                 for(i = 1, n_circles, 
27.                     num+=circles[i][pairs[i]]; 
28.                  ); 
29.                                 /* Проверяем получился ли квадрат*/ 
30.                if(issquare(num), print(num, " = ", sqrtint(num), "^2"));  
31.  
32.             ) 
33.        ); 
34.     ) 
35. ; 
36. } 
37. /* Запуск поиска */ 
38. find_special_squares(8); 

Работает так:

? find_special_squares(8);
Начинаю поиск для  1 колец
196 = 14^2
841 = 29^2
Начинаю поиск для  2 колец
38416 = 196^2
Начинаю поиск для  3 колец
Начинаю поиск для  4 колец
Начинаю поиск для  5 колец
Начинаю поиск для  6 колец
Начинаю поиск для  7 колец
631421811659664 = 25128108^2
Начинаю поиск для  8 колец
time = 37,510 ms.
?

То есть 17-значные перебирает за 38 секунд.
Ускорение примерно в 2 раза по сравнению с предыдущим кодом с вычислением степеней десятки внутри циклов "нанизывания".
Ваш код непросто разбирать, почему у вас там памяти не хватает не знаю. Пишите комментарии!