Перевод программ(ы) на язык функций (либо PARI)

[kthxbye]

Имеется последовательность A360288. В разделе программ оных есть две штуки.

Одна на Maple:

b:= proc(s, t) option remember; (m->

      `if`(m=0, x^(t/2), add(b(s minus {i}, t+

      `if`(i<m, 2^i, 0)), i=s)))(nops(s))

    end:

T:= n-> (p-> seq(coeff(p, x, i), i=0..degree(p)))(b({$1..n}, 0)):

seq(T(n), n=0..7);

Другая на языке Mathematica:

b[s_, t_] := b[s, t] = Function [m, If[m == 0, x^(t/2), Sum[b[s ~Complement~ {i}, t + If[i < m, 2^i, 0]], {i, s}]]][Length[s]];

T[n_] := CoefficientList[b[Range[n], 0], x];

Table[T[n], {n, 0, 7}]  // Flatten

Судя по комментарию к последней эти программы тождественны друг другу.

Может ли кто-нибудь перевести эти программы на язык функций, чтобы я потом попытался воспроизвести их на PARI? Вариант где программы переводятся сразу на язык PARI приветствуется еще больше.

Заранее благодарю.

[wrest]

Птичий язык и там и там.

[GAA]

Почему «птичий язык и там и там»? Синтаксис текста на Maple, на первый взгляд, понятный (Язык Mathematica я не знаю).

По тексту Maple

1. Процедура T возвращает последовательность коэффициентов многочлена, который формируется функцией b.
Если перевести с arrow на основной синтаксис процедур, то

>T := proc(n)
          (p-> seq(coeff(p, x, i), i=0..degree(p)))(b({$1..n}, 0)):
end proc;

Аргумент процедуры (n) используется для построения многочлена. Тело процедуры — вызов анонимной arrow-процедуры: в качестве единственного аргумента передаётся многочлен. Это можно переписать с arrow на основной синтаксис процедур

> pcoeff:= proc(p)
              seq(coeff(p, x, i), i=0..degree(p)):
end proc:
> T:= proc (n)
       pcoeff((b({$1..n}, 0))):
end proc:

Здесь seq — формирование последовательности значений коэффициентов при переменной x многочлена p: сначала коэффициент при нулевой степени x (т.е. свободный член), затем при первой и так до самой старшей). [Вызов coeff(p, x, i) — возвращает значение коэффициента многочлена p при переменной x в степени i]. Описание процедур coeff, degree можно найти в online справочнике: coeff, degree.
Может быть проще процедуру T воспринять в таком виде

> T:= proc (n)
  p := b({$1..n}, 0):
  seq(coeff(p, x, i), i=0..degree(p)):
end proc:

2. seq(T(n), n=0..7); — это демонстрация работы (процедур T и b). Возможно, удобней выводить в разные строки. Может быть проще понять в виде цикла:

 > for n from 1 to 4 do T(n) end do;
                   1
                  1, 1
               1, 3, 1, 1
           1, 7, 3, 7, 1, 3, 1, 1

3. Процедура b имеет два аргумента: s — множество чисел, t — число.
Тело процедуры b — вызов анонимной arrow-функции. Аргумент этой анонимной функции — число элементов множества s (nops(s) — возвращает число элементов s). Первоначально процедуре b передаётся множество, содержащее последовательность чисел от 1 до n (см. тело процедуры T). Таким образом m в теле b первоначально имеет значение n. В теле b стоит if-процедура. Если m=0 то возвращается x^(t/2) иначе выполняется суммирование вызова функции b с уменьшенным множеством и увеличенным значением t. Рано иди поздно количество элементов множества станет равно нулю и рекуррентные вызовы прекратятся.
___________________________________

kthxbye в сообщении #1585422 писал(а):

Может ли кто-нибудь перевести эти программы на язык функций

Я по ссылке (из Вашего сообщения в этой ветке) не ходил. PARI не знаю. Что значит перевести на язык функций?

[wrest]

GAA в сообщении #1585578 писал(а):

Почему «птичий язык и там и там»? Синтаксис текста на Maple, на первый взгляд, понятный

Ну там сразу же

b:= proc(s, t) option remember; (m->

`if`(m=0, x^(t/2), add(b(s minus {i}, t+

`if`(i<m, 2^i, 0)), i=s)))(nops(s))

end:

В процедуре b возникают переменные x и i которые никак не объявлены.
Как я понял ваше объяснение, x это не конкретная переменная а абстрактный символ (в том смысле что значение скажем х=0 или х=5 никого не интересует). Кажется, это назвается "литерал" или как-то так. А вот что там значит i ? На первый взгляд, первое обращение к i происходит до того, как этой переменной что-то присваивается. На второй взгляд, после присвоения (если это конечно присвоение) i=s переменная i становится того же типа что и переменная s то есть "множество". Тогда что значит например 2^i - это множество степеней двойки? Видимо есть какая-то тонкость, что просто i это количество элементов множества, а {i} само множество? Почему в s minus {i} в фигурных скобках только {i} но не s?

