text:
read ("metaMPFR_straightforwardAlgo.mpl"):

f := AiryAi(x):
deq := holexprtodiffeq(f, y(x)):
rec := diffeqtorec(deq, y(x), a(n)):
name_of_function := op(0,f):
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

f := erf(x):
deq := holexprtodiffeq(f, y(x)):
rec := diffeqtorec(deq, y(x), a(n)):
name_of_function := op(0,f):
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = -(6*n+1)*(6*n+2)*(6*n+3)*(6*n+4)*(6*n+5)*(6*n+6)*a(n)/( (n+1)^3*(3*n+1)*(3*n+2)*(3*n+3)*12288000 ), a(0)=1 }:
name_of_function := alpha:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1 }:
name_of_function := test0a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi }:
name_of_function := test1a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1 }:
name_of_function := test2a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi }:
name_of_function := test3a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=2 }:
name_of_function := test4a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi, a(1)=0}:
name_of_function := test5a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1, a(1)=Pi }:
name_of_function := test6a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(1)=0}:
name_of_function := test7a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=0, a(2)=0 }:
name_of_function := test8a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=0, a(1)=Pi, a(2)=2 }:
name_of_function := test9a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+3) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=0, a(1)=0, a(2)=1 }:
name_of_function := test10a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+7) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(4)=1, a(6)=2  }:
name_of_function := test11a:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), CONSTANT_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1 }:
name_of_function := test0b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi }:
name_of_function := test1b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1 }:
name_of_function := test2b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi }:
name_of_function := test3b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=2 }:
name_of_function := test4b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi, a(1)=0}:
name_of_function := test5b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1, a(1)=Pi }:
name_of_function := test6b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(1)=0}:
name_of_function := test7b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=0, a(2)=0 }:
name_of_function := test8b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=0, a(1)=Pi, a(2)=2 }:
name_of_function := test9b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+3) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=0, a(1)=0, a(2)=1 }:
name_of_function := test10b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+7) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(4)=1, a(6)=2  }:
name_of_function := test11b:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES, name_of_function, name_of_file, f):

rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1 }:
name_of_function := test0c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+1) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi }:
name_of_function := test1c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1 }:
name_of_function := test2c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+1) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi }:
name_of_function := test3c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=2 }:
name_of_function := test4c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+2) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=Pi, a(1)=0}:
name_of_function := test5c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=1, a(1)=Pi }:
name_of_function := test6c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+2) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(1)=0}:
name_of_function := test7c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=1, a(1)=0, a(2)=0 }:
name_of_function := test8c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+3) = -(6*n+1)*a(n)/( (n+1)^3 ), a(0)=0, a(1)=Pi, a(2)=2 }:
name_of_function := test9c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+3) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=0, a(1)=0, a(2)=1 }:
name_of_function := test10c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):

rec := { a(n+7) = (1/Pi)*(3*n+1)*a(n)/( (n+2) ), a(0)=Pi, a(4)=1, a(6)=2  }:
name_of_function := test11c:
name_of_file := sprintf("%a.c", name_of_function):
printf("\n\n\n/************************** Implémentation de %s ******************************/\n", name_of_file):
generateStraightforwardAlgo(rec, a(n), FUNCTION_SERIES_RATIONAL, name_of_function, name_of_file, f):


text:
var t1 : int := 0;
var t2 : int := 1;
var t3 : int := t1 + t2;
var nTimes : int := 0;
print "nth term? ";
read nTimes;
var x : int;
for x in 2..nTimes-1 do
  t1:=t2;
  t2:=t3;
  t3:=t1+t2;
end for;
print t3;
print "\n";

ext: mpl

max_stars_repo_name: nicohi/mini-pl-interpreter

max_issues_repo_name: nicohi/mini-pl-interpreter

max_forks_repo_name: nicohi/mini-pl-interpreter


text:
moore_net_complex := proc(n::posint)
 local T,m,i;

 m := 3*n;

 T["vertices"] := {seq(i,i=1..2*m+1)};
 T["faces"] := {
  seq({1,2+i,2+modp(i+1,m)},i=0..m-1),
  seq({2+i,2+modp(i+1,m),m+2+modp(i+1,m)},i=0..m-1),
  seq({2+i,m+2+i,m+2+modp(i+1,m)},i=0..m-1)
 };

 T["max_simplices"] := T["faces"];

 T["embedding_dim"] := 2;
 T["embedding"] := table([
  1 = [0,0],
  seq(i = [2*cos(2*Pi*(i-2)/m),2*sin(2*Pi*(i-2)/m)],i=2..m+1),
  seq(i = [3*cos(2*Pi*(i-m-2)/m),3*sin(2*Pi*(i-m-2)/m)],i=m+2..2*m+1)
 ]);

 T["plot"] :=
   `plot/simplicial_complex`(T["vertices"])(T["faces"],2,T["embedding"]);

 return eval(T);
end:

moore_complex := proc(n::posint)
 local T0,T,m,i;

 m := 3*n;

 T["vertices"] := {seq(i,i=1..m+4)};
 T["faces"] := {
  seq({1,2+i,2+modp(i+1,m)},i=0..m-1),
  seq({2+i,2+modp(i+1,m),m+2+modp(i+1,3)},i=0..m-1),
  seq({2+i,m+2+modp(i,3),m+2+modp(i+1,3)},i=0..m-1)
 };
 T["max_simplices"] := T["faces"];

 return eval(T);
end:

punctured_moore_complex := proc(n::posint)
 local T,E,F,m,a,R0,R1,tor,tor_arc,i,v;

 m := 3*n;

 T["vertices"] := {seq(i,i=2..m+4)};
 F := {
  seq({2+i,2+modp(i+1,m),m+2+modp(i+1,3)},i=0..m-1),
  seq({2+i,m+2+modp(i,3),m+2+modp(i+1,3)},i=0..m-1)
 };
 T["faces"] := F;
 T["max_simplices"] := T["faces"];
 E := map(op,map(f -> {{f[1],f[2]},{f[1],f[3]},{f[2],f[3]}},F));
 T["edges"] := E;

 R0 := 3;
 R1 := 1;

 tor := u ->
  (R0+R1*u[3]*cos(2*Pi*u[2])) *~ [cos(2*Pi*u[1]),sin(2*Pi*u[1]),0] +~
   [0,0,R1*u[3]*sin(2*Pi*u[2])];

 a := table([
  seq(  2+i = evalf([i/3,i/m,1]),i=0..m-1),
  seq(m+2+i = evalf([i/3,i/m,0]),i=0..2)
 ]);

 T["tor_embedding"] := eval(a);

 T["embedding_dim"] := 2;
 T["embedding"] := table([
  seq(v = tor(a[v]),v in T["vertices"])
 ]);

 T["plot"] :=
   `plot/simplicial_complex`(T["vertices"])(T["faces"],2,T["embedding"]);

 tor_arc := proc(u,v)
  local w;
  w := v -~ u;
  w := [w[1] - round(w[1]),w[2] - round(w[2]),w[3]];
  spacecurve(tor(u +~ t *~ w),t=0..1);
 end:

 T["curved_plot"] :=
  display(map(e -> tor_arc(a[e[1]],a[e[2]]),E),
   axes=none,scaling=constrained);

 return eval(T);
end:


text:
# 尽可能的化简解的约束

$ifndef _CON_REFINE_
$define _CON_REFINE_

$include "Condition.mpl"
$include "../seg/Seg.mpl"

# 总的简化函数
conRefine:=proc(r::RepSol)
    signRefine(r);
    epRefine(r);
    uniqueAndSort(r);# 简化之后重新对条件进行排序
    return;
end proc:

# 对于代表元v[k]删去a[k]>0.
# 对于代表元-v[k]删去a[k]<0.
# 这里只删去了acon的条件，对于osol并没有进行修改
signRefine:=proc(r::RepSol)
    local vs,k,rc;
    vs:=indets(r:-rep,name);
    if numelems(vs)<>1 then
        return;
    end if;
    k:=op([1,1],vs);
    rc:=eval({a[k]>0,a[k]<0});
    r:-acon:=map(x->x minus rc,r:-acon);
    return;
end proc:

# 尽可能删去和epsilon有关的约束
epRefine:=proc(r::RepSol)
    r:-acon:=epDeal~(r:-acon);
end proc:

# 处理每个可能情况的约束条件
epDeal:=proc(s::set)
    local ca,ce,r;
    ce,ca:=selectremove(has,s,epsilon);
    ce:=ceDeal~(ce);
    r:=ce union ca;
    r:=classifySolve(r);
    return r;
end proc:

# 单个关于epsilon的条件处理
ceDeal:=proc(e)
    local rs;
    # 有的时候solve会抽，对于这种情况不做处理
    try
        rs:=conSolve(e);
    catch :
        return e;
    end try;
    if numelems(rs)=0 then
        # 无解则原样返回
        return e;
    elif numelems(rs)=1 then
        # 只有一个解，则不论有多少条件都可以
        return rs[][];
    else
       return multiConRefine(rs);
    end if;
end proc:

multiConRefine:=proc(rs)
    local s;
    # 单变量约束的并集，只处理无约束和，x<>c的约束
    if numelems(indets(rs,name))=1 and
        andmap(x->evalb(numelems(x)=1),rs) then
        s:=`union`(seq(Seg(x),x in rs));
        if s=real then
            return NULL;
        end if;
        s:=&C s;
        if type(s:-bound,numeric) then
            return indets(rs,name)[]<>s:-bound;
        else
            return e;
        end if;
    # 删去一个共同条件之后能够构成单变量约束的并集
    elif numelems(`intersect`(rs[]))=1 then
        s:=`intersect`(rs[]);
        return s[],multiConRefine(map(x->x minus s,rs));
    else
        return e;
    end if;
end proc:

# 求解不等式，
# 在解中删除自由变量
# 删除和epsilon有关的约束
conSolve:=proc(c)
    local r;
    r:=[RealDomain:-solve(c)];
    r:=map(x->remove(t->evalb(lhs(t)=rhs(t)),x),r);
    r:=map(x->remove(has,x,epsilon),r);
    r:=[{r[]}[]];# 去重
    return r;
end proc:

$endif


text:
TITLE
    Invest3Scen;

INDEX
    invtype := (Stock, Bond) -> (S,B);

STOCHASTIC

STAGES
    period := 1..3;

SCENARIOS
    sc := 1..4;

TREE
    BINARY;

PROBABILITIES
    prob[sc] := ALLEQUAL;

RANDOM DATA
    Returns[invtype, period, sc] :=
      [Stock, 1, 1, 1.25,
       Stock, 1, 2, 1.25,
       Stock, 1, 3, 1.06,
       Stock, 1, 4, 1.06,
       Stock, 2, 1, 1.25,
       Stock, 2, 2, 1.06,
       Stock, 2, 3, 1.25,
       Stock, 2, 4, 1.06,
       Bond, 1, 1, 1.14,
       Bond, 1, 2, 1.14,
       Bond, 1, 3, 1.12,
       Bond, 1, 4, 1.12,
       Bond, 2, 1, 1.14,
       Bond, 2, 2, 1.12,
       Bond, 2, 3, 1.14,
       Bond, 2, 4, 1.12];

DATA
    InitialWealth = 55;
    Goal := 80;
    SurplusReward := 1;
    ShortPenalty  := 4;

VARIABLES
    Invest[period<LAST(period),invtype] -> x;

LASTSTAGE VARIABLES
    Shortage -> Sh;
    Excess -> Ex;

MODEL

   MAX obj = SurplusReward*Excess - ShortPenalty*Shortage;

SUBJECT TO

   InitInvest[period=1] -> InitInv:
       InitialWealth
     =
       SUM(invtype: Invest);

   InvestBal[period=2..LAST(period)-1] -> InvBal:
       SUM(invtype: Invest[period-1] * Returns[period-1])
     =
       SUM(invtype: Invest);

   EndInvest[period=LAST(period)] -> EndInv:
       SUM(invtype: Invest[period-1]* Returns[period-1])  + Shortage - Excess
     =
       Goal;

END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


text:
bar := proc()
 local n,m,i,j,u,a,v,ac,av,x,y;
 n := nargs;
 for i from 1 to n do
  u := seq(args[j],j=1..i-1);
  a := args[i];
  v := seq(args[j],j=i+1..n);
  if type(a,`+`) then
   return map(t -> bar(u,t,v),a);
  elif type(a,`*`) then
   ac,av := selectremove(type,a,rational);
   if ac <> 1 then
    return ac * bar(u,av,v);
   fi;
  elif type(a,specfunc(bar)) then
   return bar(u,op(a),v);
  elif type(a,specfunc(bracket)) then
   m := nops(a);
   if m = 0 then
    return FAIL;
   elif m = 1 then
    return bar(u,op(a),v);
   else
    x := op(1,a);
    y := bracket(op(2..-1,a));
    return bar(u,x,y,v) - bar(u,y,x,v);
   fi;
  fi;
 od;
 return 'procname'(args);
