интерактивный курс: введение в программирование на PARI/GP

[Dmitriy40]

Не помню было ли, но пусть даже и повторю.
Иногда бывает удобно/нужно изменить одновременно две (три и т.д.) переменные, например в рекуррентных формулах типа $(x_{i+1},y_{i+1})=(2x_i+3y_i,3x_i+4y_i)$, это можно сделать без введения новых переменных: [x,y]=[2*x+3*y,3*x+4*y].
Справа может стоять и массив или его часть, слева тоже можно массив, но только целиком.

[Dmitriy40]

Ещё небольшой лайфхак: многие функции принимают лишь два аргумента (gcd, lcm, chinese, и т.д.) и в то же время принимают и вектор значений, например:

? lcm([2,3,5,7])
%1 = 210
? chinese([Mod(1,2),Mod(1,3),Mod(3,7)])
%2 = Mod(31, 42)

[gris]

А вот возникла потребность в лайфхаке.
Предположим, некто хочет организовать сплошной перебор по парам натуральных чисел для проверки некоторой однородной функции. То есть результат проверки будет одинаков для пар $(n,m)\sim (kn,km) \;\forall\;n,m,k$
Например, такая последовательность пар:
$(1,1),(1,2),(2,1),(1,3),(3,1),(1,4),(2,3),(3,2),(4,1)...$ по возрастанию суммы элементов. Не нужно проверять пары $(2,2),(2,4),(4,2),(3,3),(6,2),(6,3)$Что бы такое придумать, если проверка функции очень длительная. Нет ли в PARI/GP чего-нибудь этакого?
Типа просто проверки на то, что у пары есть общий множитель?
А если это тройки? Например, для проверки ВТФ

[Dmitriy40]

gris
Вот как раз gcd([...]) прекрасно справится. Если я правильно понял.

? gcd([125,210,75,195])
%1 = 5

Для пары можно и без массива.

Но если перебор долгий, то лучше побольше вариантов исключить руками (да хоть проверкой что оба числа делятся на 2,3,5), не думаю что gcd() сильно быстрая.

[wrest]

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

не думаю что gcd() сильно быстрая.

Ну у меня на планшете так. Генерация миллиона пар случайных целых чисел в диапазоне до миллиона: 437 мс
? for(i=1,10^6,a=random(10^6);b=random(10^6))
? ##
*** last result computed in 439 ms.
То же, плюс вычисление их НОД: 729 мс
? for(i=1,10^6,a=random(10^6);b=random(10^6);c=gcd(a,b))
? ##
*** last result computed in 729 ms.
То же, но вместо вычисления НОД, вычисление суммы: 631 мс
? for(i=1,10^6,a=random(10^6);b=random(10^6);c=a+b)
? ##
*** last result computed in 631 ms.
То есть, вычисление НОД примерно в полтора раза медленнее вычисления суммы, для чисел в пределах миллиона.
Всё едят оверхеды (вызовы функций, проверки корректности аргументов и т.п.)

[gris]

Dmitriy40, спасибо. Но gcd для небольших чисел пригоден, а мне бы что-нибудь типа isvzaimnoprosty от массива трёх чисел из диапазона 1е60 - 1 е80. Хотя может быть это равносильно проверке делимости на последовательные простые, а вдруг более хитрые алгоритмы действуют

[Dmitriy40]

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

Dmitriy40, спасибо. Но gcd для небольших чисел пригоден, а мне бы что-нибудь типа isvzaimnoprosty от массива трёх чисел из диапазона 1е60 - 1 е80.

С чего бы? НОД вычисляется или делением (алгоритм Евклида), или сдвигами (лень смотреть как в PARI, по скорости похоже на сдвиги), и в том и в другом случае величина чисел может быть любой и к простоте (медленная проверка для больших чисел) они не привязаны:

? gcd([10^8000+13350,10^7000+4450,10^6000+22250])
%1 = 50

[gris]