И кстати, а почему там `if` в апострофах?

[kthxbye]

GAA, благодарю за комментарий! Перечитал его несколько раз, но ввиду моих умственных способностей до конца не понял и буду еще дополнительно разбираться.

GAA в сообщении #1585578 писал(а):

Что значит перевести на язык функций?

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

Дополнительно присоединяюсь к вопросам wrest'а.

[GAA]

Предварительно по процедуре b

3A. В Maple `if` в апострофах — это (в данном случае) имя процедуры. if без апострофов — аналогичен управляющей конструкции в таких языках как Pascal или С. См. в oline справке по ifelse. `if` появилась раньше (в более ранних версиях) ifelse.
Выражение
`if`(conditional expression, true expression, false expression)
возвращает true expression, если conditional expression приняло значение true, иначе возвращает false expression.
В следующем примере `if` вернёт 1+ x+ x^2 (переменная x не инициализирована)

> `if`(false, 1, 1+x+x^2);
                        1+ x+ x^2

3B. "Строкиовые литералы" в Maple задаются в кавычках. Например,

> st := "мама":
> type(st, string);
                                                   true

3C. Переменная i входит в add. Это переменная суммирования. В простейшем виде add(f(i), i = 1..n) — это $\sum_{i=1}^n f(i)$. В данном случае переменной i присваивается диапазон (range). Диапазон 1..n — это 1, 2, 3, ..., n.
Переменной суммирования можно присвоить строку или множество. В следующем примере переменной суммирования присваивается множество чисел от 0 до 3 (т.е. {0,1,2, 3}) и add возвращает 1+x+x^2+x^3 (в этом месте переменная x не определена инициализирована).

 > add(x^i, i ={$0..3});
                                   1+x+x^2+x^3

См. oline справку add.
В процедуре b множество s — это множество чисел. 2^i это возведение числа 2 в степень, где в показателе одно из чисел множества s. Можно грубо считать add —цикл (for) с переменной цикла i. На каждой итерации переменная i принимает одно из значений указанного в «заголовке цикла».

3D. В процедуре b выражение s minus {i} — это разность множества s и множества {i}. В этом месте i — это один из элементов множества s (в данном случае одно из чисел множества s). Фигурные скобки, грубо говоря, конструктор множества. {i} — создать множество, содержащее один элемент i (одно из чисел s). s minus {i} — возвращает множество s без элемента i (без конкретного числа, которое было присвоена на итерации переменой суммирования i).

Исправление: определена зачеркнуто и заменено на инициализирована.

[GAA]

[В большинстве языков for, while, if не являются функциями. Функции, возвращающие количество элементов вектора/множества, на разных языках называются по-разному. Также по-разному называется функция для получения коэффициентов многочлена.]
Далее перевод на Maple без add, `if`, seq и анонимных функций, но с циклами и ветвлениями. Функция b из начального сообщения переименована в b1, однако во многих языках есть анонимные функции и аналоги ifelse, add и seq. Такое «расписывание» через циклы скорее затрудняет перевод, чем помогает. На мой взгляд, проще переводить с исходного текста (начальное сообщение). Просто по online справке посмотреть синтаксис функций.

> b1:= proc(s, t) option remember;
 local m, Sum, i, dt; 
 m := nops(s):
 if m = 0
  then Sum := x^(t/2)
  else 
   Sum:= 0;
   for i in s
    do 
     if i < m
      then dt := 2^i
      else dt := 0
     end if; 
     Sum:= Sum + b1(s minus {i}, t + dt):
    end do:
 end if:
 Sum:
end proc:
Тест
> n := 1: S := {$1..n};
          S := {1}
> b1(S, 0);
            1
> n := 2: S := {$1..n};
        S := {1, 2} 
> b1(S, 0);
           1 + x
> n := 3: S := {$1..n};
        S := {1, 2, 3}
> b1(S, 0);
      3*x+x^3+x^2+1
> n := 4: S := {$1..n};
        S := {1, 2, 3, 4}
> b1(S, 0);
      7*x^3+7*x+3*x^2+3*x^5+x^7+x^6+x^4+1

> T:= proc (n)
     local p, i, L:
     p := b1({$1..n}, 0):
     L:= [];
     for i from 0 to degree(p) do L := [op(L),  coeff(p, x, i)] end do;
     op(L):