end:

`mult/free_algebra` := () -> bar(args);

`commutator/free_algebra` := () -> bar(bracket(args));

`basis/free_algebra` := (A::set) -> proc(d::nonnegint)
 local i,B,a,b;
 B := [bar()];
 for i from 1 to d do
  B := [seq(seq(bar(op(b),a),a in A),b in B)];
 od:
 return B;
end:

`vec/free_algebra` := (A::set) -> (d::nonnegint) -> proc(u)
 local n,ub,uc,T,B,i;
 global `vec_table/free_algebra`;
 n := nops(A);
 if u = 0 then
  return [0$(n^d)];
 elif type(u,`+`) then
  return map(`vec/free_algebra`(A)(d),u);
 else
  if type(u,specfunc(bar)) then
   ub := u;
   uc := 1;
  elif type(u,`*`) then
   ub,uc := selectremove(type,u,specfunc(bar));
   ub := map(op,bar(ub));
  else
   return FAIL;
  fi;
  if not(type(`vec_table/free_algebra`[A,d],table)) then
   T := table():
   B := `basis/free_algebra`(A)(d);
   for i from 1 to n^d do
    T[B[i]] := [0$(i-1),1,0$(n!-i)];
   od:
   `vec_table/free_algebra`[A,d] := eval(T);
  fi;
  return uc *~ `vec_table/free_algebra`[A,d][ub];
 fi;
end:


text:
# Proof by PDF that X <~ rayleigh(1/sqrt(2*lambda)) == X^2 <~ exponential(lambda)

PDF_rayleigh := (a, x) -> 2*x*exp(-x^2/a)/a;

PDF_transform := (x) -> PDF_rayleigh(1,sqrt(x))*diff(sqrt(x),x);
simplify(PDF_transform(x));

PDF_exponential := (l, x) -> l*exp(-l*x);
simplify(PDF_exponential(1, x));

evalb(PDF_exponential(1, x) = PDF_transform(x));

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
#@ Not autoload

# (This has been partially converted to autoloadable form, but
# needs more work.)

# This file contains formulae for elliptic curves taken from
# "Courbes Elliptiques: Formulaire (d'apr\`es J. Tate)"
# by P. Deligne (pages 53-73 of Modular Functions of One
# Variable IV, LNM 476, edited by B.J. Birch and W. Kuyk, 1975)

`is_element/weierstrass_vectors` := proc(a,t := constant)
 return type(a,[t $ 5]);
end:

`weierstrass_cubic` := (a) -> proc(xy)
 local a1,a2,a3,a4,a6,x,y;
 a1,a2,a3,a4,a6 := op(a);
 x,y := op(xy);
 return y^2 + a1 * x * y + a3 * y - x^3 - a2 * x^2 - a4 * x - a6;
end:

`homogeneous_cubic` := (a) -> proc(xyz)
 local a1,a2,a3,a4,a6,x,y,z;
 a1,a2,a3,a4,a6 := op(a);
 x,y,z := op(xyz);
 return y^2*z + a1*x*y*z + a3*y*z^2 - x^3 - a2*x^2*z - a4*x*z^2 - a6*z^3;
end:

`a_to_b/weierstrass_vectors` := proc(a)
 local a1,a2,a3,a4,a6,b2,b4,b6,b8;

 a1,a2,a3,a4,a6 := op(a);

 b2 := a1^2 + 4 * a2;
 b4 := a1 * a3 + 2 * a4;
 b6 := a3^2 + 4 * a6;
 b8 := a1^2 * a6 - a1 * a3 * a4 + 4 * a2 * a6 + a2 * a3^2 - a4^2;

 return [b2,b4,b6,b8];
end:

`b_to_c/weierstrass_vectors` := proc(b)
 local b2,b4,b6,b8,c4,c6,delta;

 b2,b4,b6,b8 := op(b);

 c4 := b2^2 - 24 * b4;
 c6 := - b2^3 + 36 * b2 * b4 - 216 * b6;

 delta := -b2^2 * b8 - 8 * b4^3 - 27 * b6^2 + 9 * b2 * b4 * b6;

 return [c4,c6,delta];
end:

`c_to_j/weierstrass_vectors` := proc(c)
 local c4,c6,delta;

 c4,c6,delta := op(c);

 return c4^3/delta;
end:

`b_to_j/weierstrass_vectors` := proc(b)
 `c_to_j/weierstrass_vectors`(`b_to_c/weierstrass_vectors`(b));
end:

`a_to_c/weierstrass_vectors` := proc(a)
 `b_to_c/weierstrass_vectors`(`a_to_b/weierstrass_vectors`(a));
end:

`a_to_j/weierstrass_vectors` := proc(a)
 `b_to_j/weierstrass_vectors`(`a_to_b/weierstrass_vectors`(a));
end:

`b_rel/weierstrass_vectors` := proc(b)
 local b2,b4,b6,b8;

 b2,b4,b6,b8 := op(b);

 return 4 * b8 - b2 * b6 + b4^2;
end:

`c_rel/weierstrass_vectors` := proc(c)
 local c4,c6,delta;

 c4,c6,delta := op(c);

 return 1728 * delta - c4^3 + c6^2;
end:

`invariant_differential0/weierstrass_vectors` :=
 (a) -> (y * dx - x * dy)/diff(homogeneous_cubic(a)([x,y,z],z);

`invariant_differential1/weierstrass_vectors` :=
 (a) -> (z * dy - y * dz)/diff(homogeneous_cubic(a)([x,y,z],x);

`invariant_differential2/weierstrass_vectors` :=
 (a) -> (x * dz - z * dx)/diff(homogeneous_cubic(a)([x,y,z],y);

######################################################################

`id/weierstrass_group` := [1,0,0,0]:

`o/weierstrass_group` := proc(urst1,urst0)
 local u0,u1,u2,r0,r1,r2,s0,s1,s2,t0,t1,t2;

 u1,r1,s1,t1 := op(urst1);
 u0,r0,s0,t0 := op(urst0);
 u2 := u0 * u1;
 s2 := s1 + s0 * u1;
 r2 := r1 + r0 * u1^2;
 t2 := t1 + t0 * u1^3 + r0 * s1 * u1^2;
 return [u2,r2,s2,t2];
end:

`inv/weierstrass_group` := proc(urst)
 local u,r,s,t;
 u,r,s,t := op(urst);
 return [1/u,-r/u^2,-s/u,(r*s-t)/u^3];
end:

`act/vector/weierstrass_group` := (urst) -> proc(a)
 local u,r,s,t,a1,a2,a3,a4,a6;
 u,r,s,t := op(urst);
 a1,a2,a3,a4,a6 := op(a);
 return [
  a1 * u - 2 * s,
  a2 * u^2 + a1 * s * u - 3 * r - s^2,
  a3 * u^3 - a1 * r * u + 2 * r * s - 2 * t,
  a4 * u^4 + a3 * s * u^3 - 2 * a2 * r * u^2 +
    a1 * (t - 2 * r * s) * u + 3 * r^2 + 2 * r * s^2 - 2 * s * t,
  a6 * u^6 - a4 * r * u^4 +
    a3 * (t - r * s) * u^3 + a2 * r^2 * u^2 +
    a1 * (r^2 * s - r * t) * u + 2 * r * s * t - t^2 - r^2 * s^2 - r^3
 ];
end:

`act/curve/weierstrass_group` := (urst) -> proc(xyz)
 local u,r,s,t,x,y,z;
 u,r,s,t := op(urst);
 x,y,z := op(xyz);
 return [u^2 * x + r * z,u^3 * y + s * u^2 * x + t * z,z];
end:

#### UNFINISHED

# Three expressions for the invariant differential
alpha0 := (y * dx - x * dy)/diff(HomogeneousCubic,z);
alpha1 := (z * dy - y * dz)/diff(HomogeneousCubic,x);
alpha2 := (x * dz - z * dx)/diff(HomogeneousCubic,y);

eta := y + (a1 * x + a3)/2;
xi  := x + b2 / 12;

WeierstrassCubic2 := eta^2 - x^3 - b2/4 * x^2 - b4/2 * x + b6/4;
WeierstrassCubic6 := eta^2 - xi^3 + c4/48 * xi + c6/864;

Pu := xi;
Pdashu := 2 * eta;
g2 := c4/12;
g3 := c6/216;
du := omega;
dxi := 2 * eta * omega;

# alpha[t_Plus]    := Expand[alpha /@ t]
# alpha[t_Times]   := Expand[alpha /@ t]
# alpha[t_^n_]     := Expand[alpha[t]^n]
# alpha[n_Integer] := n

checkalpha := {
alpha[b2] - (b2 * u^2 - 12 * r),
alpha[b4] - (b4 * u^4 - b2 * r * u^2 + 6 * r^2),
alpha[b6] - (b6 * u^6 - 2 * b4 * r * u^4 + b2 * r^2 * u^2 - 4 * r^3),
alpha[b8] - (b8 * u^8 - 3 * b6 * r * u^6 + 3 * b4 * r^2 * u^4 - b2 * r^3 * u^2 + 3 * r^4),
alpha[c4] - c4 * u^4,
alpha[c6] - c6 * u^6,
alpha[Delta] - Delta * u^12
};

CurveRule := proc(F)
  [a1 =   coeff(coeff(F,x,1),y,1),
   a2 = - coeff(coeff(F,x,2),y,0),
   a3 =   coeff(coeff(F,x,0),y,1),
   a4 = - coeff(coeff(F,x,1),y,0),
   a6 = - coeff(coeff(F,x,0),y,0)
  ]
end:

GenericCubic := y^2 + x * y - x^3 + 36/(jj - 1728) * x + 1/(jj - 1728);

mChord := proc(P,Q)
 (Q[2] - P[2])/(Q[1] - P[1]);
end:

mTangent := proc(P)
 (3 * P[1]^2 + 2 * a2 * P[1] + a4 - a1 * P[2]) / (2 * P[2] + a1 * P[1] + a3);
end:

bChord := proc(P,Q)
 (P[2] * Q[1] - Q[2] * P[1])/(Q[1] - P[1]);
end:

bTangent := proc(P)
 (- P[1]^3 + a4 * P[1] + 2 * a6 - a3 * P[2]) / (2 * P[2] + a1 * P[1] + a3);
end:

EllipticSumChord := proc(P,Q)
 local m,b,x3;
  m := mChord(P,Q);
  b := bChord(P,Q);
  x3 := m^2 + a1 * m - a2 - P[1] - Q[1];
  [x3, -(m + a1) * x3 - b - a3];
end:

EllipticSumTangent := proc(P,Q)
 local m,b,x3;
  m := mTangent(P);
  b := bTangent(P);
  x3 := m^2 + a1 * m - a2 - P[1] - Q[1];
  [x3, -(m + a1) * x3 - b - a3]
end:

Negate := proc(P)
 if nops(P) = 2 then
  return([P[1],-P[2]-a1*P[1]-a3]);
 else
  return([P[1],-P[2]-a1*P[1]-a3*P[3],P[3]]);
 fi;
end:

Double := proc(x)
 (x^4 - b4 * x^2 - 2 * b6 * x - b8)/(4 * x^3 + b2 * x^2 + 2 * b4 * x + b6);
end:

# xTate = sum[q^n u / (1 - q^n u)^2 , {n,-Infinity,Infinity}] -
#         sum[n q^n / (1 - q^n) , {n,1,Infinity}]

# yTate = sum[q^(2 n) u^2 / (1 - q^n u)^3 , {n,-Infinity,Infinity}] +
#         sum[n q^n / (1 - q^n) , {n,1,Infinity}]

# TateCubic = y^2 + x y - x^3 - aTate4 x - aTate6

# aTate4 = - 5 sum[ DivisorSigma[3,n] q^n , {n,1,Infinity}]
# aTate6 = - sum[ (5 DivisorSigma[3,n] + 7 DivisorSigma[5,n])/12 q^n,
#                  {n,1,Infinity}]

# cTate4 = E4
# cTate6 = -E6

# E4 = 1 + 240 sum[ DivisorSigma[3,n] q^n , {n,1,Infinity}]
# E6 = 1 - 504 sum[ DivisorSigma[5,n] q^n , {n,1,Infinity}]
# DeltaTate = q product[ (1 - q^n)^24 , {n,1,Infinity}]


text:
input := FileTools:-Text:-ReadFile("AoC-2021-16-input.txt" ):

binstr2dec := proc(b) convert(parse(b),decimal,2) end proc:
hex2bits := proc(input) uses StringTools; local tmp := Explode(input); Join(subs([
    "0" = "0000", "1" = "0001", "2" = "0010", "3" = "0011", "4" = "0100",
    "5" = "0101", "6" = "0110", "7" = "0111", "8" = "1000", "9" = "1001",
    "A" = "1010", "B" = "1011", "C" = "1100", "D" = "1101", "E" = "1110", "F" = "1111"],
        tmp),""); end proc:

processpacket := proc( bitinput )
local vertotal, ver, pid, rest, parts, lofbits, subpak, numpaks,
    n, val, theops, thevalue;

    if length(bitinput) < 6 then
        error "error %1 is too short to be valid", bininput;
    end if;

    PRINT("Processing %d bits", length(bitinput));

    ver := binstr2dec(bitinput[1..3]);
    pid := binstr2dec(bitinput[4..6]);
    rest := bitinput[7..];
    vertotal := ver;

    if pid = 4 then # literal
        PRINT("processing literal");
        parts := NULL;
        while rest[1] <> "0" do
            parts := parts, rest[2..5];
            rest := rest[6..];
        end do;
        parts := parts, rest[2..5];

        thevalue := binstr2dec(cat(parts));
        PRINT("literal had %d parts with value=%d", nops([parts]), thevalue);
        rest := rest[6..];
        PRINT("after literal %d bits left", length(rest));
    else
        PRINT("processing operator ID=%d", pid);
        theops := NULL;
        if rest[1] = "0" then

            lofbits := binstr2dec(rest[2..16]);
            PRINT("type-0 operator with %d bits of subpackets", lofbits);
            rest := rest[17..];
            subpak := rest[1..lofbits];
            while length(subpak) > 6 do
                (n, val, subpak) := thisproc(subpak);
                vertotal += n;
                theops := theops, val;
            end do;
            rest := rest[lofbits+1..];
            PRINT("after type-0 %d bits left", length(rest));
        else
            numpaks := binstr2dec(rest[2..12]);
            PRINT("type-1 operator with %d subpackets",numpaks);
            rest := rest[13..];
            to numpaks do
                (n, val, rest) := thisproc(rest);
                theops := theops, val;
                vertotal += n;
            end do;
            PRINT("after type-1 %d bits left", length(rest));
        end if;
        if pid = 0 then thevalue := `+`(theops);
        elif pid = 1 then thevalue := `*`(theops);
        elif pid = 2 then thevalue := min(theops);
        elif pid = 3 then thevalue := max(theops);
        elif pid = 5 then thevalue := ifelse( theops[1] > theops[2], 1, 0);
        elif pid = 6 then thevalue := ifelse( theops[1] < theops[2], 1, 0);
        elif pid = 7 then thevalue := ifelse( theops[1] = theops[2], 1, 0);
        else error "bad operator pid";
        end if;
    end if;

    PRINT("returning with %d bits left", length(rest));
    return vertotal, thevalue, rest;

end proc:

# uncomment for verbose output
#PRINT := s->printf(cat(s,"\n"), _rest);

(*
# unit tests
processpacket( hex2bits("C200B40A82") )[2]=3;
processpacket( hex2bits("04005AC33890") )[2]=54;
processpacket( hex2bits("880086C3E88112"))[2]=7;
processpacket( hex2bits( "CE00C43D881120" ) )[2]=9;
processpacket( hex2bits( "D8005AC2A8F0" ))[2]=1;
processpacket( hex2bits("F600BC2D8F") )[2]=0;
processpacket( hex2bits("9C005AC2F8F0") )[2] = 0;
processpacket( hex2bits("9C0141080250320F1802104A08"))[2] = 1;
*)

# answers
(answer1, answer2, fluff) := processpacket( hex2bits( input ) ):
answer1;
answer2;


text:
######################################################################

`eta/ord_simplex_interior` := proc(A::set)
 if nops(A) <> 1 then return FAIL; fi;

 return [`eta/ord`(A),`eta/simplex_interior`(A)];
end;

`gamma/ord_simplex_interior` := (A::set,B::set) -> (p) -> proc(U,V)
 local RU,RV,xU,xV,b;

 RU := U[1];
 RV := table([seq(b = eval(V[b][1]),b in B)]);
 xU := U[2];
 xV := table([seq(b = eval(V[b][2]),b in B)]);
 return [`gamma/ord`(A,B)(p)(RU,RV),
         `gamma/simplex_interior`(A,B)(p)(xU,xV)];
end;


text:
with(BinomSums):
with(CodeTools):

### sumtoct

S := Sum(Binomial(n,k),k=0..n);
R := sumtoct(S, u);

Test(residue(subs(n=10,R)/u[1], u[1]=0), 2^10);

### sumtores

S := Sum(t^n*Sum(Binomial(n,k),k=0..n), n=0..infinity);
R, ord := sumtores(S, u);

Test(nops(ord), 1);
Test(R, 1/(1-2*t), `normal`);

max_stars_repo_name: lairez/binomsum

max_issues_repo_name: lairez/binomsum

max_forks_repo_name: lairez/binomsum


text:
`permutation_matrix` := proc(s)
  local n,g,i;

  n := nops(s);
  g := Matrix(n,n);
  for i from 1 to n do
   g[s[i],i] := 1;
  od:
  return g;
end:

`elementary_matrix/D` := (n) -> proc(i,a)
  local g,k;
  g := Matrix(n,n);
  for k from 1 to n do g[k,k] := 1; od;
  g[i,i] := a;
  return g;
end:

`elementary_matrix/E` := (n) -> proc(i,j,a_)
  local g,k;
  g := Matrix(n,n);
  for k from 1 to n do g[k,k] := 1; od;
  g[i,j] := `if`(nargs > 2,a_,1);
  return g;
end:

`elementary_matrix/F` := (n) -> proc(i,j)
  local g,k;
  g := Matrix(n,n);
  for k from 1 to n do g[k,k] := 1; od;
  g[i,i] := 0; g[j,j] := 0;
  g[i,j] := 1; g[j,i] := 1;
  return g;
end:

make_general_linear := proc(n,p)
  local G,i,j,k,m,g,L,f;

  m := n^2;

  G := table():
  G["p"] := p;
  G["id"] := IdentityMatrix(n);
  G["o"] := (g,h) -> map(mods,g . h,G["p"]);
  G["inv"] := g -> map(mods,1/g,G["p"]);
  G["order"] := mul(p^n - p^i,i = 0..n-1);

  L := [[]];
  for i from 1 to m do
   L := [seq(seq([op(u),j],j=0..p-1),u in L)];
  od:

  f := proc(u)
   local i,j,k,g;
   g := Matrix(n,n);
   k := 1;
   for i from 1 to n do
    for j from 1 to n do
     g[i,j] := mods(u[k],p); k := k + 1;
    od:
   od:
   if mods(Determinant(g),p) <> 0 then
    return g;
   else
    return NULL;
   fi;
  end:

  G["elements"] := map(f,L);

  return eval(G);
end:

bruhat_permutation := (n,p) -> proc(g)
 local h,i,j,k,a,q,qi;

 if [Dimension(g)] <> [n,n] then
  error("Not an n x n matrix");
 fi;

 h := eval(g);
 q := [];
 for i from n to 1 by -1 do
  j := 1;
  while h[i,j] = 0 do j := j+1; od;
  q := [j,op(q)];
  a := mods(1/h[i,j],p);
  h := map(mods,`elementary_matrix/D`(n)(i,a) . h,p);
  for k from 1 to i-1 do
   h := map(mods,`elementary_matrix/E`(n)(k,i,-h[k,j]) . h,p);
  od;
 od;
 qi := table():
 for i from 1 to n do qi[q[i]] := i; od:

 return([seq(qi[i],i=1..n)]);
end:

# This returns a triple [b,w,c] where b is upper triangular,
# w is a permutation matrix, c is upper unitriangular and g = bwc.

bruhat_decomposition := (n,p) -> proc(g)
 local h,i,j,k,a,b,c,x,q,w;

 if [Dimension(g)] <> [n,n] then
  error("Not an n x n matrix");
 fi;

 h := eval(g);
 q := [];
 b := Matrix(n,n);
 for k from 1 to n do b[k,k] := 1; od:

 for i from n to 1 by -1 do
  j := 1;
  while h[i,j] = 0 do j := j+1; od;
  q := [j,op(q)];
  x := mods(1/h[i,j],p);
  b := map(mods, b . `elementary_matrix/D`(n)(i,h[i,j]),p);
  h := map(mods, `elementary_matrix/D`(n)(i,x) . h,p);
  for k from 1 to i-1 do
   b := map(mods,b . `elementary_matrix/E`(n)(k,i,h[k,j]),p);
   h := map(mods,`elementary_matrix/E`(n)(k,i,-h[k,j]) . h,p);
  od;
 od;
 w := permutation_matrix(q);
 c := w . h;
 return([b,Transpose(w),c]);
end:


text:
grassmann_count := (q) -> proc(n::nonnegint,d::nonnegint)
 local i;
 return factor(mul(q^(n-i)-1,i=0..d-1)/mul(q^(d-i)-1,i=0..d-1));
end:

splitting_count := (q) -> (n::nonnegint,m::nonnegint) ->
 grassmann_count(q)(n+m,n) * q^(n*m);

`indices/twisted_product/A` := proc(r::nonnegint)
 option remember;
 local U,V,u,v,e;

 if r = 0 then return [[]]; fi;
 if not(type(r,posint)) then return FAIL; fi;

 U := `indices/twisted_product/A`(r-1);
 V := NULL;
 for u in U do
  e := d-1-`+`(op(u));
  V := V,seq([op(u),v],v=0..e);
 od;
 return [V];
end:

`is_element/twisted_product/A` := (r::nonnegint) -> proc(f)
 if not(type(f,table)) then
  return false;
 fi;

 return evalb({indices(f)} = {op(`indices/twisted_product/A`(r))});
end:

`is_equal/twisted_product/A` := (r::nonnegint) -> proc(f,g)
 local U,u;

 U := `indices/twisted_product/A`(r);

 for u in U do
  if simplify(f[op(u)]-g[op(u)]) <> 0 then
   return false;
  fi;
 od;

 return true;
end:

`to_list/twisted_product/A` := (r::nonnegint) -> proc(f)
 local u;
 return [seq(f[op(u)],u in `indices/twisted_product/A`(r))];
end:

`from_symbol/twisted_product/A` := (r::nonnegint) -> proc(f)
 local u;
 return table([seq(op(u)=f[op(u)],u in `indices/twisted_product/A`(r))]);
end:

`zero/twisted_product/A` := proc(r::nonnegint)
 option remember;
 local u;
 return table([seq(op(u)=0,u in `indices/twisted_product/A`(r))]);
end:

`one/twisted_product/A` := proc(r::nonnegint)
 option remember;
 local u;
 return table([seq(op(u)=1,u in `indices/twisted_product/A`(r))]);
end:

`delta/twisted_product/A` := proc(r::nonnegint)
 option remember;
 local u;
 return table([seq(op(u)=`if`(u=[0$r],1,0),u in `indices/twisted_product/A`(r))]);
end:

`plus/twisted_product/A` := (r::nonnegint) -> proc(f,g)
 local u;
 return table([seq(op(u) = f[op(u)] + g[op(u)],u in `indices/twisted_product/A`(r))]);
end:

`minus/twisted_product/A` := (r::nonnegint) -> proc(f,g)
 local u;
 return table([seq(op(u) = f[op(u)] - g[op(u)],u in `indices/twisted_product/A`(r))]);
end:

`dot/twisted_product/A` := (r::nonnegint) -> proc(f,g)
 local u;
 return table([seq(op(u) = f[op(u)] * g[op(u)],u in `indices/twisted_product/A`(r))]);
end:

`scalar_times/twisted_product/A` := (r::nonnegint) -> proc(c,f)
 local u;
 return table([seq(op(u) = c * f[op(u)],u in `indices/twisted_product/A`(r))]);
end:

`mu/twisted_product/A` := (r::nonnegint) -> proc(u,v)
 option remember;
 local i,j;
 mul(splitting_count(q)(u[i],v[i]),i=1..r) *
 q^add(add(u[i]*v[j]+2*u[j]*v[i],j=i+1..r),i=1..r-1);
end:

`star/twisted_product/A` := (r::nonnegint) -> proc(f,g)
 local h,X,u,v,w,x;

 h := table([seq(op(u)=0,u in `indices/twisted_product/A`(r))]);

 X := `indices/twisted_product/A`(2*r);
 for x in X do
  u := [op(1..r,x)];
  v := [op(r+1..2*r,x)];
  w := u +~ v;
  h[op(w)] := h[op(w)] + `mu/twisted_product/A`(r)(u,v) * f[op(u)] * g[op(v)];
 od;

 return eval(h);
end:


text:
######################################################################

`is_element/PP` := (N::posint) -> (A::set) -> proc(R)
 local E,i,a,b;
 global reason;

 if not(type(R,list) and nops(R) = N) then
  reason := [convert(procname,string),"R is not a list of length N",R,N];
  return false;
 fi;

 E := {seq(seq([a,b],b in A),a in A)};

 for i from 1 to N do
  if not `is_element/preord`(A)(R[i]) then
   reason := [convert(procname,string),"R[i] is not a preorder on A",R,i,A,reason];
   return false;
  fi;

  E := E intersect R[i] intersect `op/preord`(A)(R[i]);
 od;

 if E <> {seq([a,a],a in A)} then
  reason := [convert(procname,string),"E <> Delta",E];
  return false;
 fi;

 return true;
end:

`is_equal/PP` := (N::posint) -> (A::set) -> proc(R,S)
 local i;
 global reason;

 for i from 1 to N do
  if not(`is_equal/preord`(A)(R[i],S[i])) then
   reason := [convert(procname,string),"R[i] <> S[i]",i,R[i],S[i]];
   return false;
  fi;
 od;

 return true;
end:

`is_leq/PP` := (N::posint) -> (A::set) -> proc(R,S)
 local i;
 global reason;

 for i from 1 to N do
  if not(`is_leq/preord`(A)(R[i],S[i])) then
   reason := [convert(procname,string),"R[i] is not a subset of S[i]",i,R[i],S[i]];
   return false;
  fi;
 od;

 return true;
end:

`mu/F/PP` := (N::posint) -> (A::set) -> proc(x)
 local AA,a,b,i;

 AA := {seq(seq([a,b],b in A),a in A)};

 [seq(select(ab -> x[ab[1]][i] <= x[ab[2]][i],AA),i=1..N)];
end:

`sigma/PP/F` := (N::posint) -> (A::set) -> proc(R)
 local x0,x,r,a,i;

 x0 := table();
 for i from 1 to N do
  r := eval(`rank_table/preord`(A)(R[i]));
  for a in A do x0[a,i] := r[a]; od;
 od;

 x := table();
 for a in A do x[a] := [seq(x0[a,i],i=1..N)]; od;

 return eval(x);
end;

`random_element/PP` := (N::posint) -> (A::set) -> proc()
 local R,ok,i;

 ok := false;
 while not(ok) do
  R := [seq(`random_element/preord`(A),i=1..N)];
  ok := `is_element/PP`(N)(A)(R);
 od;

 return R;
end;

`list_elements/PP` := (N::posint) -> proc(A::set)
 local RR,QQ,R,Q,i;

 RR := `list_elements/preord`(A);
 QQ := [[]];
 for i from 1 to N do
  QQ := [seq(seq([op(Q),R],R in RR),Q in QQ)];
 od;

 QQ := select(`is_element/PP`(N)(A),QQ);

 return QQ;
end;

######################################################################

`is_element/realisation/PP` := (N::posint) -> (A::set) ->
 `is_element/realisation/generic`("PP",
  eval(`is_element/PP`(N)(A)),
  eval(`is_leq/PP`(N)(A))
 );

`is_equal/realisation/PP` := (N::posint) -> (A::set) ->
 `is_equal/realisation/generic`("PP",
  eval(`is_equal/PP`(N)(A))
 );


text:
program parseformalparameters;
procedure func(a, b : integer;
    c, d : char);
    begin
        a := b;
        c := d;
    end;
var x1, x2 : integer;
var y1, y2 : char;
begin
    call func(x1, x2, y1, y2);
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing

max_forks_repo_name: Hiroya-W/lang-processing


text:
read "../ComputeIdentifiableFunctions.mpl";

cases := [
    [
        [
            diff(x1(t), t) = (a * x1(t) - b * x1(t) * x2(t) + u(t)),
            diff(x2(t), t) = (-c * x2(t) + d * x1(t) * x2(t)),
            y(t) = x1(t)
        ], [1, {a, c, d}]
    ],
    [
        [
            diff(xa(t), t) = k1 * xa(t),
            diff(xb(t), t) = k1 * xa(t) - k2 * xb(t),
            diff(xc(t), t) = k2 * xb(t),
            diff(eA(t), t) = 0,
            diff(eC(t), t) = 0,
            y1(t) = eA(t) * xa(t) + eB * xb(t) + eC(t) * xc(t),
            y2(t) = xc(t),
            y3(t) = eA(t),
            y4(t) = eC(t)
        ], [2, {k1, k2, eB}]
    ]
]:

num_passed := 0:
num_failed := 0:

for case in cases do
    input := case[1]:
    correct := case[2]:
    result := MultiExperimentIdentifiableFunctions(input, infolevel=0, simplified_generators=true):
    if CompareFields(result[3], correct[2]) and result[1] = correct[1] then
        printf("PASSED\n");
        num_passed := num_passed + 1:
    else
        num_failed := num_failed + 1:
        printf("FAILED\n");
        printf("Expected %a, got %a", correct, result);
    end if:
end do:

printf("Passed: %a, failed %a \n", num_passed, num_failed);

max_stars_repo_name: iliailmer/AllIdentifiableFunctions

max_issues_repo_name: iliailmer/AllIdentifiableFunctions

max_forks_repo_name: iliailmer/AllIdentifiableFunctions


text:
repeat process ACV_lightBridges {
// Remove unsupported bridges
/kill @e[type=area_effect_cloud,name=ACV_LightBridge,rxm=-90,rx=-1]
/kill @e[type=area_effect_cloud,name=ACV_LightBridge,rxm=1,rx=90]

/scoreboard players tag @e[type=area_effect_cloud,name=ACV_LightBridge,tag=ACV_InAir] remove ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_LightBridge] ~ ~ ~ detect ~ ~ ~ air * scoreboard players tag @e[type=area_effect_cloud,name=ACV_LightBridge,r=0,c=1] add ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_LightBridge] ~ ~ ~ detect ~ ~ ~ carpet 3 scoreboard players tag @e[type=area_effect_cloud,name=ACV_LightBridge,r=0,c=1] add ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_LightBridge] ~ ~ ~ detect ~ ~ ~ piston_extension * execute @e[type=area_effect_cloud,name=ACV_LightBridge,r=0,c=1] ~ ~ ~ detect ~ ~1 ~ carpet 3 scoreboard players tag @e[type=area_effect_cloud,name=ACV_LightBridge,r=0,c=1] add ACV_InAir

/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=!ACV_InAir] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_PortalIn,dy=0] add ACV_createLightBridges_PARAM
/scoreboard players tag @e[type=area_effect_cloud,name=ACV_LightBridge,tag=!ACV_InAir] add ACV_TpProof
/entitydata @e[type=area_effect_cloud,name=ACV_LightBridge,tag=!ACV_InAir] {CustomName:"ACV_LightBridgeEnd"}
conditional: start ACV_createLightBridges