А вот что даже совсем в начале:(

{M=1000; N=0; for(i=1,M,  for(j=i,M, for(k=j,M,
   if(gcd([i,j,k])==1,N++);
  ))); print(floor(100*N/(M*(M+1)*(M+2)/6)))}
83

И это у меня считалось 140 секунд. Увы... Я уж боюсь посчитать долю упорядоченных троек с не взаимно простыми числами не большими миллиона. Ясно, что она маленькая.

[wrest]

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

И это у меня считалось 140 секунд. Увы... Я уж боюсь посчитать долю упорядоченных троек с не взаимно простыми числами не большими миллиона. Ясно, что она маленькая.

Я полагаю, вам надо вычислять в два захода. Сначала факторизовать всё до миллиона, а потом работать с факторизациями. Думаю где-то на порядок будет быстрее (два порядка не обещаю).
Вам не надо хранить показатели, только основания. Можно сделать что-то вроде битовой подписи числа: маладший бит равен 1 если четное, следующий равен 1 если делится на три и т.п. У взаимнопростых побитовое И будет ноль. Возможно стоит ограничиться например первой сотней простых, уже на этом этапе почти все отсеется.

[Dmitriy40]

gris
Реализация идеи wrest с дополнительными оптимизациями:

? M=1000; N=0; for(i=1,M, for(j=i,M, for(k=j,M, if(gcd([i,j,k])==1, N++;); ))); print(floor(100*N/(M*(M+1)*(M+2)/6)))\\Ваш исходный вариант
83
time = 1min, 18,266 ms.

? M=1000; N=0; for(i=1,M, for(j=i,M, t=gcd(i,j); if(t==1, N+=M-j+1, for(k=j,M, if(gcd(t,k)==1, N++;);); ))); print(floor(100*N/(M*(M+1)*(M+2)/6)))\\Добавим ему оптимизации
83
time = 17,364 ms.

? M=2000; N=0; for(i=1,M, for(j=i,M, t=gcd(i,j); if(t==1, N+=M-j+1, for(k=j,M, if(gcd(t,k)==1, N++;);); ))); print(floor(100*N/(M*(M+1)*(M+2)/6)))\\Увеличим предел вдвое
83
time = 2min, 21,540 ms.

? M=1000; b=vector(M); for(x=2,M, t=0; foreach(factor(x)[,1],p, t=bitor(t,2^p)); b[x]=t;);\\Формирование факторизаций, единицы мс
? N=0; for(i=1,M, for(j=i,M, t=bitand(b[i],b[j]); for(k=j,M, if(bitand(t,b[k])==0, N++;); ))); print(floor(100*N/(M*(M+1)*(M+2)/6)))\\Используем
83
time = 45,959 ms.

? M=1000; b=vector(M); for(x=2,M, t=0; foreach(factor(x)[,1],p, t=bitor(t,2^p)); b[x]=t;);\\Формирование факторизаций, единицы мс
? N=0; for(i=1,M, for(j=i,M, t=bitand(b[i],b[j]); if(t==0, N+=M-j+1, for(k=j,M, if(bitand(t,b[k])==0, N++;);); ))); print(floor(100*N/(M*(M+1)*(M+2)/6)))\\С оптимизацией лишнего цикла
83
time = 14,572 ms.

? M=2000; b=vector(M); for(x=2,M, t=0; foreach(factor(x)[,1],p, t=bitor(t,2^p)); b[x]=t;);\\Формирование факторизаций, единицы мс
? N=0; for(i=1,M, for(j=i,M, t=bitand(b[i],b[j]); if(t==0, N+=M-j+1, for(k=j,M, if(bitand(t,b[k])==0, N++;);); ))); print(floor(100*N/(M*(M+1)*(M+2)/6)))\\Увеличим предел вдвое
83
time = 2min, 5,254 ms.

Выходит смысла в предрасчёте факторизаций в общем и нет, достаточно не делать лишних циклов.
Рост времени от M практически кубический, как и ожидалось (три вложенных цикла). Значит до миллиона не досчитаете, даже до сотни тысяч. И дело не только в медленности gcd или bitand, а в кубической зависимости.

-- 29.10.2022, 20:11 --

gris
С другой стороны, значения N для каждого M уже посчитаны в A015631. Там и программы на PARI приведены, намного быстрее:

? M=1000000; N=sum(k=1, M, sumdiv(k, d, moebius(k/d)*binomial(d+1, 2))); print("N=",N);print(floor(100*N/(M*(M+1)*(M+2)/6)));
N=138651542411519020
83
time = 5,836 ms.

Почему работает правильно — меня не спрашивайте.

-- 29.10.2022, 20:13 --

Т.е. в таких случаях полезно найти самому несколько первых (или близких к началу) значений самой простой функции из требуемых — и поискать их в OEIS.

[gris]

Спасибо! Мне неловко своей ерундой забивать тему, но всё-таки тут есть, чему поучиться именно в ПАРИ. Да, я подробно изучу. Собственно, я уже понял, что общие делители 2,3,5,7 и 11 уже дают немалый отсев троек, по-моему даже независимо от положения проверяемого интервала. Этакое скользящее решето Э:) То есть тройки можно отправлять на финальную проверку. А немного лишних замедлят ненадолго.