end proc:
> for i from 1 to 4 do T(i) end do;
                           1
                         1, 1
                     1, 3, 1, 1
               1, 7, 3, 7, 1, 3, 1, 1

После перезагрузки Maple 15 Classic Worksgeet всё вместе

> b1:= proc(s, t) option remember;
   local m, Sum, i, dt; 
   m := nops(s):
   if m = 0
    then Sum := x^(t/2)
    else 
     Sum:= 0;
     for i in s
      do 
       if i < m
        then dt := 2^i
        else dt := 0
       end if; 
       Sum:= Sum + b1(s minus {i}, t + dt):
      end do:
   end if:
   Sum:
  end:
> T:= proc (n)
    local p, i, L:
    p := b1({$1..n}, 0):
    L:= [];
    for i from 0 to degree(p) do L := [op(L),  coeff(p, x, i)] end do;
    op(L):
  end proc:
> for i from 1 to 5 do T(i) end do;
                          1
                         1, 1
                     1, 3, 1, 1
               1, 7, 3, 7, 1, 3, 1, 1
   1, 15, 7, 31, 3, 17, 7, 15, 1, 7, 3, 7, 1, 3, 1, 1

[wrest]

GAA
Уже яснее (мне по крайней мере). А вот это
for i in s
выполняется в случайном порядке? Ну и, если в s остались скажем {4,5,7} то и i будет итерироваться только принимая значения 4,5,7
При этом m:=nops(s) будет равно 3 и соответсвенно
if i < m
then dt := 2^i
else dt := 0
end if;
все три раза пойдёт по ветке dt := 0

Я просто начал эту рекурсию писать, но она у меня бесконечная получается пока что.

-- 16.03.2023, 19:50 --

kthxbye
Ну вот кажется функция b написалась...