// To avoid race conditions, if there are multiple lightbridges in the same block one needs to be processed (setblock) before the others, so that the others summon all markers correctly
/scoreboard players tag @e[type=area_effect_cloud,name=ACV_LightBridge] add ACV_Multiple
/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=ACV_Multiple] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LightBridge,dy=0,c=-1] remove ACV_Multiple

/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=!ACV_Multiple] ~ ~ ~ detect ~ ~-1 ~ carpet 3 summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_BridgeBelow",Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=!ACV_Multiple] ~ ~ ~ detect ~ ~ ~ carpet 3 summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_BridgeCross",Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=!ACV_Multiple] ~ ~ ~ detect ~ ~1 ~ carpet 3 summon area_effect_cloud ~ ~1 ~ {CustomName:"ACV_BridgeBelow",Duration:2147483647}

/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=!ACV_Multiple] ~ ~ ~ setblock ~ ~-1 ~ piston_extension 0 keep
/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=!ACV_Multiple] ~ ~ ~ setblock ~ ~ ~ carpet 3

/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=ACV_Multiple] ~ ~ ~ detect ~ ~-1 ~ carpet 3 summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_BridgeBelow",Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=ACV_Multiple] ~ ~ ~ detect ~ ~ ~ carpet 3 summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_BridgeCross",Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_LightBridge,tag=ACV_Multiple] ~ ~ ~ detect ~ ~1 ~ carpet 3 summon area_effect_cloud ~ ~1 ~ {CustomName:"ACV_BridgeBelow",Duration:2147483647}

/tp @e[type=area_effect_cloud,name=ACV_LightBridge,rym=45,ry=135] ~-1 ~ ~
/tp @e[type=area_effect_cloud,name=ACV_LightBridge,rym=135,ry=-135] ~ ~ ~-1
/tp @e[type=area_effect_cloud,name=ACV_LightBridge,rym=-135,ry=-45] ~1 ~ ~
/tp @e[type=area_effect_cloud,name=ACV_LightBridge,rym=-45,ry=45] ~ ~ ~1

// Wake up all idle light bridges
/entitydata @e[type=area_effect_cloud,name=ACV_LightBridge_Idle] {CustomName:"ACV_LightBridge"}
/scoreboard players tag @e[type=area_effect_cloud,name=ACV_LightBridge] add ACV_TpProof

/testfor @e[type=area_effect_cloud,name=ACV_LightBridge]
invert: stop
}


text:
program typeconv;
var a : array[3] of integer;
var b : integer;
var c : char;
var d : boolean;
begin
    b := integer(a[0]);
    c := char(a[0]);
    d := boolean(a[0]);
end.


text:
###############################################################################
# Part 1: Algorithms for computation with subfields of rational functions
###############################################################################

with(PolynomialIdeals):
#------------------------------------------------------------------------------

IdealsEq := proc(A, B)
    # Checks whether polynomial ideals are equal
    return evalb((A subset B) and (B subset A)):
end:

#------------------------------------------------------------------------------

FieldToIdeal := proc(gens)
    # Input: generators of a subfield of the field of rational functions
    # Computes the MQS ideal of the field with the new variables of the form x_aux
    # See: https://mediatum.ub.tum.de/doc/685465/document.pdf Definition 2.16
    #      https://doi.org/10.1016/j.jsc.2005.09.010          Lemma 2.1
    local all_vars, subs_dupl, all_dupl, common_denom, polys, f, gb:
    all_vars := indets(gens):
    subs_dupl := map(v -> v = cat(v, _aux), all_vars):
    all_dupl := map(v -> subs(subs_dupl, v), all_vars):
    common_denom := 1:
    polys := []:
    for f in gens do
        common_denom := lcm(common_denom, denom(f)):
        polys := [op(polys), numer(f) * subs(subs_dupl, denom(f)) - subs(subs_dupl, numer(f)) * denom(f)]:
    end do:
    #gb := Groebner[Basis]([op(polys), common_denom * t - 1], plex(t, op(all_dupl))):

    gb := Groebner[Basis]([op(polys), common_denom * t - 1], tdeg(t, op(all_dupl))):
    gb := Groebner[Walk](gb, tdeg(t, op(all_dupl)), lexdeg([t], [op(all_dupl)])):

    gb := select(p -> not (t in indets(p)), gb):
    return PolynomialIdeal(gb, variables=all_dupl):
end proc:

#------------------------------------------------------------------------------