[Dmitriy40]

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

общие делители 2,3,5,7 и 11 уже дают немалый отсев троек, по-моему даже независимо от положения проверяемого интервала.

Разумеется независимо, по модулю же $\operatorname{lcm}(2,3,5,7,11)=2310$ начиная с любого выбранного числа все свойства делимости (и взаимной простоты соответственно) будут повторяться.

А насчёт троек, возможно 83% вообще предел на бесконечности, раз уж держится от 331 до сотни миллионов без исключений.

Кстати если нужны количества для M больше сотен миллионов (что займёт уже часы счёта), то можно заменить binomial() на выборку из предрасчитанной таблички, ведь оно вызывается много-много раз для достаточно малых аргументов и соответственно прилично тормозит счёт, примерно так: binomial(d+1,2) -> if(d<=#table, table[d], binomial(d+1,2)). Да и moebius() можно аналогично оптимизировать.

[wrest]

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

Значит до миллиона не досчитаете, даже до сотни тысяч. И дело не только в медленности gcd или bitand, а в кубической зависимости.

Да, надо признаться, я крупно погорячился. Просто пустой цикл (N=1000) идет у меня 47 секунд. Т.е. просто перебрать все упорядоченные тройки ничего с ними не делая будет 47 секунд, т.е. где-то 3,6 млн. троек в секунду. Проверка троек на взаимную простоту через gcd() добавляет ещё 40 секунд. Ну это же примерно следует и из моего предыдущего поста с таймингами. Таким образом, перебор всех троек не подходит вне зависимости от того, что с ними делается в переборе, да. Другое дело если чисел немного, но они большие, тогда, возможно именно gcd() будет тормозить. Это я тоже проверил, на тройках от $10^{60}$ до $10^{60}+10^3$, время увеличилось с 90 секунд до 140 секунд. Так что gcd() тут не главное что тормозит. Результат, для справки, тоже 83 процента.

-- 29.10.2022, 22:59 --

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

Почему работает правильно — меня не спрашивайте.

OMG! Гениальные люди... И правда, хоть бы написали почему.

[gris]

Вот чисто по PARI вопросы. Как наиболее просто считается модуль числа по основанию? Нет ли другого?

{print(lift(Mod(15,6)));}
3

Как считается If c || нескольких условий? Они по очереди проверяются до первого выполнения? Равносильно ли в смысле времени выполнения

if (a==1 || b==2 || c==3, N++) 
vs
if(a==1, if (b==2, if (c==3, N++)))

Столкнулся с непонятным замедлением. Переписал программку подсчёта доли выкидываемых троек с однородностью по 2,3,5 без функций со скользящими модулями.

{MA=1;MB=1000;N=0;
i2=0;i3=0;i5=0;
for (i=MA,MB,i2++;if(i2==2, i2=0);i3++;if(i3==3, i3=0);i5++;if(i5==5, i5=0);j2=0;j3=0; j5=0;
  for (j=1,i,j2++;if(j2==2, j2=0);j3++;if(j3==3, j3=0);j5++;if(j5==5, j5=0);k2=0;k3=0;k5=0;
    for (k=1,j,k2++;if(k2==2, k2=0);k3++;if(k3==3, k3=0);k5++;if(k5==5, k5=0);
       if(i2+j2+k2==0 || i3+j3+k3==0 || i5+j5+k5==0,N++);
)));print(floor(100*N/(MB*(MB+1)*(MB+2)/6)))}
16

Без предварительного расчёта массивов. Конечно, много лишних сложений и if-ов на простое равенство. Это тормозит?
Кстати при М=1000 из
167 167 000 троек
27 530 308 удалено по 2,3,5
28 188 452 по НОДу
Спасибо за 15%6

[wrest]

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

Как наиболее просто считается модуль числа по основанию?

Если "модуль числа по основанию" это остаток от [целочисленного] деления числа на основание, то так
? print(15%6)
3

-- 30.10.2022, 11:31 --

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

Они по очереди проверяются до первого выполнения? Равносильно ли в смысле времени выполнения

|| это ИЛИ, && это И :) так что написанное вами не эквивалентно логически
По идее, должно идти последовательно и не проверять второе если уже невыполнено первое, но я программистам не доверяю и всегда пишу вложенные if.