? b(s,t)=if(#s==0,return(x^(t/2)));sum1=0;foreach(s,i,sum1=sum1+b(setminus(s,[i]),t+if(i<#s,2^i,0)));return(sum1);
? b([1,2,3,4],0)
%9 = x^7 + x^6 + 3*x^5 + x^4 + 7*x^3 + 3*x^2 + 7*x + 1
?

Тут не используется мемоизация, поэтому для больших множеств будет работать, возможно, долго.

P.S. Да, для множества из 10 считает 30 секунд...

[GAA]

wrest, в общем случае порядок элементов в множестве не фиксирован. Если нужен порядок, то вместо множества нужно использовать список (list). Но в данном случае выполняется в Maple 15 по порядку.
Более точно смотри ниже (я добавил в b1 вывод: print (S = s) )

 > b1:= proc(s, t) option remember;
 local m, Sum, i, dt; 
 print(S = s):
 m := nops(s):
 if m = 0
  then Sum := x^(t/2)
  else 
   Sum:= 0;
   for i in s
    do 
     if i < m
      then dt := 2^i
      else dt := 0
     end if; 
     Sum:= Sum + b1(s minus {i}, t + dt):
    end do:
 end if:
 Sum:
end:

 > for i from 1 to 3 do T(i) end do;
   S = {1}
   S = {}
     1
   S = {1, 2}
   S = {2}
   S = {}
     1, 1
   S = {1,2,3}
   S = {2,3}
   S = {3}
   S = {1, 3}
   S = {3}
   S = {}
   S = {1}
   S = {}
     1, 3, 1, 1

В Maple 15 зацикливания для n = 7 нет

> T(7);
   1, 63, 31, 391, 15, 245, 115, 675, 7, 145, 69, 445, 31, 261, 115, 391, 3, 77, 37, 261, 17, 155, 69, 245, 7, 81, 37, 145, 15, 77, 31, 63, 1, 31, 15, 115, 7, 69, 31, 115, 3, 37, 17, 69, 7, 37, 15, 31, 1, 15, 7, 31, 3, 17, 7, 15, 1, 7, 3, 7, 1, 3, 1, 1

Аналогично для n =9, n = 11, n = 13. Дальше не проверял.

[wrest]

GAA в сообщении #1585677 писал(а):

В Maple 15 зацикливания для n = 7 нет

Да, это я там неправильно понял.
В pari есть аналог цикла "для каждого из" for i in s do f end do, синтаксис похожий, foreach(s,i,f).

[kthxbye]

GAA, еще раз благодарю за столь подробное объяснение! Извините, что назвал функциями то, что функциями, как оказалось, отнюдь не является.

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

[wrest]

kthxbye в сообщении #1585689 писал(а):

Вы не могли бы написать мемоизацию и просто оценить скорость работы программы уже с нею?

Мог бы. Множество из 10 обсчитывается 0,3с вместо 30с без мемоизации.

-- 16.03.2023, 21:34 --

kthxbye в сообщении #1585689 писал(а):

с которой я так до сих пор и не разобрался

Самое время!

b1(s,t) =
  { local(M = Map());
    my(memo_b = (S,T) -> my(res,sum1=0);
       if(#S==0, return(x^(T/2)));
       if(mapisdefined(M,[S,T],&res),  return(res));
       foreach(S,i,sum1=sum1+self()(setminus(S,[i]),T+if(i<#S,2^i,0)));
       res=sum1;
       mapput(M,[S,T],res);
       return(res));
    memo_b(s,t);
  }

Функция выше b1(s,t) возвращает многочлен. Что с ним делать дальше разберётесь?

-- 16.03.2023, 21:47 --

GAA
Спасибо большое! Ваши комментарии очень помогли разобраться в птичьем языке Maple.
Думаю что приведенная выше функция на pari/gp выглядит не менее птичей, с учетом того что пришлось применить костыли в виду отсутствия встроенной мемоизации (то что определяется option remember в Maple).
Так-то концепции схожие и там и там есть (даже arrow), но тонкости синтаксиса не позволяют понять написанное сходу.

[wrest]

Немного статистики по работе функции.
Считался результат для входных параметров: s это множество из $n=10$ чисел (от 1 до 10),t=0
В варианте честной рекурсии без мемоизации, количество вызовов составило 9846101 (т.е. $e\cdot 10!$), из них количество вызовов с пустым множеством равно 3628800 (т.е. $10!$)
В варианте с мемоизацией, общее количество вызовов составило 267385, количество попаданий в память 187240, количество вызовов с пустым множеством 7168 (т.е. $7\cdot 2^{10}$)
То есть, мемоизация в этом случае сократила количество вызовов в 37 раз, и при этом 70% вызовов - попадания в память.

Кстати как ни удивительно, но для множества из 12 чисел все количества (вызовов функции с мемоизацией) те же, что и для 10 чисел.

[wrest]

kthxbye
Функция ниже T(n)вызывает функцию b1(s,t) и работает аналогично приведенной в OEIS: возвращает вектор из $2^{n-1}$ значений n-й строки треугольника.

? T(n)=vector(2^(n-1),i,polcoef(b1(Set(vector(n,i,i)),0),i-1));
? print(T(5))
[1, 15, 7, 31, 3, 17, 7, 15, 1, 7, 3, 7, 1, 3, 1, 1]

Функция T_fl(n) вызывает функцию T(n) и возвращает собственно последовательность (первые $2^n$ членов), составленную из n первых строк треугольника. Тут обращу ваше внимание, что в OEIS нумеруют с нуля, так что мной вначале мной добавлена единица как нулевой член последовательности, чтобы вывод соответствовал OEIS.

? T_fl(n)=my(v=T(n),v1=[1]);for(i=1,n,v1=concat(v1,v[(2^(n-1)-2^(i-1)+1)..2^(n-1)]));return(v1);
? print(T_fl(5))
[1, 1, 1, 1, 1, 3, 1, 1, 1, 7, 3, 7, 1, 3, 1, 1, 1, 15, 7, 31, 3, 17, 7, 15, 1, 7, 3, 7, 1, 3, 1, 1]
?

В функции T_fl(n) коэффициенты реально вычисляются только для n-й строки треугольника, так как там видно, что верхняя т.е. $(n-1)$-я строка в треугольнике это правая половина строки под ней, т.е. $n$-й, но доказательство этого за ваш счёт.

[kthxbye]

wrest, еще раз огромное вам спасибо! Статистика интересная, замечание про 10 и 12 действительно очень удивительное. Вы уверены что корректно вычислили значения, а не посчитали дважды для 10 или для 12? Я вам доверяю, но все-таки хочется уточнить этот нюанс.

Синтаксис функций T(n) и T_fl(n) абсолютно понятен.

Скорость работы программы - вот все что мне нужно было узнать. Дело в том, что для тех же значений я нашел очень интересный на мой взгляд алгоритм, быстрее которого, вероятно, ничего нет. Выглядит он так:

a(n)=my(n=(n+1)/2^valuation(n+1, 2), A=n, B, C, v=[], v1); while(A>0, B=valuation(A, 2); v=concat(v, B+1); A\=2^(B+1)); v=Vecrev(v); A=#v; if(n==1, 1, v1=vector(A, i, A-i+1); for(i=1, A-1, B=A-i; for(j=1, B, C=B-j+2; v1[j]=v1[j]*C^v[B] - v1[j+1]*(C-1)^v[B])); v1[1])

Нумерация начинается с нуля. Алгоритм позволяет с высокой скоростью получать подпоследовательности A110501 через $a(\frac{4^n-1}{3})$.