FieldCoefficientRestriction := proc(J, msq_for_subfield)
    # Input: J - a polynomial ideals over a field of rational functions
    #        msq_for_subfield - the MQS ideal for a subfield E of coefficients (see FieldToIdeal)
    # Computes the radical of the restriction of the ideal to the subfield E
    # in the sense of https://doi.org/10.1016/j.jsc.2005.09.010 (MSQ-paper in what follows)
    #
    # NOTE: unlike the algorithm in the MQS-paper, we compute the radical, not the restriction itself
    # one can obtain the algorithm MQS-paper by replacing PrimeDecomposition with PrimaryDecomposition
    # in the code below
    local poly_vars, coeff_vars, subs_aux, coeff_aux, gens, subs_aux_msq, gens_msq, msq_ideal_aux,
    msq_components, J_ext, components, primes_to_keep, P, elim_P, comp, cleaned_ideal:

    poly_vars := IdealInfo[Variables](J):
    coeff_vars := IdealInfo[Parameters](J) union IdealInfo[Parameters](msq_for_subfield):
    subs_aux := map(v -> v = parse(cat(v, _aux_aux)), coeff_vars):
    coeff_aux := subs(subs_aux, coeff_vars):

    # list F of polynomials in the notation of the MSQ-paper, page 375
    gens := map(p -> subs(subs_aux, p * denom(p)), IdealInfo[Generators](J)):

    # Substitution to avoid names clashing between the aux variables in the msq ideal and the variables in J
    subs_aux_msq := map(v -> v = parse(cat(v, _aux)), IdealInfo[Variables](msq_for_subfield)):
    gens_msq := map(p -> subs(subs_aux_msq, p), IdealInfo[Generators](msq_for_subfield)):
    msq_ideal_aux := PolynomialIdeal(gens_msq, variables=map(s -> rhs(s), subs_aux_msq)):
    msq_components := [PrimeDecomposition(msq_ideal_aux)]:

    J_ext := PolynomialIdeal([op(gens), op(gens_msq)], variables=poly_vars union coeff_aux):
    components := [PrimeDecomposition(J_ext)]:

    # Selecting prime components as in Remark on page 377 in MSQ-paper
    primes_to_keep := []:
    for P in components do
        # Going through the prime decomposition of the MSQ-deal mimics the
        # working over the restricted field (which is what has been done in the MSQ-paper)
        elim_P := EliminationIdeal(P, coeff_aux):
        for comp in msq_components do
            if elim_P subset comp then
                primes_to_keep := [op(primes_to_keep), P]:
            end if:
        end do:
    end do:
    if nops(primes_to_keep) > 0 then
        cleaned_ideal := Intersect(op(primes_to_keep)):
    else
        cleaned_ideal := PolynomialIdeal([0], variables = poly_vars):
    end if:

    # Applying Lemma 2.2 from the MSQ-paper
    return EliminationIdeal(cleaned_ideal, poly_vars):
end proc:

#------------------------------------------------------------------------------

FilterGenerators := proc(ideal)
    # Input: ideal over a rational function field
    # Computes a simplified set of generators of the field of definition
    local gb, gens, p, cf, lc, gsorted, result, big_ideal, cur_ideal, g, new_ideal:
    gb := Groebner[Basis](ideal, tdeg(op(IdealInfo[Variables](ideal)))):
    gens := {}:
    for p in gb do
        cf := sort([coeffs(p, IdealInfo[Variables](ideal))], (a, b) -> length(convert(a, string)) < length(convert(b, string))):
        lc := cf[1]:
        cf := map(c -> c / lc, cf):
        gens := {op(gens), op(cf[2..nops(cf)])}:
    end do:
    gsorted := sort([op(gens)], (a, b) -> length(convert(a, string)) < length(convert(b, string)));
    result := {}:
    big_ideal := FieldToIdeal(gens):
    cur_ideal := FieldToIdeal(result):
    for g in gsorted do
        if big_ideal = cur_ideal then
            return result:
        end if:
        new_ideal := FieldToIdeal({op(result), g}):
        if new_ideal <> cur_ideal then
        	  # a dirty hack to transform -1/a to a
            if convert(g, string)[1] = "-" then
                g := -g:
            end:
            if convert(g, string)[1] = "1" then
                g := 1 / g:
            end:
            result := {op(result), g}:
            cur_ideal := new_ideal:
        end if:
    end do:
    return result:
end proc:

#------------------------------------------------------------------------------

FieldIntersection := proc(gens_left, gens_right)
    # Input: gens_left and gens_right - generators of a subfields E and F of a field of rational functions
    # If terminates, resturns an ideal with field of definition being the intersection of generators of E and F
    # Is guaranteed to terminate if at least one of E and F is algebraically closed in rational functions (see REF)
    # Implementation below is a version of Algorithm 2.38 from https://mediatum.ub.tum.de/doc/685465/document.pdf
    local msq_left, msq_right, poly_vars, coeff_vars, Ii, Ji, gb, result, p, cf, lc;

    msq_left := FieldToIdeal(gens_left):
    msq_right := FieldToIdeal(gens_right):
    poly_vars := IdealInfo[Variables](msq_left) union IdealInfo[Variables](msq_right):
    coeff_vars := IdealInfo[Parameters](msq_left) union IdealInfo[Parameters](msq_right):

    Ii := PolynomialIdeal([1], variables=poly_vars):
    Ji := msq_left:

    while not IdealsEq(Ii, Ji) do
        Ii := FieldCoefficientRestriction(Ji, msq_right):
        Ji := FieldCoefficientRestriction(Ii, msq_left):
    end do:

    return Ji:
end proc:

#------------------------------------------------------------------------------

CompareFields := proc(gens_l, gens_r)
    # Checks whether gens_l and gens_r generate the same subfield of rational functions
    return IdealsEq(FieldToIdeal(gens_l), FieldToIdeal(gens_r)):
end proc:

#------------------------------------------------------------------------------

###############################################################################
# Part 2: Algorithm for computing  identifiable functions
###############################################################################

with(DifferentialAlgebra):

#===============================================================================

ExtractDenominator := proc(model)
    # Input: model a list of rational functions
    # returns the function multiplied by their denominators and
    # an inequality corresponding to the LCM of the denominators
    local common_denom, r;
    common_denom := 1:
    for r in model do
        common_denom := lcm(common_denom, denom(r)):
    end do:
    return [op(map(p -> denom(p) * p, model)), common_denom <> 0]:
end proc:

#===============================================================================
CheckReducibilitySet := proc(polys, charset)
    # Input: polys - list of differential polynomials
    #        charset - a characteristic set
    # Returns true iff all the polynomials are reduced to zero wrt to the charset
    local result, e:
    result := true:
    for e in polys do
        if NormalForm(e, charset) <> 0 then
            result := false:
            break:
        end if:
    end do:
    return result:
end proc:
#===============================================================================

GetIOEquations := proc(model, states, inputs, outputs, params, use_brackets, infolevel, target)
    # Input: model - list of differential polynomials defining the model
    #        states, ios, params - list of names of state variables, input-output variables
    #                              and parameter, respectively
    # Computes a list of input-output equations of the model
    local Relim, Rorig, charsets, chset_orig, general_comps, general, c, e, gen_comp, io_eqs:
    Relim := DifferentialRing(blocks = [[op(states)], [op(outputs)], [op(inputs)]], derivations = [t], arbitrary = params):
    if not use_brackets then
        Relim := DifferentialRing(blocks = [[op(states)], op(outputs), [op(inputs)]], derivations = [t], arbitrary = params):
    end if:
    Rorig := DifferentialRing(blocks = [[op(outputs)], [op(states)], [op(inputs)]], derivations = [t], arbitrary = params):
           # DifferentialRing(blocks = [[op(inputs)], [op(outputs)], [op(states)]], derivations = [t], arbitrary = params):
    chset_orig := RosenfeldGroebner(model, Rorig)[1]:

    if infolevel > 0 then
        LogText(sprintf("    Computing the characteristic set (singsol = none)\n"), target):
    end if:
    charsets := RosenfeldGroebner(model, Relim, singsol = none):
    if CheckReducibilitySet(Equations(charsets[1]), chset_orig) then
        gen_comp := charsets[1]:
    else
        if infolevel > 0 then
            LogText(sprintf("    Did not pick the right component. Using singsol = all\n"), target):
        end if:
        charsets := RosenfeldGroebner(model, Relim):
        if infolevel > 0 then
            LogText(sprint("    Selecting the general component\n"), target):
        end if:
        general_comps := []:
        for c in charsets do
            general := true:
            for e in Equations(c) do
                if NormalForm(e, chset_orig) <> 0 then
                    general := false:
                    break:
                end if:
            end do:
            if general then
                general_comps := [op(general_comps), c]:
            end if:
        end do:
        if nops(general_comps) > 1 then
            LogText(sprintf("More than one component is considered general!", general_comps), target):
        end if:
        gen_comp := general_comps[1]:
    end if:
    io_eqs := Equations(gen_comp, leader < parse(cat(states[-1], "(t)"))):
    return io_eqs:
end proc:

#===============================================================================

# Adapted from
# https://www.mapleprimes.com/questions/36772-Extract-Specific-Coefficients-Of-A-Multivariate
# by Kitonum 15364
coefff:=proc(P, t)
    local L, H, i, k:
    L:=[coeffs(P, indets(P), 'h_aux_for_coef')]: H:=[h_aux_for_coef]: k:=0:
    for i from 1 to nops(H) do
        if H[i]=t then k:=L[i] fi:
    end do:
    return k;
end proc:

#===============================================================================

DecomposePolynomial := proc(p, vars_main, vars_coef, infolevel, target)
    # Input: p - polynomial in two groups of variables: vars_main and vars_coef
    # Computes a decomposition of minimal length of p as a linear combination
    # of products A * B, where A is a polynomial in vars_main and B
    # is a polynomial in vars_coef return two lists: list of A's and list of B's
    local cf, monoms, result_cf, result_monom, i, c, m, j, lc, lm, coeff_in_c:
    cf := [coeffs(collect(p, vars_main, 'distributed'), vars_main, 'monoms')]:
    monoms := [monoms]:
    result_cf := []:
    result_monom := Array([]):
    if infolevel > 0 then
        LogText(sprintf("        Number of monomials: %a\n", nops(cf)), target):
    end:
    for i from 1 to nops(cf) do
        c := cf[i]:
        m := monoms[i]:
        for j from 1 to nops(result_cf) do
            lc, lm := Groebner[LeadingTerm](result_cf[j], plex(op(vars_coef))):
            coeff_in_c := coefff(c, lm):
            c := c - coeff_in_c / lc * result_cf[j]:
            result_monom[j] := result_monom[j] + coeff_in_c / lc * m:
        end do:
        if c <> 0 then
            result_cf := [op(result_cf), c]:
            ArrayTools[Append](result_monom, m):
        end if:
    end do:
    if infolevel > 0 then
        LogText(sprintf("        Reduced to: %a\n", nops(result_cf)), target):
    end:
    return [result_cf, convert(result_monom, list)]:
end proc:

#===============================================================================

ConstructWronskian := proc(io_eq, model, states, inputs, outputs, params, state_space, infolevel, target)
    # Input - the same as for GetIOEquations + one IO-equation + flag subs_param
    # Computes the Wronskian for this equation using the representation
    # given by DecomposePolynomial. Return a pair of the Wronskian
    # reduced modulo the original system and a list of coefficients
    # of the compressed io_eq
    local getNormalForm, diff_to_ord, v, vt, h, v_ord, vd, p, decomp, diff_polys, Rorig, chset_orig,
    M, yus, yus_reduced, yus_list, M_sub, roll, params_sub, ios, ps_sol, ps, i:

    diff_to_ord := {}:
    ios := [op(inputs), op(outputs)]:
    for v in ios do
        vt := parse(cat(v, "(t)")):
        diff_to_ord := {op(diff_to_ord), vt = v}:
        for h from 1 to nops(states) + 1 do
            v_ord := parse(cat(v, "_", h)):
            vd := diff(vt, t$h):
            diff_to_ord := {op(diff_to_ord), vd = v_ord}:
        end do:
    end do:
    p := subs(diff_to_ord, io_eq):
    if infolevel > 0 then
        LogText(sprintf("    Combining monomials to reduce the dimension\n"), target):
    end if:
    decomp := DecomposePolynomial(p, map(e -> rhs(e), diff_to_ord), params, infolevel, target):
    diff_polys := map(p -> subs(map(e -> rhs(e) = lhs(e), diff_to_ord), p), decomp[2]):
    Rorig := DifferentialRing(blocks = [[op(outputs)], [op(states)], [op(inputs)]], derivations = [t], arbitrary = params):
           # DifferentialRing(blocks = [[op(inputs)], [op(outputs)], [op(states)]], derivations = [t], arbitrary = params):
    chset_orig := RosenfeldGroebner(model, Rorig)[1]:

    if infolevel > 0 then
        LogText(sprintf("    Computing the Wronskian\n"), target):
    end if:
    M := VectorCalculus[Wronskian](diff_polys, t):
    yus := indets(M) minus {t}:

    if infolevel > 0 then
        LogText(sprintf("    Reducing the Wronskian\n"), target):
    end if:
    if state_space <> [] then
        ps_sol := GetPSSolution(state_space, nops(yus)):
        yus_reduced := []:
        for v in ios do
            vt := parse(cat(v, "(t)")):
            ps := subs(ps_sol, v(t)):
            for i from 0 to nops(yus) - 1 do
                yus_reduced := [op(yus_reduced), vt = subs({t = 0}, ps)]:
                vt := diff(vt, t):
                ps := diff(ps, t):
            end do:
        end do:
    else
        getNormalForm := proc (p, chset_orig)
            return p=DifferentialAlgebra:-NormalForm(p, chset_orig):
        end proc:
        Rorig := DifferentialRing(blocks = [[op(outputs)], [op(states)], [op(inputs)]], derivations = [t], arbitrary = params):
               # DifferentialRing(blocks = [[op(inputs)], [op(outputs)], [op(states)]], derivations = [t], arbitrary = params):
        chset_orig := RosenfeldGroebner(model, Rorig)[1]:
        yus_list := convert(yus, list):
        yus_reduced:=Threads:-Seq[tasksize=nops(yus_list)](getNormalForm(yus_list[i], chset_orig), i =1..nops(yus_list)):
        # yus_reduced := map(p -> p = NormalForm(p, chset_orig), yus):
    end if:

    # iofun:=map(y->parse(cat(convert(y, string), "(t)")), ios):
    # derivative_orders:=table([seq(p=0, p in iofun)]):
    # for each in yus_list do
    #     indets_ := indets(each) minus {t}: # get indets of derivative to compute order
    #     function_ := op(indets_ intersect {op(iofun)}): # for which function?
    #     order:=numelems(indets_): # compute order
    #     if derivative_orders[function_]<order then
    #         derivative_orders[function_]:=order;
    #     fi:
    # end do;
    # LogExpression(derivative_orders):
    # yus_reduced:=[]:
    # count:=0:
    # # yus_reduced:=Threads:-Seq[tasksize=nops(yus_list)](getNormalForm(yus_list[i], chset_orig), i =1..nops(yus_list)):
    # for y in ios do
    #     y:=parse(cat(convert(y, string), "(t)")):
    #     y_rem:=y:
    #     fun:=y:
    #     while count<derivative_orders[fun] do
    #         yus_reduced:=[op(yus_reduced), y=DifferentialAlgebra:-NormalForm(y_rem, chset_orig)]:
    #         count := count+1:
    #         y_rem:=diff(y_rem, t):
    #         y:=diff(y, t):
    #         LogText("Done Processing %a, order: %a\n", fun, count):
    #     end do:
    #     count:=0:
    # end do:

    M_sub := subs(yus_reduced, M):
    M_sub := subs(map(x -> parse(cat(x, "(t)")) = x, states), M_sub):
    return [M_sub, decomp[1]]:
end proc:

#===============================================================================

SingleExperimentIdentifiableFunctions := proc(model, output_targets, {infolevel := 1})
    # Input: model - ODE model represented as a list of differential polynomials
    # Computes generators of the field of single-identifiable functions
    local result, start, finish, states, inputs, outputs, model_eq, ios, params, io_eqs, si_gens, eq, wrnsk, echelon_form, model_denomfree, target:
    states, inputs, outputs, params, model_eq := op(ParseInput_(model)): # states, ios, params := op(ParseInput_(model)):
    ios := [op(inputs), op(outputs)]:

    # Step 1
    if infolevel > 0 then
        LogText(sprintf("SE Step 1: Computing input-output equations\n"), output_targets[log]):
    end if:
    model_denomfree := ExtractDenominator(model_eq):
    io_eqs := GetIOEquations(model_denomfree, states, inputs, outputs, params, false, infolevel, output_targets[log]):
    if infolevel > 0 then
        LogText(sprintf("SE Total number of io-equations: %a\n", nops(io_eqs)), output_targets[log]):
    end if:

    	    si_gens := {}:
	    for eq in io_eqs do
	        # Step 2
	        if infolevel > 0 then
	            LogText(sprintf("SE Step 2: Constructing the Wronskian\n"), output_targets[log]):
	        end if:
	        wrnsk := ConstructWronskian(eq, model_denomfree, states, inputs, outputs, params, [], infolevel, output_targets[log])[1]:
	        # Step 3
	        if infolevel > 0 then
	            LogText(sprintf("SE Step 3: Computing the reduced row echelon form of the Wronskian\n"), output_targets[log]):
	        end if:
	        echelon_form := LinearAlgebra[ReducedRowEchelonForm](wrnsk):
	        si_gens := {op(si_gens), op(select(x -> not type(x, numeric), convert(echelon_form, list)))}:
	    end do:

    # Step 4
    if infolevel > 0 then
        LogText(sprintf("SE Step 4: restricting to the field of parameters"), output_targets[log]):
    end if:
    	result:=map(simplify, convert(FilterGenerators(FieldIntersection(si_gens, params)), list)):
	DocumentTools:-SetProperty(output_targets[single], expression, result, 'refresh'):
	return result:
end proc:

#===============================================================================
# Adapted from https://github.com/pogudingleb/SIAN
FunctionToVariable_ := proc(f):
    convert(convert(f, string)[1..-4], symbol):
end proc:

ParseInput_ := proc(model)
   local all_symbols, x_functions, io_functions, params, states, ios, diff_polys, lhss, out_functions, in_functions, inputs, outputs:
   diff_polys := map(eq -> lhs(eq) - rhs(eq), model):
   all_symbols := foldl(`union`, op( map(e -> indets(e), diff_polys) )) minus {t}:
   x_functions := map(f -> int(f, t), select( f -> type(int(f, t), function(name)), all_symbols )):
   io_functions := select( f -> not type(int(f, t), function(name)) and type(f, function(name)) and not f in x_functions, all_symbols ):
   lhss := map(eq -> lhs(eq), model):
   out_functions := select(f -> f in lhss, io_functions):
   in_functions := select(f -> not (f in lhss), io_functions):
   params := [op(select(f -> not type(f, function(name)) and not type(int(f, t), function(name)), all_symbols))]:
   states := [op(map(FunctionToVariable_, x_functions))]:
   inputs := [op(map(FunctionToVariable_, in_functions))]:
   outputs := [op(map(FunctionToVariable_, out_functions))]:
   return [states, inputs, outputs, params, diff_polys]:
end proc:

#===============================================================================

MultiExperimentIdentifiableFunctions := proc(model, simplified_generators, no_bound, simplify_bound, max_perms, output_targets): # {simplified_generators := false, no_bound := false})
    # Input: model - ODE model in the state-space form
    # Computes the bound from Theorem REF. Returns a list containing
    # 1. the bound
    # 2. the list of lists of coefficients of the IO-equations (after compression)
    # 3. (optional) simplified set of generators of the field of multi-experiment identifiable functions
    #
    # Optional keyword arguments:
    # 1. simplified_generators - if False, then just the coefficients of the IO-equations are returned.
    #                            if True, they are being simplified (maybe be a substantial simplification)
    #                            but this takes time
    # 2. no_bound - if True, then bound for the number of experiments is not computed (may save a lot of time)

    local states, inputs, outputs, ios, params, model_eqs, io_eqs, io_coeffs, io_coef, wrnsk, s, roll, wrnsk_sub, r, bound,
    generators, i, eq, result, model_denomfree, target, start, finish, infolevel, use_brackets, output_permutations,
    outputs_, io_coeffs_sb, io_eqs_sb, bound_sb, skip_simplify, result_sb, best_output_ordering, returned_generators, idx, old_bound:

    states, inputs, outputs, params, model_eqs := op(ParseInput_(model)):
    use_brackets:=false:
    output_permutations := combinat[permute](outputs):
    output_permutations:= output_permutations[..min(nops(output_permutations),max_perms)];

    infolevel := 1:
    model_denomfree := ExtractDenominator(model_eqs):
    skip_simplify := false:
    old_bound := 10: # this will be replaced by a smaller boud every time
    idx:=0:
    for use_brackets in [true, false] do
        for outputs in output_permutations do
        	LogText(sprintf("Permutation:\t%a\n", outputs),  output_targets[log]):
        	# the first iteration is default permutation
	    	   LogText(sprintf("Use Brackets?\t%a\n", use_brackets), output_targets[log]):
	        if infolevel > 0 then
	            LogText(sprintf("ME Computing input-output equations\n"), output_targets[log]):
	   	   end if:
	   	   io_eqs := GetIOEquations(model_denomfree, states, inputs, outputs, params, use_brackets, infolevel, output_targets[log]):
	   	   if infolevel > 0 then
	            LogText(sprintf("ME Total number of io-equations: %a\n", nops(io_eqs)), output_targets[log]):
	        end if:
	        io_coeffs := []:

	        for eq in io_eqs do
	            io_coef := DecomposePolynomial(eq, indets(eq) minus {op(params)}, params, infolevel, output_targets[log])[1]:
	            io_coeffs := [op(io_coeffs), io_coef]:
	        end do:

	        generators := {}:
    	        for io_coef in io_coeffs do
                 io_coef := sort(io_coef, (a, b) -> length(convert(a, string)) < length(convert(b, string)));
                 for i from 2 to nops(io_coef) do
                     generators := {op(generators), io_coef[i] / io_coef[1]}:
                 end do:
              end do:
              bound := 0:
		    if no_bound = false then
                  for eq in io_eqs do
                      if infolevel > 0 then
                          LogText(sprintf("ME Constructing the Wronskian\n"), output_targets[log]):
                      end if:
                      wrnsk, io_coef := op(ConstructWronskian(eq, model_denomfree, states, inputs, outputs, params, model, infolevel, output_targets[log])):
                      # in the notation of the theorem
                      s := nops(io_coef) - 1:
                      # substitution does not increase the rank, so the resulting bound will be correct
                      roll := rand(1..15):
                      wrnsk_sub := map(v -> v = roll(), indets(wrnsk)):
                      r := LinearAlgebra[Rank](subs(wrnsk_sub, wrnsk)):
                      bound := max(bound, s - r + 1):
                  end do:
               end if:
               if bound<old_bound or idx=0 then # idx makes sure we assign a new bound at the first run
               	old_bound:=bound:
               end if:
               idx := idx +1;
               LogText(sprintf("New Bound:\t%a\n", bound),  output_targets[log]):
               result := [old_bound, generators]:
               if simplified_generators then
                   if infolevel > 0 then
		             LogText(sprintf("ME WARNING: Entering simplification of generators! if this takes too long, try unchecking \"Simplified Generators\"\n"), output_targets[log]):
                   end if:
                   result := [op(result), FilterGenerators(FieldToIdeal(generators))]:
                   DocumentTools:-SetProperty(output_targets[multi], expression, map(simplify, convert(result[3], list)), 'refresh'):
                   returned_generators := result[3]:
               else
                   DocumentTools:-SetProperty(output_targets[multi], expression, map(simplify, convert(result[2], list)), 'refresh'):
                   returned_generators := result[2]:
               end if:
               if old_bound > 0 then
    	              DocumentTools:-SetProperty(output_targets[bound_], expression, old_bound, 'refresh'):
    	              if old_bound=1 then
    	  	             skip_simplify := true:
    	  	             if simplify_bound then
    	  		            LogText(sprintf("Bound on the number of experiments = 1, \"Try to refine bound\" was selected but it will be skipped"), output_targets[log]):
    	  	             end if:
	        	        if simplified_generators then
       		            DocumentTools:-SetProperty(output_targets[single], expression, map(simplify, convert(result[3], list)), 'refresh'):
   		             else
    	   		            DocumentTools:-SetProperty(output_targets[single], expression, map(simplify, convert(result[2], list)), 'refresh'):
   		             end if:
    	              end if:
               else
	              DocumentTools:-SetProperty(output_targets[bound_], expression, "",'refresh'):
               end if:
               if skip_simplify or not simplify_bound then
               	break:
              	else
              		DocumentTools:-SetProperty("being_refined", caption, "being refined", 'refresh'):
               end if:
         end do: # loop over outputs
         if skip_simplify or not simplify_bound then
             # only breaks if we do not simplify. the code runs at most once
             break:
         end if:
    end do:
    DocumentTools:-SetProperty("being_refined", caption, "", 'refresh'):
    return [bound, returned_generators];#table([bd=bound]):
end proc:

GetPSSolution := proc(model, ord)
    # Input: model and integer ord
    # Output: a truncated power series solution with random parameter values
    # and initial conditions of order ord
    local roll, states, inputs, outputs, params, model_eqs, model_sub,
    x_funcs, x_sols, x_eqs, cur_ord, i, rhs_eval, total_sub, y_funcs, y_sols,
    params_subs, input_subs:
    roll := rand(1..15):
    states, inputs, outputs, params, model_eqs := op(ParseInput_(model)):
    params_subs := map(p -> p = roll(), params):
    input_subs := map(u -> parse(cat(u, "(t)")) = add([seq(roll() * t^i, i=0..ord)]), inputs):
    model_sub := subs(params_subs, model):
    model_sub := subs(input_subs, model_sub):
    x_funcs := map(x -> parse(cat(x, "(t)")), states):
    x_sols := map(x -> x = roll(), x_funcs):
    x_eqs := map(x -> subs(model_sub, diff(x, t)), x_funcs):

    for cur_ord from 1 to ord do
        for i from 1 to nops(x_funcs) do
            rhs_eval := series(subs(x_sols, x_eqs[i]), t, ord + 1):
            x_sols[i] := (lhs(x_sols[i]) = (rhs(x_sols[i]) + t^cur_ord * coeff(rhs_eval, t, cur_ord - 1) / cur_ord)):
        end do:
    end do:

    total_sub := [op(x_sols), op(params_subs), op(input_subs)]:
    y_funcs := map(y -> parse(cat(y, "(t)")), outputs):
    y_sols := map(y -> y = subs(total_sub, subs(model_sub, y)), y_funcs):

    return [op(y_sols), op(input_subs), op(params_subs), op(x_sols)]:
end proc:

max_stars_repo_name: iliailmer/sian-web-app

max_issues_repo_name: iliailmer/sian-web-app

max_forks_repo_name: iliailmer/sian-web-app


text:
(*
 Copyright (C) 2020 Susi Lehtola

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)
(* prefix:
  gga_c_ccdf_params *params;

  assert(p->params != NULL);
  params = (gga_c_ccdf_params * )(p->params);
*)

(* Equation (26) *)
f_ccdf  := (rs, z, xt, xs0, xs1) ->
  params_a_c1 / (1 + params_a_c2*n_total(rs)^(-1/3)) * (1 - params_a_c3 / (1 + exp(-params_a_c4*(2^(1/3)*X2S*xt - params_a_c5)))):

f  := (rs, z, xt, xs0, xs1) -> f_ccdf(rs, z, xt, xs0, xs1):


text:
NAME program var array of begin end if then else procedure return
call break
while do not or div and char integer boolean readln
writeln true false NUMBER STRING + - * = <>
< <= > >= ( ) [ ] := .
, : ;
1 2 3 4 5 6 7 8 9 0 -1 -2 -3 -4 -5 -6 -7 -8 -9
'string'  '''s''t''r''i''n''g''''''s'''  /* string */
read write const
a=b<=c:=-d;

 break


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_gga_x *)
(* prefix:
  gga_k_dk_params *params;

  assert(p->params != NULL);
  params = (gga_k_dk_params * )(p->params);
*)

f := x ->
  add(1*params_a_aa[i]*x^(2*(i-1)), i=1..5) /
  add(1*params_a_bb[i]*x^(2*(i-1)), i=1..5):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
TITLE
    AircraftAlloc;

INDEX
    aircraft := (A,B,C,D);

STOCHASTIC

INDEP
    route := (NL1, NL2, ND0, ND1, NB0) -> (r1,r2,r3,r4,r5);

OUTCOME
   out5 := 1..5;

EVENT
    RouteOut[route,out5] :=
       (NL1, 1..5,
        NL2, 1..2,
        ND0, 1..5,
        ND1, 1..5,
        NB0, 1..3);

PROBABILITIES
    p[route, out5 IN RouteOut] :=
[NL1, 1, 0.2,
 NL1, 2, 0.05,
 NL1, 3, 0.35,
 NL1, 4, 0.2,
 NL1, 5, 0.2,
 NL2, 1, 0.3,
 NL2, 2, 0.7,
 ND0, 1, 0.1,
 ND0, 2, 0.2,
 ND0, 3, 0.4,
 ND0, 4, 0.2,
 ND0, 5, 0.1,
 ND1, 1, 0.2,
 ND1, 2, 0.2,
 ND1, 3, 0.3,
 ND1, 4, 0.2,
 ND1, 5, 0.1,
 NB0, 1, 0.1,
 NB0, 2, 0.8,
 NB0, 3, 0.1];

RANDOM DATA
    Demand[route, out5 IN RouteOut] :=
[NL1, 1, 200,
 NL1, 2, 220,
 NL1, 3, 250,
 NL1, 4, 270,
 NL1, 5, 300,
 NL2, 1,  50,
 NL2, 2, 150,
 ND0, 1, 140,
 ND0, 2, 160,
 ND0, 3, 180,
 ND0, 4, 200,
 ND0, 5, 220,
 ND1, 1,  10,
 ND1, 2,  50,
 ND1, 3,  80,
 ND1, 4, 100,
 ND1, 5, 340,
 NB0, 1, 580,
 NB0, 2, 600,
 NB0, 3, 620];

DATA
    AircraftAvail[aircraft] := (10, 19, 25, 15);
    !Demand[route]  := (250, 120, 180, 90, 600);
    PriceTicket[route] := (13, 13, 7, 7, 1);

    PotentialPass[aircraft, route]
       := (16, 15, 28, 23, 81,
             , 10, 14, 15, 57,
             ,  5,   ,  7, 29,
            9, 11, 22, 17, 55);

    MonthlyCost[aircraft, route]
       := (18, 21, 18, 16, 10,
             , 15, 16, 14,  9,
             , 10,   ,  9,  6,
           17, 16, 17, 15, 10);

DECISION VARIABLES
   Assign[aircraft, route]  -> x
      WHERE (PotentialPass > 0);

STAGE2 VARIABLES
   EmptySeats[route] -> y1;
   TurnedAway[route] -> y2;

MACRO
   OperatingCost := SUM(aircraft, route: MonthlyCost * Assign);
   LostRevenue   := SUM(route: PriceTicket * TurnedAway);

MODEL

    MIN TotalCost = OperatingCost + LostRevenue;

SUBJECT TO
   AircraftCap[aircraft]:
       SUM(route: Assign)
    <=
       AircraftAvail;

   PassBal[route]:
       SUM(aircraft: PotentialPass * Assign)
       + TurnedAway[route]
       - EmptySeats[route]
     =
       Demand[route];
END


text:
 func $foo1 (
  var %i i32, var %j i32
  ) i32 {
   return (
     add i32(dread i32 %i, dread i32 %j))}

 func $addi32I (
  var %i i32
  ) i32 {
   return (
     add i32(dread i32 %i,
#constval i32 0xfff))}
#       constval i32 0xffff))}
       constval i32 0xfffff))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: mgga_exc *)
(* prefix:
  mgga_x_tpss_params *params;

  assert(p->params != NULL);
  params = (mgga_x_tpss_params * ) (p->params);
*)

$ifdef mgga_x_revtpss_params
params_a_b      := 0.40:
params_a_c      := 2.35203946:
params_a_e      := 2.16769874:
params_a_kappa  := 0.804:
params_a_mu     := 0.14:
params_a_BLOC_a := 3.0:
params_a_BLOC_b := 0.0:
$endif

tpss_ff     := z -> params_a_BLOC_a + params_a_BLOC_b*z:
tpss_kappa := (x, t) -> params_a_kappa:

$include "tpss_x.mpl"

(* Equation (5) *)

tpss_a1  := (x, t) -> tpss_kappa(x, t)/(tpss_kappa(x, t) + tpss_fx(x, t)):
tpss_f   := (x, u, t) -> 1 + tpss_kappa(x, t)*(1 - tpss_a1(x, t)):

f := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) -> mgga_exchange(tpss_f, rs, z, xs0, xs1, u0, u1, t0, t1):


text:
is_numerical_poly := proc(g,u)
 local m,h;

 if not(type(g,polynom(rational,u))) then
  return false;
 fi;

 m := degree(g,u);
 h := g;
 while m >= 0 do
  if not(type(subs(u=0,h),integer)) then
   return(false);
  fi;
  h := expand(subs(u=u+1,h)-h);
  m := m-1;
 od:
 return(true);
end:

`b/numerical_poly` := (k::nonnegint,u) -> expand(binomial(u,k));

`alpha/numerical_poly` := proc(i::nonnegint,j::nonnegint,k::nonnegint)
 if i <= k and j <= k and k <= i+j then
  return k!/((k-i)!*(k-j)!*(i+j-k)!);
 else
  return 0;
 fi;
end:

`beta/numerical_poly` := proc(i::nonnegint,j::nonnegint,k::nonnegint)
 local m,n;

 add(add(
  (-1)^(m+n)*binomial(i,m)*binomial(j,n)*binomial((i-m)*(j-n),k),
   n=0..j),m=0..i);

end:

`Delta/numerical_poly` := proc(f,u)
 subs(u = u+1,f) - f;
end:

`gamma/numerical_poly` := proc(p,n)
 local N;
 N := floor(n/(p-1));
 N + padic_val(N!,p);
end:

`Gamma/numerical_poly` := proc(n::posint)
 local g,i,p;
 g := 1;
 i := 1;
 p := ithprime(i);
 while p <= 2*n do
  g := g * p^`gamma/numerical_poly`(p,n);
  i := i+1;
  p := ithprime(i);
 od;
 return(g);
end:

`Gamma_alt/numerical_poly` := proc(n::posint)
 local k;

 if modp(n,2) <> 0 then return FAIL; fi;

 return mul(2*denom(bernoulli(2*k)/(2*k)),k=1..n/2);
end:

# There is a canonical basis for the ring of stably numerical
# polynomials, as discussed in the document ImJ.tex.
# The first few basis elements are as follows:

`f/stably_numerical_poly` := table([
  0 = ((u) -> 1),
  1 = ((u) -> (u+1)/2),
  2 = ((u) -> (u^2-1)/24),
  3 = ((u) -> (u^3+u^2-u-1)/48),
  4 = ((u) -> (u^4+70*u^2-71)/5760),
  5 = ((u) -> (u^5+u^4+70*u^3+70*u^2-71*u-71)/11520)
]);

`f_witness/stably_numerical_poly` := table([
 "max_deg" = 5,
 "u_shift" = table([0 = 0,1 = 1,2 = 3,3 = 4,4 = 7,5 = 8]),
 "eval_points"  = [1,-1,7,11,13,17,23,37],
 "split_matrix" =
 [[0, 1, 0, 0, 0, 0, 0, 0],
  [1, -1, 0, 0, 0, 0, 0, 0],
  [-2497, 1127, 5994, -13853, 10734, -1511, 2, 4],
  [-628, 283, 1513, -3505, 2720, -385, 1, 1],
  [-4976, 2245, 11943, -27574, 21348, -2996, 2, 8],
  [-4334, 1955, 10406, -24026, 18599, -2608, 1, 7]]
]):

`find_split_matrix/stably_numerical_poly` := proc()
 local f,W,d,A,m,B,eqs,sol,i,j;

 f := `f/stably_numerical_poly`;
 W := `f_witness/stably_numerical_poly`;
 d := W["max_deg"];

 A := Transpose(Matrix(
  [seq(map(f[i],W["eval_points"]),i=0..d)]));

 m := nops(W["eval_points"]);

 B := Matrix([seq([seq(b[i,j],j=1..m)],i=0..d)]):

 eqs := {op(map(r -> r = 0,map(op,convert(B.A - IdentityMatrix(d+1),listlist))))}:
 sol := isolve(eqs);

 if sol = NULL then return FAIL; fi;

 B := subs(sol,B);
 B := subs(map(u -> u=0,indets(B)),B);

 return B;
end:

######################################################################

`check_witness/stably_numerical_poly` := proc()
 global Z;
 local f,W,d,s,i,j,M,A,B,p;

 f := `f/stably_numerical_poly`;
 W := `f_witness/stably_numerical_poly`;
 d := W["max_deg"];

 _ASSERT(type(d,posint),"max_deg is a positive integer");

 for i from 0 to d do
  s := W["u_shift"][i];
  _ASSERT(is_numerical_poly(u^s*f[i](u),u),
          sprintf("f_%d(u) is stably numerical",i));
 od:

 for i from 0 to d do
  M[i] := 1/coeff(f[i](u),u,i);
  _ASSERT(type(M[i],posint),sprintf("Leading coefficient of f_%d(u) is egyptian",i));
  _ASSERT(
   `and`(seq(-1/(2*M[i]) < coeff(f[j](u),u,i) and
             coeff(f[j](u),u,i) <= 1/(2*M[i]),j=i+1..d)),
   sprintf("Coefficient bounds for f_%d(u)",i));
 od:

 _ASSERT(type(W["eval_points"],list(integer)),"eval_points is a list of integers");

 for p in W["eval_points"] do
  _ASSERT(p = 1 or p = -1 or (p > d and isprime(p)),
          sprintf("evaluation point %d is +/- 1 or prime and greater than max_deg",p));
 od:

 A := Transpose(Matrix(
  [seq(map(f[i],W["eval_points"]),i=0..d)]));

 B := Matrix(W["split_matrix"]);
 Z := [A,B];

 _ASSERT(Equal(B.A,IdentityMatrix(d+1)),"split_matrix splits");

end:

######################################################################

check_numerical_poly := proc()
 local u,v,i,j,k,n,ok,err;

 _ASSERT(
  `b/numerical_poly`(0,u) = 1,
  "b_0(u) = 1"
 );

 _ASSERT(
  `b/numerical_poly`(1,u) = u,
  "b_1(u) = u"
 );

 ok := true;
 for k from 0 to 10 do
  err :=
   `b/numerical_poly`(k,u+v) -
    add(`b/numerical_poly`(i,u) * `b/numerical_poly`(k-i,v),i=0..k);

  err := expand(err);
  if err <> 0 then ok := false; break; fi;
 od:
 _ASSERT(ok,"Formula for b_k(u+v)");

 ok := true;
 for k from 0 to 10 do
  err :=
   `b/numerical_poly`(k,u*v) -
    add(add(
     `beta/numerical_poly`(i,j,k) *
      `b/numerical_poly`(i,u) *
       `b/numerical_poly`(j,v),
    j=0..k),i=0..k);

  err := expand(err);
  if err <> 0 then ok := false; break; fi;
 od:
 _ASSERT(ok,"Formula for b_k(uv)");

 ok := true;
 for i from 0 to 10 do
  for j from 0 to 10 do
   err := `b/numerical_poly`(i,u) * `b/numerical_poly`(j,u) -
    add(
     `alpha/numerical_poly`(i,j,k) * `b/numerical_poly`(k,u),
      k = 0 .. i+j);
   err := expand(err);

   if err <> 0 then ok := false; break; fi;
  od:

  if not(ok) then break; fi;
 od:
 _ASSERT(ok,"Formula for b_i(u) b_j(u)");

 _ASSERT(`Delta/numerical_poly`(`b/numerical_poly`(0,u),u) = 0,"Delta(b_0) = 0");

 _ASSERT(
  `and`(seq(
   expand(`Delta/numerical_poly`(`b/numerical_poly`(i+1,u),u) - `b/numerical_poly`(i,u))=0,
    i=0..10)),
   "Delta(b_{i+1}(u)) = b_i(u)"
 );

 _ASSERT(
  `and`(seq(`Gamma_alt/numerical_poly`(2*n) =
            `Gamma/numerical_poly`(2*n),n=1..10)),
  "Gamma in terms of Bernoulli numbers"
 );
end:


text:
#You can also play with deg
deg:=2:
for n from 4 to 10 do
  for s from 2 to n-1 do
    vars:=[seq(cat(x,i),i=1..n)]:
    f:=expand(add(randpoly(vars,degree=deg)^2,i=1..s)):

    #To decompose
    F:=[f, seq(diff(f,v), v in vars[2..n])]:

#  rr:=rand(1..512):
#  l:=add(rr()*v, v in vars):

#  J:=linalg:-jacobian([f,l], vars):
#  J:=convert(J, Matrix):
#  with(LinearAlgebra):
#  lcols:=combinat:-choose(n, 2):

#  #To decompose
#  F1:=[f, seq(Determinant(SubMatrix(J, [1,2], cols)), cols in lcols)]:

   str:=cat("./soscritpoints_deg_", deg, "_nvars_", n, "_s_",s,".jl");
   fd:=fopen(str, WRITE):
   fprintf(fd, "using SignatureGB\n");
   fprintf(fd, "using Singular\n");
   fprintf(fd, "P, (%a, ", vars[1]);
   for j from 2 to n-1 do
   fprintf(fd, "%a, ", vars[j]);
   od:
   fprintf(fd, "%a)  = Singular.PolynomialRing(Fp(65521), [", vars[n]);
   for j from 1 to n-1 do
   fprintf(fd, "\"%a\", ", vars[j]);
   od:
   fprintf(fd, "\"%a\"]);\n", vars[n]);
   fprintf(fd, "gb = SignatureGB.f5(%a);", vars[1..n]);

   fprintf(fd, "F = %a;\n", F);
   fprintf(fd, "sort!(F, by = p -> leading_monomial(p))\n");
   fprintf(fd, "println(@elapsed SignatureGB.f5(F, verbose=true));");
   fprintf(fd, "exit()\n");
   fclose(fd):
  od:
od:


text:
`is_element/binary_standard_stasheff_trees` := (n::posint) -> proc(TT)
 global reason;

 if not(`is_element/standard_stasheff_trees`(n)(TT)) then
  reason := [convert(procname,string),"TT is not a standard Stasheff tree",TT,reason];
  return false;
 fi;

 if not(`is_binary/standard_stasheff_trees`(n)(TT)) then
  reason := [convert(procname,string),"TT is not binary",TT,reason];
  return false;
 fi;

 return true;
end:

`is_equal/binary_standard_stasheff_trees` :=
 eval(`is_equal/standard_stasheff_trees`);

`is_leq/binary_standard_stasheff_trees` :=
 eval(`is_leq/standard_stasheff_trees`);

`list_elements/binary_standard_stasheff_trees` := proc(n::posint)
 option remember;
 local L,A,UUU,VVV,UU,VV,p,i;

 if n = 0 then
  return [];
 elif n = 1 then
  return [{{1}}];
 fi;

 L := NULL;
 A := {seq(i,i=1..n)};

 for p from 1 to n-1 do
  UUU := `list_elements/binary_standard_stasheff_trees`(p);
  VVV := `list_elements/binary_standard_stasheff_trees`(n-p);
  VVV := map(VV -> map(V -> map(v -> v+p,V),VV),VVV);
  L := L,seq(seq({op(UU),op(VV),A},VV in VVV),UU in UUU);
 od;
 return [L];
end:

`random_element/binary_standard_stasheff_trees` := (n::posint) -> proc()
 local L;

 L := `list_elements/binary_standard_stasheff_trees`(n);
 L[rand(1..nops(L))()];
end:

`count_elements/binary_standard_stasheff_trees` :=
 (n::posint) -> binomial(2*(n-1),n-1)/n;


text:
var $iconst1 <[4] [4] i32> = [ [1007, 707, -273, 75], [0113, 0x4b, 75u, 75l], [75ul, 75lu, 31425926, 60223], [60223, 1619, 30, 314159]]
var $fconst1 <[4] [4] f64> = [ [1007.0, 707.0, -273.0, 75.0], [113.1, 4.0, 75.0, 75.1], [75.0, 75.1, 3.1425926, 6.02e23], [6.02e+23, 1.6e-19, 3.0, 3.14159]]
func &printf (var %p1 <* i8>)void
func &main ( ) i32 {
   call &printf (addrof a32 $fconst1)
   return (constval i32 0) }

 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
# ;, / and | are allowed name characters, but they cannot be the starting char
 func $foo (
  var %i; i32,
  var %j|/ i32
  ) i32 {
   return (
     add i32(dread i32 %i;, dread i32 %j|/))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)

lg93_ad  := 1e-8:
lg93_b   := 0.024974:
lg93_a2  := (lg93_ad + 0.1234)/lg93_b:
lg93_a4  := 29.790:
lg93_a6  := 22.417:
lg93_a8  := 12.119:
lg93_a10 := 1570.1:
lg93_a12 := 55.944:

lg93_f0 := s-> 1 + lg93_a2*s^2 + lg93_a4*s^4
        + lg93_a6*s^6 + lg93_a8*s^8 + lg93_a10*s^10 + lg93_a12*s^12:

lg93_f1 := s-> lg93_f0(s)^lg93_b/(1 + lg93_ad*s^2):

lg93_f  := x->lg93_f1(X2S*x):

f := (rs, zeta, xt, xs0, xs1) -> gga_exchange(lg93_f, rs, zeta, xs0, xs1):


text:
`basis/ass_operad` := (A::set) -> map(u -> bar(op(u)),`list_elements/ord`(A));

`vec/ass_operad` := (A::set) -> proc(u)
 local n,ub,uc,T,B,i;
 global `vec_table/ass_operad`;
 n := nops(A);
 if u = 0 then
  return [0$(n!)];
 elif type(u,`+`) then
  return map(`vec/ass_operad`(A),u);
 else
  if type(u,specfunc(bar)) then
   ub := u;
   uc := 1;
  elif type(u,`*`) then
   ub,uc := selectremove(type,u,specfunc(bar));
   ub := map(op,bar(ub));
  else
   return FAIL;
  fi;
  if not(type(`vec_table/ass_operad`[A],table)) then
   T := table():
   B := `basis/ass_operad`(A);
   for i from 1 to n! do
    T[B[i]] := [0$(i-1),1,0$(n!-i)];
   od:
   `vec_table/ass_operad`[A] := eval(T);
  fi;
  return uc *~ `vec_table/ass_operad`[A][ub];
 fi;
end:


text:
process ACV_loadChunks {
/summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_Chunkloader",Duration:1}
/spreadplayers 5 5 0 1 false @e[type=area_effect_cloud,name=ACV_Chunkloader]
/kill @e[type=area_effect_cloud,name=ACV_Chunkloader]
}


text:
flavor 1
srclang 3
id 65535
numfuncs 2
import "HelloWorld.mplt"
import "/home/bravewtz/Desktop/openarkcompiler/libjava-core/libjava-core.mplt"
entryfunc &LHelloWorld_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V
fileinfo {
  @INFO_filename "HelloWorld.jar"}
srcfileinfo {
  1 "HelloWorld.java"}
javaclass $LHelloWorld_3B <$LHelloWorld_3B> public
func &LHelloWorld_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$LHelloWorld_3B>>) void
func &LHelloWorld_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg2_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void
var $_C_STR_d9bd4a0c2c2117158ed933ab7468a461 fstatic <[4] u64> readonly = [0, 0x1900000000, 0x6c6c6548c63cb61d, 0x21646c726f57206f]
var $Ljava_2Flang_2FSystem_3B_7Cout extern <* <$Ljava_2Fio_2FPrintStream_3B>> final public static
var $__cinf_Ljava_2Flang_2FString_3B <$__class_meta__>
func &MCC_GetOrInsertLiteral () <* <$Ljava_2Flang_2FString_3B>>
func &LHelloWorld_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$LHelloWorld_3B>>) void {
  funcid 48153
  var %Reg1_R43694 <* <$LHelloWorld_3B>>
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>>

  dassign %Reg1_R43694 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43694))
  callassigned &Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  return ()
}
func &LHelloWorld_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg2_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void {
  funcid 48154
  var %Reg0_R562 <* <$Ljava_2Fio_2FPrintStream_3B>>
  var %Reg1_R43 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161334 <* <$Ljava_2Flang_2FString_3B>>

  intrinsiccallwithtype <$LHelloWorld_3B> JAVA_CLINIT_CHECK ()
  #INSTIDX : 0||0000:  getstatic
  intrinsiccallwithtype <$Ljava_2Flang_2FSystem_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R562 0 (dread ref $Ljava_2Flang_2FSystem_3B_7Cout)
  #INSTIDX : 3||0003:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_d9bd4a0c2c2117158ed933ab7468a461) { dassign %L_STR_161334 0 }
  dassign %Reg1_R43 0 (dread ptr %L_STR_161334)
  #INSTIDX : 5||0005:  invokevirtual
  virtualcallassigned &Ljava_2Fio_2FPrintStream_3B_7Cprintln_7C_28Ljava_2Flang_2FString_3B_29V (dread ref %Reg0_R562, dread ref %Reg1_R43) {}
  #INSTIDX : 8||0008:  return
  return ()
}

max_stars_repo_name: Bravewtz/OpenArkCompiler

max_issues_repo_name: Bravewtz/OpenArkCompiler

max_forks_repo_name: Bravewtz/OpenArkCompiler


text:
(*
 Copyright (C) 2017 M.A.L. Marques
               2019 Susi Lehtola

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)
(* prefix:
  gga_x_am05_params *params;

  assert(p->params != NULL);
  params = (gga_x_am05_params * )(p->params);
*)

am05_d     := 28.23705740248932030511071641312341561894: (* POW(CBRT(4/3) * 2*M_PI/3, 4) *)

am05_csi  := s -> (3/2 * LambertW(s^(3/2) / (2*sqrt(6))))^(2/3):
am05_fb   := s -> Pi/3 * s/(am05_csi(s) * (am05_d + am05_csi(s)^2)^(1/4)):
am05_flaa := s -> (1 + params_a_c*s^2)/(1 + params_a_c*s^2/am05_fb(s)):
am05_XX   := s -> 1 - params_a_alpha*s^2/(1 + params_a_alpha*s^2):

am05_f    := x ->  am05_XX(X2S*x) + (1 - am05_XX(X2S*x))*am05_flaa(X2S*x):

f := (rs, zeta, xt, xs0, xs1) -> gga_exchange(am05_f, rs, zeta, xs0, xs1):