text:
n_max := 6; # Number of generators to use in the dual Steenrod algebra

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,integer) and modp(a,2) = 0 then
   return 0;
  elif type(a,`*`) then
   ac,av := selectremove(type,a,rational);
   ac := modp(ac,2);
   if ac <> 1 then
    return ac * bar(u,av,v);
   fi;
  elif type(a,specfunc(bar)) then
   return bar(u,op(a),v);
  fi;
 od;
 return 'procname'(args);
end:

steenrod_length := proc(u::specfunc(nonnegint,Sq)) nops(u); end:
steenrod_deg    := proc(u::specfunc(nonnegint,Sq)) `+`(op(u)); end:
steenrod_excess := proc(u::specfunc(nonnegint,Sq))
 add(op(i,u) - 2*op(i+1,u), i=1..nops(u)-1);
end:

is_admissible := proc(u::specfunc(nonnegint,Sq))
 local m,i;
 m := nops(u);
 for i from 1 to m-1 do
  if op(i,u) < 2*op(i+1,u) then return false; fi;
 od;
 return true;
end:

admissible_basis := proc(n::nonnegint)
 option remember;
 local L,i,v;

 if n = 0 then return [Sq()]; fi;
 L := NULL;
 for i from 1 to n do
  for v in admissible_basis(n-i) do
   if i = n or i >= 2 * op(1,v) then
    L := L,Sq(i,op(v));
   fi;
  od;
 od;

 return [L];
end:

admissible_lt := proc(u,v)
 local i;
 if steenrod_deg(u) < steenrod_deg(v) then return true;  fi;
 if steenrod_deg(u) > steenrod_deg(v) then return false; fi;

 for i from 1 to min(nops(u),nops(v)) do
  if op(i,u) < op(i,v) then return true; fi;
  if op(i,u) > op(i,v) then return false; fi;
 od:

 if nops(u) > nops(v) then return true; fi;
 return false;
end:

mate := proc(u::specfunc(nonnegint,Sq))
 local m,i,t;

 if not(is_admissible(u)) then return FAIL; fi;

 m := nops(u);
 i := [op(u),0];
 return mul(a[t]^(i[t] - 2 * i[t+1]),t=1..m);
end:

comate := proc(v)
 local r,alpha;
 r := n_max;
 while r > 0 and degree(v,a[r]) = 0 do r := r - 1; od;
 if r = 0 then return Sq(); fi;
 alpha := [seq(degree(v,a[t]),t=1..r)];
 return Sq(seq(add(2^(i-t)*alpha[i],i=t..r),t=1..r));
end;

adem_rhs := proc(k::nonnegint,j::nonnegint)
 modp(add(binomial(j-m-1,k-2*m) * Sq(j+k-m,m),m=0..floor(k/2)),2);
end:

adem_relation := (k,j) -> modp(Sq(k,j) + adem_rhs(k,j),2);

adem_reduce0 := proc(u)
 local v,v0,v1,m,r,i;
 if type(u,list) or type (u,set) or type(u,specfunc(anything,bar)) then
  return(map(adem_reduce0,u));
 elif type (u,`+`) then
  return modp(expand(map(adem_reduce0,u)),2);
 elif type (u,`*`) then
  return modp(expand(map(adem_reduce0,u)),2);
 elif type(u,specfunc(nonnegint,Sq)) then
  v := remove(i -> i=0, [op(u)]);
  m := nops(v);
  if m < 2 then return Sq(op(v)); fi;
  i := 1;
  while i < m and v[i] >= 2 * v[i+1] do
   i := i + 1;
  od;
  if i = m then return Sq(op(v)); fi;
  r := adem_rhs(v[i],v[i+1]);
  if r = 0 then return 0; fi;
  if type(r,`+`) then r := [op(r)] else r := [r]; fi;
  v0 := op(1..i-1,v);
  v1 := op(i+2..-1,v);
  r := map(w -> Sq(v0,op(w),v1),r);
  return modp(`+`(op(r)),2);
 else
  return modp(u,2);
 fi;
end:

adem_reduce := proc(u)
 local v,w;
 v := u;
 w := adem_reduce0(u);
 while w <> v do
  v := w;
  w := adem_reduce0(v);
 od;
 return w;
end:

bullett_macdonald_series := (n,s,t) ->
  modp(expand(add(add(
    ((t^2+s*t)^i*s^(2*j) + (s^2+s*t)^i*t^(2*j)) * Sq(i,j),
      j=0..n-i),i=0..n)),2);

bullett_macdonald_term := proc(i::nonnegint,j::nonnegint)
 local k0;
 k0 := modp(i,2);
 if modp(j,2) <> k0 then return 0; fi;
 return add(
   modp(binomial((j+k)/2,k),2) * Sq((j+k)/2,(i-k)/2),
   k=k0..min(i,j),2);
end:

bullett_macdonald_relation := proc(i::nonnegint,j::nonnegint)
 modp(bullett_macdonald_term(i,j) + bullett_macdonald_term(j,i),2);
end:

singer_relation := proc(p::nonnegint,q::nonnegint)
 if p > 2*q then return 0; fi;

 return add(modp(binomial(p,i),2) * Sq(2*q-i+1,q-p+i+1),
            i = max(0,p-q-1)..p);
end:

steenrod_mul := proc()
 local u,v;
 if nargs = 0 then
  return Sq();
 elif nargs = 1 then
  return modp(args[1],2);
 elif nargs > 2 then
  return steenrod_mul(steenrod_mul(args[1],args[2]),args[3..-1]);
 fi;

 u := modp(args[1],2);
 v := modp(args[2],2);
 if type(v,`+`) then
  return modp(expand(map(t -> steenrod_mul(u,t),v)),2);
 fi;
 if type(u,`+`) then
  return modp(expand(map(t -> steenrod_mul(t,v),u)),2);
 fi;

 if type(u,specfunc(nonnegint,Sq)) and type(v,specfunc(nonnegint,Sq)) then
  return Sq(op(u),op(v));
 else
  return u * v;
 fi;
end:

adem_mul := proc() adem_reduce(steenrod_mul(args)); end:

two_power_Sq := proc(n::nonnegint)
 option remember;
 local m,p,s,t,rel;

 if n = 0 then return Sq(); fi;

 m := n;
 p := 1;
 while modp(m,2) = 0 do m := m/2; p := 2*p; od;
 if p = n then return Sq(n); fi;

 rel := bullett_macdonald_series(n,s,t);
 rel := coeff(coeff(rel,s,p),t,2*n-p) - Sq(n,0);

 return two_power_reduce(rel);
end:

two_power_reduce := proc(u)
 if type(u,list) or type(u,set)  or type(u,specfunc(anything,bar)) then
  return map(two_power_reduce,u);
 elif type(u,`+`) or type(u,`*`) then
  return modp(expand(map(two_power_reduce,u)),2);
 elif type(u,specfunc(nonnegint,Sq)) then
  return steenrod_mul(op(map(two_power_Sq, [op(u)])));
 fi;
end:

milnor_Q := proc(n::nonnegint)
 option remember;
 local u,v;

 if n = 0 then
  return Sq(1);
 else
  u := Sq(2^n);
  v := milnor_Q(n-1);
  return modp(adem_mul(u,v) + adem_mul(v,u),2);
 fi;
end:

cartan_psi := proc(u)
 local v,f,m,i,j,k;
 if type(u,list) or type(u,set) then
  return map(cartan_psi,u);
 elif type(u,`+`) or type(u,`*`) then
  return modp(expand(map(cartan_psi,u)),2);
 elif type(u,specfunc(nonnegint,Sq)) then
  v := [bar(Sq(),Sq())];
  f := (i,j) -> v -> bar(Sq(op(op(1,v)),i),Sq(op(op(2,v)),j));
  for k from 1 to nops(u) do
   m := op(k,u);
   v := [seq(seq(f(i,m-i)(x),x in v),i=0..m)];
  od:
  return `+`(op(v));
 else
  return u;
 fi;
end:

a[0] := 1;
aL[0] := 1;
aR[0] := 1;

dual_psi := (u) ->
 subs({seq(a[i] = add(aL[j] * aR[i-j]^(2^j),j=0..i),i=0..n_max)},u);

steenrod_pairing := proc(u,v)
 local m,c,i;

 if u = 0 or v = 0 then return 0; fi;

 if type(u,`+`) then
  return modp(expand(map(t -> steenrod_pairing(t,v),u)),2);
 fi;
 if type(v,`+`) then
  return modp(expand(map(t -> steenrod_pairing(u,t),v)),2);
 fi;

 if type(u,specfunc(nonnegint,Sq)) then
  m := nops(u);
  if m = 0 then
   c := v;
   for i from 1 to n_max do c := coeff(c,a[i],0); od;
   return modp(c,2);
  elif m = 1 then
   c := coeff(v,a[1],op(u));
   for i from 2 to n_max do c := coeff(c,a[i],0); od;
   return modp(c,2);
  else
   c := coeff(dual_psi(v),aR[1],op(1,u));
   for i from 2 to n_max do c := coeff(c,aR[i],0); od;
   c := eval(subs(aL = a,c));
   c := steenrod_pairing(Sq(op(2..-1,u)),c);
   return c;
  fi;
 else
  return FAIL;
 fi;
end:

dickson_fW := proc(n)
 local t;
 return unapply(add(W[n,i] * t^(2^i),i = 0..n-1) + t^(2^n),t);
end:

dickson_fx := proc(n)
 option remember;
 local t,g,h,i;

 if n = 0 then
  return unapply(t,t);
 else
  g := dickson_fx(n-1)(t);
  h := modp(expand(g * add(coeff(g,t,2^i) * (t^(2^i) + x[n-1]^(2^i)),i=0..n-1)),2);
  return unapply(h,t);
 fi;
end:

dickson_W_deg := (n,i) -> 2^n - 2^i;

dickson_Wx := proc(n,k)
 option remember;
 return coeff(dickson_fx(n)(t),t,2^k);
end:

mui_Vx := (k) -> dickson_fx(k)(x[k]);

mui_V_deg := (i) -> 2^i;

mui_gV := (k) -> unapply(V[k] * t + t^2,t);

mui_gx := (k) -> unapply(mui_Vx(k) * t + t^2,t);

dickson_WV := proc(n,k)
 local L;
 if k > n then return 0; fi;

 L := combinat[choose]([seq(j,j=0..n-1)],n-k):
 return add(mul(V[i[t]]^(2^(t-i[t]+k-1)),t=1..n-k),i in L);
end:

dickson_WV_alt := proc(n,k)
 option remember;

 if k > n or k < 0 then return 0; fi;
 if k = n then return 1; fi;
 if k = 0 then return mul(V[i],i=0..n-1); fi;

 return modp(expand(V[n-1] * dickson_WV(n-1,k) + dickson_WV(n-1,k-1)^2),2);
end:

TSq_x := proc(n) local t; unapply(x[n] + t * x[n]^2, t); end;

TSq_mui_V := proc(n)
 local t,u;
 u := t^(2^n) * V[n]^2 + add(t^(2^n - 2^i) * V[n] * dickson_WV(n,i),i=0..n);
 return unapply(u,t);
end:

TSq_dickson_W := proc(n,k)
 local t,u;
 u := subs(W[n,n] = 1,
       W[n,k]^2 * t^(2^n - 2^k) +
       add(add(W[n,i]*W[n,j]* t^(2^n+2^k-2^i-2^j),j=k+1..n),i=0..k));
 return unapply(u,t);
end:

dickson_w_deg := (n,k) -> 1 - 2^k;

dickson_wW := proc(n::nonnegint,i::nonnegint)
 if i < n then
  return W[n,i]/W[n,0];
 elif i = n then
  return 1/W[n,0];
 else
  return 0;
 fi;
end:

dickson_Ww := proc(n,i)
 if i = 0 then
  return 1/w[n,n];
 elif i <= n then
  return w[n,i]/w[n,n];
 else
  return 0;
 fi;
end:

dickson_wx := (n,i) -> dickson_Wx(n,i)/dickson_Wx(n,0);

dickson_wv := proc(n,k)
 local L;
 L := combinat[choose]([seq(j,j=0..n-1)],k):
 return add(mul(v[i[t]]^(-2^(k-t)),t=1..k),i in L);
end:

dickson_wv_alt := proc(n,k)
 option remember;
 local u;

 if k > n or k < 0 then return 0; fi;
 if k = 0 then return 1; fi;

 u := dickson_wv(n-1,k) + dickson_wv(n-1,k-1)^2/v[n-1];
 u := expand(modp(factor(u),2));
 return u;
end:

dickson_fw := (n) -> unapply(t + add(w[n,i] * t^(2^i),i=1..n),t);

mui_v_deg := (i) -> 1;

mui_vV := (i) -> V[i] / mul(V[j],j=0..i-1);
mui_Vv := (i) -> mul(v[j] ^ (2^max(i-j-1,0)),j=0..i);

mui_vx := (i) -> mui_Vx(i) / mul(mui_Vx(j),j=0..i-1);

mui_gv := (i) -> unapply(t + t^2/v[i],t);

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

check_admissible := proc(m::nonnegint)
 local B,BB,n;

 B := admissible_basis(m);

 _ASSERT(
  `and`(op(map(is_admissible,B))),
  "Admissible moomials are admissible"
 );

 _ASSERT(
  B = sort(B,admissible_lt),
  "Admissible monomials are sorted"
 );

 _ASSERT(
  B = map(comate,map(mate,B)),
  "comate o mate = 1"
 );

 n := nops(B);
 _ASSERT(
  {1,op(map(u -> steenrod_pairing(u,mate(u)),B))} = {1},
  "<A,mate(A)> = 1"
 );

 _ASSERT(
  {0,seq(seq(steenrod_pairing(B[j],mate(B[i])),j=i+1..n),i=1..n)} = {0},
  "<A,mate(B)> = 0 for A > B"
 );

 BB := {seq(seq([u,v],v in B),u in B)};
 BB := select(uv -> steenrod_length(uv[1]) < steenrod_length(uv[2]),BB);
 BB := map(uv -> steenrod_pairing(uv[1],mate(uv[2])), BB);
 _ASSERT(
  BB = {0},
  "<A,mate(B)> = 0 for len(A) < len(B)"
 );

end:

check_adem := proc(m::nonnegint)
 local err;

 err :=
 {0, seq(seq(
      steenrod_pairing(adem_relation(i,m-i),mate(u)),
     i=0..m),u in admissible_basis(m))};

 _ASSERT(err = {0},"Adem relations are compatible with pairing");
end:

check_milnor_Q := proc(m)
 local i,Q;

 for i from 0 to m do Q[i] := milnor_Q(i); od;

 _ASSERT(
  {0,seq(adem_mul(Q[i],Q[i]),i=0..m)} = {0},
  "Q_i^2 = 0"
 );

 _ASSERT(
  {0,seq(seq(modp(adem_mul(Q[i],Q[j])+adem_mul(Q[j],Q[i]),2),j=i+1..m),i=0..m)} = {0},
  "Q_i Q_j = Q_j Q_i"
 );

 _ASSERT(
  {0,seq(modp(eval(adem_reduce(cartan_psi(Q[i]) + bar(Sq(),Q[i]) + bar(Q[i],Sq()))),2),i=0..m)} = {0},
  "Q_i is primitive"
 );
end:

check_two_power := proc(m)
 local B,err;

 B := [seq(seq(seq(Sq(i,j,k),k=0..m),j=0..m),i=0..m)];

 err := {0,op(map(u -> adem_reduce(u + two_power_reduce(u)),B))};

 _ASSERT(err = {0},"Two-power reduction is consistent");
end:

check_bullett_macdonald := proc(m)
 local s,t;

 _ASSERT(
  adem_reduce(bullett_macdonald_series(m,s,t)) = 0,
  "Bullett-Macdonald series"
 );

 _ASSERT(
  {0,seq(seq(adem_reduce(bullett_macdonald_relation(i,j)),j=0..20),i=0..20)} = {0},
  "Bullett-Macdonald relations"
 );

 _ASSERT(
  {0,seq(seq(adem_reduce(singer_relation(p,q)),p=0..2*q),q=0..15)} = {0},
  "Singer relations"
 );
end:

check_dickson := proc(m::nonnegint)
 local t,err;

 _ASSERT(
  modp(expand(dickson_fx(m)(t) - add(dickson_Wx(m,i) * t^(2^i),i=0..m)),2) = 0,
  "The Dickson polynomial is additive"
 );

 _ASSERT(
  modp(expand(dickson_Wx(m,0) + mul(mui_Vx(i),i=0..m-1)),2) = 0,
  "W_{m0} = prod_{i<m} V_i"
 );

 _ASSERT(
  {0,seq(modp(expand(
      dickson_Wx(m,k) + mui_Vx(m-1) * dickson_Wx(m-1,k) + dickson_Wx(m-1,k-1)^2
     ),2),k=1..m-1)} = {0},
  "W_{mk} = V_{m-1} W_{m-1,k} + W_{m-1,k-1}^2"
 );

 _ASSERT(
  modp(expand(dickson_fx(m)(t) + mui_gx(m-1)(dickson_fx(m-1)(t))),2) = 0,
  "f_m(t) = g_{m-1}(f_{m-1}(t))"
 );

end:

check_rules := proc(m)
 local Wx_rule, wx_rule, WV_rule, wv_rule,
       Ww_rule, wW_rule, Vx_rule, vx_rule, Vv_rule, vV_rule;

 Wx_rule := {seq(W[m,i] = dickson_Wx(m,i),i=0..m+1)}:
 wx_rule := {seq(w[m,i] = dickson_wx(m,i),i=0..m+1)}:
 WV_rule := {seq(W[m,i] = dickson_WV(m,i),i=0..m  )}:
 wv_rule := {seq(w[m,i] = dickson_wv(m,i),i=0..m+1)}:
 Ww_rule := {seq(W[m,i] = dickson_Ww(m,i),i=0..m+1)}:
 wW_rule := {seq(w[m,i] = dickson_wW(m,i),i=0..m+1)}:
 Vx_rule := {seq(V[i] = mui_Vx(i),i=0..m+1)}:
 vx_rule := {seq(v[i] = mui_vx(i),i=0..m+1)}:
 Vv_rule := {seq(V[i] = mui_Vv(i),i=0..m+1)}:
 vV_rule := {seq(v[i] = mui_vV(i),i=0..m+1)}:

 _ASSERT(
  {seq(subs(Vv_rule,subs(vV_rule,v[i])) - v[i],i=0..m-1),
   seq(subs(vV_rule,subs(Vv_rule,V[i])) - V[i],i=0..m-1),
   seq(subs(Ww_rule,subs(wW_rule,w[m,i])) - w[m,i],i=1..m),
   seq(subs(wW_rule,subs(Ww_rule,W[m,i])) - W[m,i],i=0..m-1),
   seq(subs(Ww_rule,subs(wW_rule,w[m,i])) - w[m,i],i=1..m),
   seq(subs(vx_rule,subs(Vv_rule,V[i])) - subs(Vx_rule,V[i]),i=1..m),
   seq(subs(wx_rule,subs(Ww_rule,W[m,i])) - subs(Wx_rule,W[m,i]),i=0..m-1),
   seq(modp(factor(subs(vx_rule,subs(wv_rule,w[m,i])) - subs(wx_rule,w[m,i])),2),i=1..m),
   seq(modp(factor(subs(Vx_rule,subs(WV_rule,W[m,i])) - subs(Wx_rule,W[m,i])),2),i=0..m-1),
   0} = {0},
  "translation rules"
 );
end:


text:
x:=1:y:=2:z:=3:w:=4:
fun:=proc()
    global x,y,z;
    print(x,y,z,w);
end proc:
Grid:-Launch(fun,imports=[':-x'=x,'y',"z"]);

ext: mpl


text:
# An n-simplex in R^d is represented as a length (n+1) list of
# length d lists of real numbers.

`is_element/simplices` := (d::posint) -> (n::nonnegint) -> proc(a)
 global reason;
 local i,A,indep;

 if not (type(a,list) and nops(a) = n+1) then
  reason := ["is_element/simplices","a is not a list of length n+1",a,n];
  return false;
 fi;

 for i from 1 to n+1 do
  if not `is_element/R`(d)(a[i]) then
   reason := ["is_element/simplices","a[i] is not in R^d",a,i,d];
   return false;
  fi;
 od;

 A := Matrix([seq(a[i] -~ a[n+1],i=1..n)]);
 indep := evalb(NullSpace(Transpose(A)) = {});
 if not indep then
  reason := ["is_element/simplices","a is not affinely independent",a];
  return false;
 fi;

 return true;
end:

`is_equal/simplices` := (d::posint) -> (n::nonnegint) -> proc(a,b)
 return evalb(simplify(a -~ b) = [[0$d]$(n+1)]);
end:

`is_leq/simplices` := NULL;
`list_elements/simplices` := NULL;
`count_elements/simplices` := NULL;

`random_element/simplices` := (d::posint) -> (n::nonnegint) -> proc(m := 100)
 local indep,a,A,i,j;

 if n > d then
  error "n cannot be bigger than d";
 fi;

 indep := false;
 while not(indep) do
  a := [seq([seq(rand(-m..m)(),j=1..d)],i = 1..n+1)];
  A := Matrix([seq(a[i] -~ a[n+1],i=1..n)]);
  indep := evalb(NullSpace(Transpose(A)) = {});
 od:

 return a;
end:

`barycentre/simplex` := (d::posint)  -> proc(a)
 local n,u,i;

 n := nops(a) - 1;
 u := a[n+1];
 for i from 1 to n do u := u +~ a[i]; od;
 u := u /~ (n+1);
 return u;
end:

`radius/simplex` := (d::posint)  -> proc(a)
 local u,r,i;

 u := `barycentre/simplex`(d)(a);
 r := max(map(v -> add((v[i]-u[i])^2,i=1..d),a));
 return sqrt(r);
end:

# The primal matrix has the vertex vectors as its columns,
# but with an extra row of 1's at the bottom.

`primal_matrix/simplex` := (d::posint)  -> proc(a)
 local n,u,i;

 n := nops(a) - 1;
 return Transpose(Matrix([seq([op(a[i]),1],i=1..n+1)]));
end:

# For a d-simplex in R^d, the dual matrix is the inverse
# of the primal matrix.

`dual_matrix/simplex` := (d::posint)  -> proc(a)
 if nops(a) <> d+1 then return FAIL; fi;

 return 1/`primal_matrix/simplex`(d)(a);
end:

# With two arguments, this returns true iff the intersection
# of the simplices spanned by a and b is the simplex spanned
# by some subset c of a intersect b.  With a third argument c_,
# we require in addition that c = c_.

`intersect_nicely/simplices` := (d::posint) -> proc(a,b,c_)
 local c,aa,bb,cc,ab,n,m,p,q,e,W_cand,W_new,W_inc,W,WA,WB,A,B,r,C,w,z,i,k;

 c := {op(a)} intersect {op(b)};
 if nargs > 2 and c <> {op(c_)} then
  return false;
 fi;

 aa := [seq([op(z),1,1],z in {op(a)} minus c)];
 bb := [seq([op(-~z),-1,1],z in {op(b)} minus c)];
 cc := [seq([op(z),1,0],z in c)];
 ab := [op(aa),op(bb)];
 n := nops(aa);
 m := nops(bb);
 p := nops(cc);
 e := Vector([0$d,0,1]);
 W_cand := combinat[choose](n+m,min(d+2-p,n+m));
 W_cand := {seq({op(W)},W in W_cand)};
 W_cand := select(
  W -> W <> {} and min(op(W)) <= n and max(op(W)) > n,
  W_cand
 );

 while W_cand <> {} do
  W_new := NULL;
  W_inc := NULL;
  for W in W_cand do
   q := nops(W);
   A := Transpose(Matrix([seq(ab[w],w in W),
                          op(cc)]));
   B := ReducedRowEchelonForm(<A|e>);
   if Equal(SubMatrix(B,1..d+2,1..p+q),
            <IdentityMatrix(p+q),Matrix(d+2-p-q,p+q)>) then
    r := [seq(B[i,p+q+1],i=1..p+q)];
    if (p+q=d+2 or B[p+q+1,p+q+1] = 0) then
     if min(seq(r[k],k=1..q)) >= 0 then
      return false;
     fi;
    else
     W_inc := W_inc,W;
    fi;
   else
    WA,WB := selectremove(w -> w <= n,W);
    if nops(WA) > 1 then
     W_new := W_new,seq({op(W)} minus {w},w in WA);
    fi;
    if nops(WB) > 1 then
     W_new := W_new,seq({op(W)} minus {w},w in WB);
    fi;
   fi;
  od;
  W_cand := {W_new};
  for W in W_inc do
   W_cand := remove(U -> U minus W = {},W_cand);
  od:
 od;
 return true;
end:


text:
ctfm_expr :=  exp(x) - exp(lv) + 1/(1-x/rv) - 1/(1-lv/rv);
ctfmp_expr := ((diff(ctfm_expr,x)));
octfmp_expr := -(diff(1/ctfmp_expr,x));

ss:=eval(diff(ctfm_expr,x),{lv=-12,rv=3});
evalf(eval(ss,x=2))/evalf(eval(ss,x=-12));
tval_ctfm:=eval((ctfm_expr),{lv=-12,rv=3});
(evalf(eval(tval_ctfm,x=2)));
read ("/home/arman/FD/FD.mpl"):
MFD(); grid_functions:={};
Gen_Eval_Code(ctfm_expr,input="c",proc_name="init_ctfm");
Gen_Eval_Code(ctfmp_expr,input="c",proc_name="init_ctfmp");
Gen_Eval_Code(octfmp_expr,input="c",proc_name="init_octfmp");

grid_functions:={};
ctfm_expr :=  x/(1-x/rv);
ctfmp_expr := simplify(diff(ctfm_expr,x));
octfmp_expr := -simplify(diff(1/ctfmp_expr,x));
Gen_Eval_Code(ctfm_expr,input="c",proc_name="init_ctfm2");
Gen_Eval_Code(ctfmp_expr,input="c",proc_name="init_ctfmp2");
Gen_Eval_Code(octfmp_expr,input="c",proc_name="init_octfmp2");


text:
# Berechnung und Projektion der Dynamikgleichungen in Regressorform
# Einleitung
# Berechnung und Projektion der Dynamikgleichungen in Regressorform
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# para -> Berechnung für einen parallelen Roboter
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# projection -> Die Dynamikgleichungen werden auf EE-Koordinaten projiziert
# dynamics -> Berechnung der Dynamik
# regressor -> Regressorform (parameterlinear)
#
# TODO
# Die Regressorform ist noch nicht in Minimaldarstellung. Außerdem treten Dynamikparameter auf, die bei reduzierten EE-FG (z.B. 2T1R-Bewegung) keinen Einfluss haben.
#
# Autor
# Tim Job (Studienarbeit bei Moritz Schappler), 2018-12
# Moritz Schappler, moritz.schappler@imes.uni-hannover.de
# (C) Institut für Mechatronische Systeme, Universität Hannover
# Sources
# [GautierKhalil1990] Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots
# [KhalilDombre2002] Modeling, Identification and Control of Robots
# [Ortmaier2014] Vorlesungsskript Robotik I
# [Job2018_S759] Job, T. (Studienarbeit; Betreuer Moritz Schappler): Implementierung einer strukturunabhängigen Dynamikmodellierung für parallelkinematische Maschinen (2018)
# Initialization
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
interface(rtablesize=100):
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
with(StringTools): # Für Zeitausgabe
;
# Einstellungen für Code-Export: Optimierungsgrad (2=höchster) und Aktivierung jedes Terms.
#codegen_act := true: # noch nicht implementiert
codegen_opt := 2:
codeexport_invdyn := true:
codeexport_actcoord := false: # Generierung der Dynamik in Antriebskoordinaten nicht standardmäßig (hoher Rechenaufwand)
;
read "../helper/proc_MatlabExport":
# Roboter-Definitionen laden
read "../robot_codegen_definitions/robot_env_par":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", leg_name):
read "../robot_codegen_definitions/robot_env_par":
# Ergebnisse der Plattform-Dynamik in Regressorform laden

dynamicsfile := sprintf("../codeexport/%s/tmp/floatb_%s_platform_dynamic_reg_maple.m", robot_name, base_method_name):
if FileTools[Exists](dynamicsfile) then
  read dynamicsfile:
else
  printf("%s. Plattform-Dynamik (Regressor) konnte nicht geladen werden. Abbruch der Berechnung.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  quit: # Funktioniert in GUI nicht richtig...
  robot_name := "": # ...Daher auch Löschung des Roboternamens.
end if:
# Neu-Definieren, damit Variablen im Workspace auftauchen
paramVecP := paramVecP:
paramVecP_M := paramVecP_M:
A_E := A_E:
H := H:
dH := dH:
M_regmin := M_regmin:
c_regmin := c_regmin:
g_regmin := g_regmin:
# Ergebnisse der Dynamik der Gelenkkette in Regressorform laden

dynamicsfile := sprintf("../codeexport/%s/tmp/invdyn_%s_%s_maple.m", leg_name, "fixb", "regressor_minpar"):
if FileTools[Exists](dynamicsfile) then
  read dynamicsfile:
else
  printf("%s. PKM-Dynamik konnte nicht geladen werden. Abbruch der Berechnung.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  quit: # Funktioniert in GUI nicht richtig...
  robot_name := "": # ...Daher auch Löschung des Roboternamens.
end if:
tau_regressor_s := tau_regressor_s:
tauMM_regressor_s := tauMM_regressor_s:
MMjj_regressor_s := MMjj_regressor_s:
tauC_regressor_s := tauC_regressor_s:
taug_regressor_s := taug_regressor_s:
read "../robot_codegen_definitions/robot_env_par":

# Ergebnisse des Minimalparametervektors der Gelenkkette laden
read sprintf("../codeexport/%s/tmp/minimal_parameter_vector_fixb_maple", leg_name):
Paramvec2 := Paramvec2:
read "../robot_codegen_definitions/robot_env_par": # Nochmal laden, um Standard-Einstellungen überschreiben zu können.
;
# Ergebnisse der zusätzlichen Definitionen für parallele Roboter laden
read sprintf("../codeexport/%s/tmp/para_definitions", robot_name):
printf("Generiere Parameterlineare Form der Dynamik für PKM %s mit Minimal-Parametersatz\n", robot_name):
# Lade "robotics_repo_path"-File mit Link zum "imes-robotics-matlab"-Repo
read("../robotics_repo_path"):
# Lade die Funktionen aus dem "imes-robotics-matlab"-Repo
read(sprintf("%s/transformation/maple/proc_eul%s2r", robotics_repo_path, angleConvLeg)):
read(sprintf("%s/transformation/maple/proc_eul%sjac", robotics_repo_path, "zyx")):

# Berechne Dynamik-Matrizen für alle Beine
# Alle Basisgeschwindigkeiten und -winkel aus Berechnung der seriellen Kette zu null setzen.
omegaxs_base := 0:
omegays_base := 0:
omegazs_base := 0:
alphaxs_base := 0:
betays_base := 0:
gammazs_base := 0:
vxs_base := 0:
vys_base := 0:
vzs_base := 0:
alphaDx_base := 0:
betaDy_base := 0:
gammaDz_base := 0:
vDxs_base := 0:
vDys_base := 0:
vDzs_base := 0:
alphaDDx_base := 0:
betaDDy_base := 0:
gammaDDz_base := 0:

# Ergebnisse der Kinematik für parallelen Roboter laden

kinematicsfile := sprintf("../codeexport/%s/tmp/kinematics_%s_platform_maple.m", robot_name, base_method_name):
if FileTools[Exists](kinematicsfile) then
  read kinematicsfile:
else
  printf("%s. PKM-Kinematik konnte nicht geladen werden. Abbruch der Berechnung.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  quit: # Funktioniert in GUI nicht richtig...
  robot_name := "": # ...Daher auch Löschung des Roboternamens.
end if:
# Variablen neu laden
pivotMat := pivotMat:
pivotMatMas := pivotMatMas:
P_i := P_i:
Jinv := Jinv:
JB_i := JB_i:
JBD_i := JBD_i:
JBinv_i := JBinv_i:
JBDinv_i := JBDinv_i:
U_i := U_i:
UD_i := UD_i:

# Dynamik-Parameter für virtuelle Segmente nach den Plattform-Koppelgelenken entfernen
NQ := NQ - (NQJ-NQJ_parallel):
Paramvec3:=Paramvec2: # Dynamik-Minimalparametervektor der Beinkette
;
for i from NQJ_parallel+1 to NQJ do
  Paramvec3 := subs({XXC||i = 0}, Paramvec3):
  Paramvec3 := subs({XYC||i = 0}, Paramvec3):
  Paramvec3 := subs({XZC||i = 0}, Paramvec3):
  Paramvec3 := subs({YYC||i = 0}, Paramvec3):
  Paramvec3 := subs({YZC||i = 0}, Paramvec3):
  Paramvec3 := subs({ZZC||i = 0}, Paramvec3):
  Paramvec3 := subs({XX||i  = 0}, Paramvec3):
  Paramvec3 := subs({XY||i  = 0}, Paramvec3):
  Paramvec3 := subs({XZ||i  = 0}, Paramvec3):
  Paramvec3 := subs({YY||i  = 0}, Paramvec3):
  Paramvec3 := subs({YZ||i  = 0}, Paramvec3):
  Paramvec3 := subs({ZZ||i  = 0}, Paramvec3):
  Paramvec3 := subs({SX||i  = 0}, Paramvec3):
  Paramvec3 := subs({SY||i  = 0}, Paramvec3):
  Paramvec3 := subs({SZ||i  = 0}, Paramvec3):
  Paramvec3 := subs({MX||i  = 0}, Paramvec3):
  Paramvec3 := subs({MY||i  = 0}, Paramvec3):
  Paramvec3 := subs({MZ||i  = 0}, Paramvec3):
  Paramvec3 := subs({M||i   = 0}, Paramvec3):
end do:
Paramvec2: # Bein-Dynamik-Parametervektor vor Vereinfachungen
;
Paramvec3: # Bein-Dynamik-Parametervektor nach Vereinfachungen
;
paramVecP: # Plattform-Dynamik-Parametervektor ohne Zusammenfassungen
;
paramVecP_M: # Plattform-Dynamik-Parametervektor mit Zusammenfassung mit Bein-Dynamikparametern
;
# Alle Einträge des Minimalparametervektors, die jetzt Null sind, entfernen
ParamvecNew := Matrix(RowDimension(Paramvec3), 1): # Initialisierung mit maximaler Größe
k := 0:
for j from 1 to RowDimension(Paramvec3) do
  if not Paramvec3(j,1) = 0 then
    k := k + 1:
    ParamvecNew(k,1) := Paramvec3(j,1):
  end if:
end do:
# Neuen Parametervektor für die Beinsegmenten (ohne Plattform)
ParamvecNew := ParamvecNew(1..k,1):
# Dupliziere alle berechneten Matrizen. i steht für den Index des jeweiligen Beines
tauReg := tau_regressor_s(7..NQ,..):
#MatrixEnd := index_symmat2vec(NQJ_parallel,NQJ_parallel,NQJ_parallel):
#MMjjReg := MMjj_regressor_s(1..MatrixEnd,..):
tauMMReg := tauMM_regressor_s(7..NQ,..):
for i to NQ do
  tauMMReg := subs({qD_s(i, 1) = 0}, tauMMReg):
end do:
tauCReg := tauC_regressor_s(7..NQ,..):
taugReg := taug_regressor_s(7..NQ,..):

# Regressor aus Berechnung für serielle Beinketten. Dort erste sechs Einträge für Basis
g := <g1;g2;g3>:
tmp := <tmp1;tmp2;tmp3>:
Rmat := Transpose(parse(sprintf("eul%s2r",angleConvLeg))(frame_A_i(1..3,1))):
gtmp := Rmat.g:
for i to NQJ_parallel do
  for k to ColumnDimension(tauReg) do
    for j to 3 do
      tauReg(i,k) := subs({g(j)=tmp(j)},tauReg(i,k)):
      taugReg(i,k) := subs({g(j)=tmp(j)},taugReg(i,k)):
    end do:
    for j to 3 do
      tauReg(i,k) := subs({tmp(j)=gtmp(j)},tauReg(i,k)):
      taugReg(i,k) := subs({tmp(j)=gtmp(j)},taugReg(i,k)):
    end do:
  end do:
end do:

for i to N_LEGS do
  tau_regressor_s||i := Copy(tauReg):
  tauMM_regressor_s||i := Copy(tauMMReg):
  #MMjj_regressor_s||i := Copy(MMjjReg):
  tauC_regressor_s||i := Copy(tauCReg):
  taug_regressor_s||i := Copy(taugReg):
end do:
COLUMNreg := ColumnDimension(tau_regressor_s):
#
# Substituiere in jeder Matrix den Winkel Alpha (Verdrehung in der Basis) und die Gelenkkoordinaten und -geschwindigkeiten
for k from 1 by 1 to N_LEGS do
  for i to NQJ_parallel do
    for j to COLUMNreg do
      for p from NQJ_parallel+1 to NQJ do
        tau_regressor_s||k(i,j) := subs({qJDD||p||s=0,qJD||p||s=0,qJ||p||s=0},tau_regressor_s||k(i,j)):
        tauMM_regressor_s||k(i,j) := subs({qJDD||p||s=0,qJD||p||s=0,qJ||p||s=0},tauMM_regressor_s||k(i,j)):
        tauC_regressor_s||k(i,j) := subs({qJDD||p||s=0,qJD||p||s=0,qJ||p||s=0},tauC_regressor_s||k(i,j)):
        taug_regressor_s||k(i,j) := subs({qJDD||p||s=0,qJD||p||s=0,qJ||p||s=0},taug_regressor_s||k(i,j)):
      end do:
      for l to 3 do
        tau_regressor_s||k(i,j) := subs({frame_A_i(l,1)=frame_A_i(l,k)},tau_regressor_s||k(i,j)):
        tauMM_regressor_s||k(i,j) := subs({frame_A_i(l,1)=frame_A_i(l,k)},tauMM_regressor_s||k(i,j)):
        tauC_regressor_s||k(i,j) := subs({frame_A_i(l,1)=frame_A_i(l,k)},tauC_regressor_s||k(i,j)):
        taug_regressor_s||k(i,j) := subs({frame_A_i(l,1)=frame_A_i(l,k)},taug_regressor_s||k(i,j)):
      end do:
      for m to NQJ_parallel do
        n := (m + (k-1)*NQJ_parallel):
        tau_regressor_s||k(i,j) := subs({qJDD||m||s=qJDD||n||s,qJD||m||s=qJD||n||s,qJ||m||s=qJ||n||s},tau_regressor_s||k(i,j)):
        tauMM_regressor_s||k(i,j) := subs({qJDD||m||s=qJDD||n||s,qJD||m||s=qJD||n||s,qJ||m||s=qJ||n||s},tauMM_regressor_s||k(i,j)):
        tauC_regressor_s||k(i,j) := subs({qJDD||m||s=qJDD||n||s,qJD||m||s=qJD||n||s,qJ||m||s=qJ||n||s},tauC_regressor_s||k(i,j)):
        taug_regressor_s||k(i,j) := subs({qJDD||m||s=qJDD||n||s,qJD||m||s=qJD||n||s,qJ||m||s=qJ||n||s},taug_regressor_s||k(i,j)):
      end do:
    end do:
  end do:
end do:

printf("%s. Dynamik-Matrizen bestimmt. Beginne Projektion/Addition.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
# Berechnung, Projektion und Addition der Dynamikgleichungen
# Berechnung der Kräfte/Momente an den Gelenken der jeweiligen Beine und Projektion auf EE-Plattform
A_E := A_E:
U_i := U_i:
JBinv_i := JBinv_i:
for i to N_LEGS do
   A_||i:= simplify(Transpose(U_i(..,..,i)).Transpose(JBinv_i(..,..,i))).simplify(tau_regressor_s||i):
   Mreg_||i:= simplify(Transpose(U_i(..,..,i)).Transpose(JBinv_i(..,..,i))).simplify(tauMM_regressor_s||i):
   Creg_||i:= simplify(Transpose(U_i(..,..,i)).Transpose(JBinv_i(..,..,i))).simplify(tauC_regressor_s||i):
   greg_||i:= simplify(Transpose(U_i(..,..,i)).Transpose(JBinv_i(..,..,i))).simplify(taug_regressor_s||i):
end do:
# Aufsummieren aller Kräfte, projiziert auf EE-Plattform
Tmp := 0:
TmpM := 0:
TmpC := 0:
Tmpg := 0:
for i to N_LEGS do
  Tmp := Tmp + A_||i:
  TmpM := TmpM + Mreg_||i:
  TmpC := TmpC + Creg_||i:
  Tmpg := Tmpg + greg_||i:
end do:
A_j := Tmp:
Mreg_j := TmpM:
Creg_j := TmpC:
greg_j := Tmpg:
# Addiere Inverse Dynamik der Plattform
ROWS := RowDimension(ParamvecNew):
#
# Neuer Parametervektor für Beine und Plattform
paramMin := <ParamvecNew;paramVecP>: # TODO: Hier stand vorher paramVecP_M. Damit war die Regressorform aber nicht für alle Systeme konsistent. Jetzt ist die parameterlineare Form aber nicht mehr so stark zusammengefasst, wie sie es sein könnte.
;
# Regressormatrix zusammenstellen
A := pivotMat.<A_j(1..6,1..ROWS)|A_E>:
Mreg := pivotMat.<Mreg_j(1..6,1..ROWS)|M_regmin>:
Creg := pivotMat.<Creg_j(1..6,1..ROWS)|c_regmin>:
greg := pivotMat.<greg_j(1..6,1..ROWS)|g_regmin>:
NotNullEntries := 0:
for i to RowDimension(paramMin) do
	if not(paramMin(i,1) = 0 or Equal(Column(A,i), ZeroVector[column](RowDimension(A)))) then
		NotNullEntries := NotNullEntries + 1:
	end if:
end do:
NotNullEntries:
paramMinRed := Matrix(NotNullEntries,1):
ARed := Matrix(RowDimension(A),NotNullEntries):
MregRed := Matrix(RowDimension(A),NotNullEntries):
CregRed := Matrix(RowDimension(A),NotNullEntries):
gregRed := Matrix(RowDimension(A),NotNullEntries):
j := 1:
for i to RowDimension(paramMin) do
	if not(paramMin(i,1) = 0 or Equal(Column(A,i), ZeroVector[column](RowDimension(A)))) then
		paramMinRed(j,1) := paramMin(i,1):
		for k to RowDimension(A) do
			ARed(k,j) := A(k,i):
			MregRed(k,j) := Mreg(k,i):
			CregRed(k,j) := Creg(k,i):
			gregRed(k,j) := greg(k,i):
		end do:
		j := j + 1:
	end if:
end do:
printf("%s. Beginne Ersetzung der Geschwindigkeiten und Bestimmung der Plattform-Terme.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
# Replace Joint Velocities
# Substituiere die Gelenkgeschwindigkeiten über H-, Ui- und JBi-Matrix mit EE-Geschwindikeiten
Tmp := 0:
Tmp2 := 0:

for i to N_LEGS do
  Tmp := Multiply(H,xED_s):
  Tmp := Multiply(U_i(..,..,i),Tmp):
  z||i := simplify(Multiply(JBinv_i(..,..,i),Tmp)):
  Tmp2 := JBinv_i(..,..,i).(Multiply(U_i(..,..,i),H).xEDD_s+Multiply(U_i(..,..,i),dH).xED_s+Multiply(UD_i(..,..,i),H).xED_s):
  A||i := simplify(Multiply(JBinv_i(..,..,i),JBD_i(..,..,i))):
  B||i := Multiply(A||i,Matrix(qJD_i_s(..,i))):
  C||i := Tmp2 - B||i:
end do:

RowParamMin := RowDimension(paramMinRed):
MtoC := copy(MregRed):
for i to NX do
  for k to RowParamMin do
    for j to N_LEGS do
      for l to NQJ_parallel do
        ARed(i,k) := subs({qJDD_i_s(l,j)=C||j(l)},ARed(i,k)):
        MregRed(i,k) := subs({qJDD_i_s(l,j)=C||j(l)},MregRed(i,k)):
        MtoC(i,k) := subs({qJDD_i_s(l,j)=C||j(l)},MtoC(i,k)):
      end do:
    end do:
    for j to N_LEGS do
      for l to NQJ_parallel do
        ARed(i,k) := subs({qJD_i_s(l,j)=z||j(l)},ARed(i,k)):
        MregRed(i,k) := subs({qJD_i_s(l,j)=z||j(l)},MregRed(i,k)):
        MtoC(i,k) := subs({qJD_i_s(l,j)=z||j(l)},MtoC(i,k)):
        CregRed(i,k) := subs({qJD_i_s(l,j)=z||j(l)},CregRed(i,k)):
      end do:
    end do:
    for j to N_LEGS do
      for l to 6 do
        MregRed(i,k) := subs({xED_s(l,1)=0},MregRed(i,k)):
        MtoC(i,k) := subs({xEDD_s(l,1)=0},MtoC(i,k)):
      end do:
    end do:
  end do:
end do:
CregRed := CregRed+MtoC:
tauGes := ARed:#.paramMinRed:
# Dynamik-Regressor in Plattform-Koordinaten (oben berechnet)
tau_x := tauGes:
# Dynamik-Regressor in Antriebs-Koordinaten umrechnen. Nur machen, wenn die Jacobi-Matrix einfach genug ist.
if RowDimension(Jinv) < 5 and codeexport_actcoord then
  read sprintf("../codeexport/%s/tmp/jacobian_maple.m", robot_name): # Annahme, dass die Jacobi-Matrix (aus Invertierung) vorher berechnet wurde. Setzt aktuell das normale Dynamik-Skript voraus.
  J:=J: # Neudefinition, damit Variable im Workspace ist.
  tau_qa:=Transpose(J).tauGes:
  Mreg_qa:=Transpose(J).MregRed:
  Creg_qa:=Transpose(J).CregRed:
  greg_qa:=Transpose(J).gregRed:
end if:

# Erzeuge Minimalkoordinatenvektor zur Erstellung der Massenmatrix durch Ableitungen im folgenden Schritt
xEDD_s_dummy := Matrix(NX,1):
zaehler := 1:
for i to 6 do
	if not(xEDD_s(i) = 0) then
		xEDD_s_dummy(zaehler) := xEDD_s(i):
		zaehler := zaehler + 1:
	end if:
end do:
# Erzeugung der Massenmatrizen
MMregRed := Matrix(NX*NX, RowDimension(paramMinRed)):
MMreg_qa := Matrix(NX*NX, RowDimension(paramMinRed)):
i_rr := 0:
for i to NX do # Zeilenindex der Massenmatrix
  for j to NX do  # Spaltenindex der Massenmatrix
    #if j > i then
    #  next: # rechte obere Seite der symmetrischen Matrix. Keine neue Information. Nicht berechnen oder speichern.
    #end if:
    i_rr := i_rr + 1: # Gehe zeilenweise durch den unteren linken Teil der Massenmatrix (inkl. Diagonale)
    for k to RowDimension(paramMinRed) do # Spaltenindex der Regressormatrix
      MMregRed[i_rr, k] := diff(MregRed[i, k], xEDD_s_dummy[j, 1]):
      MMreg_qa[i_rr, k] := diff(Mreg_qa[i, k], xEDD_s_dummy[j, 1]):
    end do:
  end do:
end do:
# Dynamik in Plattform-Koordinaten (oben berechnet)
tau_x := ARed:
MMreg_x := MMregRed:
Creg_x := CregRed:
greg_x := gregRed:
# PKM-Dynamik-Funktionen mit MPV bereits eingesetzt
# Erzeuge MPV-Vektor
nMDP := RowDimension(paramMinRed):
PV := Matrix(nMDP,1):
for i from 1 to nMDP do
  PV(i,1) := parse(sprintf("MDP%d%", i)):
end do:

tau_x_mdp := tau_x.PV:
MMreg_x_mdp := MMreg_x.PV:
Creg_x_mdp := Creg_x.PV:
greg_x_mdp := greg_x.PV:

tau_qa_mdp := tau_qa.PV:
MMreg_qa_mdp := MMreg_qa.PV:
Creg_qa_mdp := Creg_qa.PV:
greg_qa_mdp := greg_qa.PV:

# Export
# Maple-Export (zur eventuellen späteren Verarbeitung in Maple)

save tau_x,   sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_reg_maple.m",     robot_name):
save MMreg_x, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_MMreg_maple.m",   robot_name):
save Creg_x,  sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_tauCreg_maple.m", robot_name):
save greg_x,  sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_taugreg_maple.m", robot_name):

# Matlab Export: Floating base (TODO: Was heißt das? Hier ist doch kein Floating Base?)
# Berechnung der Basis-Belastung ist für manche Basis-Darstellungen falsch (siehe oben unter Gravitationslast). (TODO: Ist das noch aktuell?)

if codeexport_invdyn then
  printf("%s. Beginne Code-Export Inverse Dynamik (Regressor) in Plattform-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(tau_x, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_reg_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Massenmatrix (MDP eingesetzt) in Plattform-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(MMreg_x, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_MMreg_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Coriolis-Vektor (MDP eingesetzt) in Plattform-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(Creg_x, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_tauCreg_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Gravitations-Vektor (MDP eingesetzt) in Plattform-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(greg_x, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_taugreg_matlab.m", robot_name), codegen_opt):
end if:

if codeexport_invdyn and RowDimension(Jinv) < 5 and codeexport_actcoord then
  printf("%s. Beginne Code-Export Inverse Dynamik (Regressor) in Antriebs-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(tau_qa, sprintf("../codeexport/%s/tmp/invdyn_para_actcoord_reg_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Massenmatrix (Regressor) in Antriebs-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(MMreg_qa, sprintf("../codeexport/%s/tmp/invdyn_para_actcoord_MMreg_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Coriolis-Vektor (Regressor) in Antriebs-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(Creg_qa, sprintf("../codeexport/%s/tmp/invdyn_para_actcoord_tauCreg_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Gravitations-Vektor (Regressor) in Antriebs-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(greg_qa, sprintf("../codeexport/%s/tmp/invdyn_para_actcoord_taugreg_matlab.m", robot_name), codegen_opt):
end if:

# Maple-Export

save paramMinRed, sprintf("../codeexport/%s/tmp/minimal_parameter_parrob_maple.m", robot_name):
save RowParamMin, sprintf("../codeexport/%s/tmp/RowMinPar_parallel_maple.m",       robot_name):

MatlabExport(paramMinRed, sprintf("../codeexport/%s/tmp/minimal_parameter_parrob_matlab.m", robot_name), codegen_opt):
MatlabExport(RowParamMin, sprintf("../codeexport/%s/tmp/RowMinPar_parallel.m", robot_name), 2);

# PKM-Dynamik-Funktionen mit MPV bereits eingesetzt

# Maple-Export

save tau_x_mdp,   sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_reg_mdp_maple.m",     robot_name):
save MMreg_x_mdp, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_MMreg_mdp_maple.m",   robot_name):
save Creg_x_mdp,  sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_tauCreg_mdp_maple.m", robot_name):
save greg_x_mdp,  sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_taugreg_mdp_maple.m", robot_name):

# Matlab Export:
if codeexport_invdyn then
  printf("%s. Beginne Code-Export Inverse Dynamik (MDP eingesetzt) in Plattform-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(tau_x_mdp, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_reg_mdp_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Massenmatrix (MDP eingesetzt) in Plattform-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(MMreg_x_mdp, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_MMreg_mdp_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Coriolis-Vektor (MDP eingesetzt) in Plattform-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(Creg_x_mdp, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_tauCreg_mdp_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Gravitations-Vektor (MDP eingesetzt) in Plattform-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(greg_x_mdp, sprintf("../codeexport/%s/tmp/invdyn_para_plfcoord_taugreg_mdp_matlab.m", robot_name), codegen_opt):
end if:

if codeexport_invdyn and RowDimension(Jinv) < 5 and codeexport_actcoord then
  printf("%s. Beginne Code-Export Inverse Dynamik (MDP eingesetzt) in Antriebs-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(tau_qa_mdp, sprintf("../codeexport/%s/tmp/invdyn_para_actcoord_reg_mdp_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Massenmatrix (MDP eingesetzt) in Antriebs-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(MMreg_qa_mdp, sprintf("../codeexport/%s/tmp/invdyn_para_actcoord_MMreg_mdp_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Coriolis-Vektor (MDP eingesetzt) in Antriebs-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(Creg_qa_mdp, sprintf("../codeexport/%s/tmp/invdyn_para_actcoord_tauCreg_mdp_matlab.m", robot_name), codegen_opt):
  printf("%s. Beginne Code-Export Gravitations-Vektor (MDP eingesetzt) in Antriebs-Koordinaten.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  MatlabExport(greg_qa_mdp, sprintf("../codeexport/%s/tmp/invdyn_para_actcoord_taugreg_mdp_matlab.m", robot_name), codegen_opt):
end if:

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


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 *)
(* prefix:
  gga_x_hjs_params *params;

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

$include "gga_x_hjs.mpl"

hjs2_xi := 1/(exp(20) - 1):
hjs2_fs := s -> -log((exp(-s) + hjs2_xi)/(1 + hjs2_xi)):

hjs_fx := (rs, z, x) -> hjs_f1(rs, z, hjs2_fs(X2S*x)):


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 *)

mkappa := KAPPA_PBE:
mmu    := 0.00361218645365094697:
malpha := 0.52:

f  := x -> 1 + mkappa*(1 - (1 + mmu*x^2/(malpha*mkappa))^(-malpha)):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
$ifndef _CLOSURE_
$define _CLOSURE_

# 寻找变换矩阵的闭包
# 代入的条件只能是等于0的条件，目的在于删去一些系数之后再求闭包
getClosure:=proc(A::Matrix,sol:={})
    local _a,_b,n;
    n:=LinearAlgebra[RowDimension](A);
    _a:=Matrix([seq(a[i],i=1..n)]);
    _a:=subs(sol[],_a);
    _b:=_a.A;
    _b:=map(x->sym2ind(indets(x,specindex(a))),convert(_b,list));
    return map(findClosure,[seq(1..n)],_b);
end proc:

# 提取a[k]的下标
sym2ind:=proc(s::set)
    return map(x->op(1,x),s);
end proc:

# 寻找某个元素所在的最小闭包
findClosure:=proc(k::integer,A::list)
    local aset,bset;
    aset:={};
    bset:=A[k];
    while aset<>bset do
        aset:=bset;
        bset:=`union`(bset,A[convert(bset,list)][]);
    end do;
    return aset;
end proc:

# 获取最小封闭
# 输入参数为封闭的输出结果
# 忽略空集
getMinClosure:=proc(c::list(set))
    local ns,ind;
    ns:=map(x->`if`(x={},infinity,numelems(x)),c);
    ind:=min[index](ns);
    return c[ind];
end proc:

$endif


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

 func $liori64 (
  var %i i64, var %j i64
  ) i32 {
   return (
     lior i64(dread i64 %i, dread i64 %j))}
 # 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:
with(Groebner):
with(Ore_algebra):

find_koszul_homology := proc(rels,vars_,p_)
 local vars,char,A,T,P,ix_P,S,B0,B1,d,i,j,k,u,e_K,e_Z,e_B,de_K,ee_K,dd,K,
  K_basis,K_rank,Z_basis,Z_rank,B_basis,B_rank,BZ_basis,BZ_rank;

 if nargs > 1 then
  vars := vars_;
 else
  vars := [op(indets(rels,symbol))];
 fi;

 if nargs > 2 then
  char := (characteristic = p_);
 else
  char := NULL;
 fi;

 A := poly_algebra(op(vars),e_K,e_Z,e_B,char);
 T := MonomialOrder(A,lexdeg([e_K,e_B],vars,[e_Z]),{e_K,e_Z,e_B});

 d := nops(rels);
 # P = list of all subsets of 1,..,d
 P := combinat[powerset]([seq(i,i=1..d)]);

 # Now build index for P e.g. if {1,4,6} appears as 10th
 # element of P then if_P[{1,4,6}] should be 10.
 ix_P := table():
 for i from 1 to 2^d do ix_P[P[i]] := i; od:

 # de_K[i] will be the differential applied to the i'th basis
 # element of the exterior algebra (for 1 <= i <= 2^d).
 de_K := table():
 for i from 1 to 2^d do
  S := P[i];
  de_K[i] := add((-1)^(j-1) * e_K^ix_P[subsop(j=NULL,S)] * rels[S[j]],j=1..nops(S));
  ee_K[i] := e_K^i + e_K^(2^d+1) * de_K[i];
 od:

 # dd = the differential = linear extension of the map
 # sending e_K^i to de_K[i]
 dd := (u) -> add(coeff(u,e_K,i) * de_K[i],i=1..2^d);

 # Basis() calculates Groebner basis
 B0 := Basis([seq(ee_K[i],i=1..2^d)],T,char);
 K_basis := [seq(e_K^i,i=1..2^d)];
 K_rank  := 2^d;
 Z_basis := select(u -> degree(u,e_K) <= 2^d,B0);
 Z_rank  := nops(Z_basis);
 B_basis := select(u -> u <> 0,expand(map(quo,B0,e_K^(2^d+1),e_K)));
 B_rank  := nops(B_basis);

 B1 := Basis([seq(e_B^i + B_basis[i],i=1..B_rank),
              seq(e_Z^i + Z_basis[i],i=1..Z_rank)],T,char);

 BZ_basis := [seq(NormalForm(e_B^i,B1,T,char),i=1..B_rank)];
 BZ_rank  := nops(BZ_basis);

 K := table():
 K["num_rels"] := d;
 K["A"] := A;
 K["T"] := T;
 K["P"] := P;
 K["ix_P"] := eval(ix_P);
 K["de_K"] := eval(de_K);
 K["K_basis"] := K_basis;
 K["K_rank"]  := K_rank;
 K["Z_basis"] := Z_basis;
 K["Z_rank"]  := Z_rank;
 K["B_basis"] := B_basis;
 K["B_rank"]  := B_rank;
 K["BZ_basis"] := BZ_basis;
 K["BZ_rank"]  := BZ_rank;

 K["Z_basis_matrix"] :=
  Matrix([seq([seq(coeff(u,e_K,i),i=1..2^d)],u in Z_basis)]);
 K["B_basis_matrix"] :=
  Matrix([seq([seq(coeff(u,e_K,i),i=1..2^d)],u in B_basis)]);
 K["BZ_basis_matrix"] :=
  Matrix([seq([seq(coeff(u,e_Z,i),i=1..Z_rank)],u in BZ_basis)]);

 K["B0"] := B0;
 K["B1"] := B1;

 return eval(K);
end:


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 *)

$define lda_c_pw_params
$include "lda_c_pw.mpl"

$define lda_x_params
$include "lda_x.mpl"

cc06_cnst  := (3/(4*Pi))^(2/3):

cc06_alpha := -0.0007:
cc06_beta  :=  0.0080*cc06_cnst:
cc06_gamma :=  0.026 *cc06_cnst:

cc06_f := (rs, z, us0, us1) ->
  (f_lda_x(rs, z) + f_pw(rs, z))*(1 +
    (cc06_alpha + cc06_beta*u_total(z, us0, us1))/(1 + cc06_gamma*u_total(z, us0, us1))
  ):

f := (rs, z, xt, xs0, xs1, us0, us1, ts0, ts1) ->
  cc06_f(rs, z, us0, us1):


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 *)

$define gga_x_pbe_tca_params
$include "gga_x_pbe.mpl"

$define gga_x_pw91_params
$include "gga_x_pw91.mpl"

bpccac_malpha :=  1:
bpccac_mbeta  := 19:

bpccac_fab := x -> 1/(1 + exp(-bpccac_malpha*(x - bpccac_mbeta))):
bpccac_f   := x -> (1 - bpccac_fab(x))*pbe_f(x) + bpccac_fab(x)*pw91_f(x):

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


text:
`is_element/tree_Fbar_alt` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 local TT1,P,children,T,UU,UU1;
 global reason;

 TT1 := select(T -> nops(T) > 1,TT);

 if not(type(x,table)) then
  reason := [convert(procname,string),"x is not a table",x];
  return false;
 fi;

 P := sort(map(sort,map(op,{indices(x)})));

 if P <> TT1 then
  reason := [convert(procname,string),"x is not indexed by the big sets in TT",P,TT1];
  return false;
 fi;

 children := children_map(A)(TT);

 for T in TT1 do
  UU := children[T];
  UU1 := select(U -> (nops(U) > 1),UU);

  if not(`is_element/D/Fbar`(N)(UU1,UU)(x[T])) then
   reason := [convert(procname,string),"x[T] is not in D(C1,C)",x[T],N,T,reason];
   return false;
  fi;
 od;

 return true;
end;

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

`is_equal/tree_Fbar_alt` := (N::posint) -> (A::set) -> (TT) -> proc(x,y)
 local TT1,P,children,T,UU,UU1;

 TT1 := select(T -> nops(T) > 1,TT);

 children := children_map(A)(TT);

 for T in TT1 do
  UU := children[T];
  UU1 := select(U -> (nops(U) > 1),UU);

  if not(`is_equal/D/Fbar`(N)(UU1,UU)(x[T],y[T])) then
   return false;
  fi;
 od;

 return true;
end;

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

`random_element/tree_Fbar_alt` := (N::posint) -> (A::set) -> (TT) -> proc()
 local x,TT1,children,T,UU,UU1;

 x := table();

 TT1 := select(T -> nops(T) > 1,TT);
 children := children_map(A)(TT);

 for T in TT1 do
  UU := children[T];
  UU1 := select(U -> (nops(U) > 1),UU);
  x[T] := `random_element/D/Fbar`(N)(UU1,UU)();
 od;

 return eval(x);
end;

######################################################################
# From tree_Fbar_alt to tree_Fbar

`theta/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(ty)
 local TT1,children,z,t,y,C,n,pi,w,w0,a,m,T,U;

 TT1 := select(T -> nops(T) > 1,TT);
 TT1 := sort([op(TT1)],(U1,U2) -> (nops(U1) < nops(U2)));
 children := children_map(A)(TT);
 z := table();
 t := table();
 y := table();

 for T in TT1 do
  z[T] := table();
  t[T],y[T] := op(ty[T]);
  C := children[T];
  n := combine(simplify(sqrt(add(`norm_2/R`(N)(y[T][a])^2,a in C))));
  pi := table();
  for U in C do
   for a in U do
    pi[a] := U;
   od;
  od;
  w := table();
  w0 := [0$N];
  for a in T do
   w[a] := y[T][pi[a]];
   w0 := combine(simplify(w0 +~ w[a]));
  od;
  w0 := combine(simplify(w0 /~ nops(T)));
  for a in T do
   w[a] := combine(simplify(w[a] -~ w0));
  od;
  m := combine(simplify(sqrt(add(`norm_2/R`(N)(w[a])^2,a in T))));
  if is(m > 0) then
   for a in T do
    w[a] := combine(simplify((n/m) *~ w[a]));
   od;
  fi;
  for U in C do
   if nops(U) > 1 then
    for a in U do
     w[a] := combine(simplify(w[a] +~ t[T][U] *~ z[U][a]));
    od;
   fi;
  od;
  z[T] := eval(w);
 od:

 return eval(z);
end;

######################################################################
# From tree_Fbar to tree_Fbar_alt

`phi/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 local TT1,children,x1,t,y0,y,ty,C,C1,n,m,T,U,a;

 TT1 := select(T -> nops(T) > 1,TT);
 children := children_map(A)(TT);
 x1 := table();
 t  := table();
 y0 := table();
 y  := table();
 ty := table();

 for T in TT1 do
  x1[T] := `normalise_2/SW`(N)(T)(x[T]);
  t[T]  := table();
  y0[T] := table();
  y[T]  := table();
  C := children[T];
  C1 := select(U -> nops(U) > 1,C);
  n := 0;
  for U in C do
   y0[T][U] := simplify(expand(`sum/vector_functions`(N)(U)(x1[T]) /~ nops(U)));
   n := simplify(expand(n + `norm_2/R`(N)(y0[T][U])^2));
   if nops(U) > 1 then
    t[T][U]  := simplify(expand(sqrt(add(`norm_2/R`(N)(x1[T][a] -~ y0[T][U])^2,a in U))));
   fi;
  od;
  if is(n > 0) then
   y[T] := `normalise_2/SW`(N)(C)(y0[T]);
   m := combine(simplify(expand(rationalize(sqrt(1 - add(t[T][U]^2,U in C1))))));
   for U in C do
    y[T][U] := combine(simplify(expand(m *~ y[T][U])));
   od;
  else
   for U in C do
    y[T][U] := [0$N];
   od;
  fi;
  ty[T] := [eval(t[T]),eval(y[T])];
 od:

 return eval(ty);
end:


text:
# This file sets up the theory of the simplicial category, with simplicial
# maps encoded as lists of values.  This is an alternative to simplicial.mpl,
# where simplicial maps are encoded as tables.  This version is more efficient,
# but suffers from the fact that Maple uses 1-based indexing, so we need a
# lot of extra +1's.

`is_element/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> proc(f)
 evalb(type(f,list(nonnegint)) and nops(f) = n+1 and max(op(f)) <= m)
end:

`is_leq/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> proc(f,g)
 evalb( min(op(g -~ f)) >= 0 );
end:

# This assumes that
# *  p is a subset of {1,...,n} of size k
# *  q is a subset of {0,...,m} of size k+1
# It returns the nondecreasing map f : {0,...,n} -> {0,...,m} such that
# *  The image of f is q
# * { i > 0 : f(i) > f(i-1) } = p.

`build/simplicial_maps_alt` :=  (n::nonnegint,m::nonnegint) -> proc(k,p,q)
 local pp,i,j,f;
 pp := [0,op(p),n+1];
 f := NULL;
 i := 0;
 for i from 0 to k do
  for j from pp[i+1] to pp[i+2]-1 do
   f := f,q[i+1];
  od;
 od;

 return [f];
end:

`random_element/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> proc()
 local p,q,i,j,k,l,f;

 k := rand(0..min(n,m))();
 p := combinat[randcomb](n,k);
 q := combinat[randcomb](m+1,k+1) -~ 1;
 return `build/simplicial_maps_alt`(n,m)(k,p,q);
end:

`list_elements/simplicial_maps_alt` := proc(n::nonnegint,m::nonnegint)
 option remember;
 local i;

 if n = 0 then
  return [seq(`C/simplicial_maps_alt`(0,i),i=0..m)];
 else
  if m = 0 then
   return [`C/simplicial_maps_alt`(n,0)];
  else
   return [op(map(`I/simplicial_maps_alt`(n,m-1),
                   `list_elements/simplicial_maps_alt`(n,m-1))),
           op(map(`T/simplicial_maps_alt`(n-1,m),
	           `list_elements/simplicial_maps_alt`(n-1,m)))];
  fi;
 fi;
end:

`count_elements/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> binomial(n+m+1,m);

# Constant maps
`C/simplicial_maps_alt` := proc(n::nonnegint,k::nonnegint)
 [k$(n+1)]
end:

# Inclusion Delta(n,m) -> Delta(n,m+1)
`I/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> proc(f)
 f;
end:

# Map Delta(n,m) -> Delta(n+1,m): extend by sending n+1 to m
`T/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> proc(f)
 [op(f),m];
end:

# `delta/simplicial_maps_alt`(n)(i) : [n] >-> [n+1]; image omits i
`delta/simplicial_maps_alt` := proc(n::nonnegint,i::nonnegint)
 local j;
 if i > n+1 then return FAIL; fi;
 return [seq(j,j=0..i-1),seq(j+1,j=i..n)];
end:

# `sigma/simplicial_maps_alt`(n)(i) : [n] ->> [n-1]; takes the value i twice
`sigma/simplicial_maps_alt` := proc(n::nonnegint,i::nonnegint)
 local j;
 if i > n-1 then return FAIL; fi;
 return [seq(j,j=0..i),seq(j-1,j=i+1..n)];
end:

`id/simplicial_maps_alt` := proc(n::nonnegint)
 local i;
 return [seq(i=i,i=0..n)];
end:

`eval/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> (f) -> (i) -> f[i+1];

# `compose/simplicial_maps`(n,m,p)(f,g) assumes that f:[n] -> [m] and g:[m] -> [p],
# and it returns the composite g o f : [n] -> [p].

`compose/simplicial_maps_alt` := (n::nonnegint,m::nonnegint,p::nonnegint) -> proc(f,g)
 local i;
 [seq(g[f[i+1]+1],i=0..n)];
end:

`is_wide/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> proc(f)
 return evalb(f[1] = 0 and f[n+1] = m);
end:

`is_mono/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> proc(f)
 evalb(nops({op(f)}) = n+1);
end:

`is_epi/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> proc(f)
 local j;
 evalb({op(f)} = {seq(j,j=0..m)});
end:

# This assumes that f : [n] -> [m], and it returns [k,g,h] such that
# g : [n] ->> [k] and h : [k] >-> [m] and f = h o g

`factor/simplicial_maps_alt` := (n::nonnegint,m::nonnegint) -> proc(f)
 local K,k,g,h,hi,i;

 h := sort([op({op(f)})]);
 k := nops(h) - 1;
 hi := table();
 for i from 0 to k do hi[h[i+1]] := i; od;
 g := [seq(hi[f[i+1]],i=0..n)];
 return [k,g,h];
end:

`describe/simplicial_maps_alt` := (n::nonnegint) -> proc(f)
 local i;
 cat(seq(nat_code[f[i]],i=1..n+1));
end:

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

`is_element/simplicial_epi_alt` := (n::nonnegint,m::nonnegint) -> proc(f)
 global reason;

 if not(`is_element/simplicial_maps_alt`(n,m)(f)) then
  reason := [convert(procname,string),"not a simplicial map from [n] to [m]",reason];
  return false;
 fi;

 return `is_epi/simplicial_maps_alt`(n,m)(f);
end:

`is_leq/simplicial_epi_alt` :=  (n::nonnegint,m::nonnegint) -> proc(f,g)
 `is_leq/simplicial_maps_alt`(n,m)(f,g);
end:

`random_element/simplicial_epi_alt` := (n::nonnegint,m::nonnegint) -> proc()
 local p,q,k,f,i;

 if n < m then return FAIL; fi;

 k := m;
 p := combinat[randcomb](n,k);
 q := [seq(i,i=0..m)];
 return `build/simplicial_maps_alt`(n,m)(k,p,q);
end:

`list_elements/simplicial_epi_alt` := proc(n::nonnegint,m::nonnegint)
 local P,L,S,i,j;

 P := combinat[choose]([seq(i,i=1..n)],m);
 L := map(S -> [0$S[1],seq(j$(S[j+1]-S[j]),j=1..m-1),m$(n+1-S[m])],P);

 return L;
end:

`count_elements/simplicial_epi_alt` := (n::nonnegint,m::nonnegint) -> binomial(n,m);

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

`is_element/simplicial_mono_alt` := (n::nonnegint,m::nonnegint) -> proc(f)
 global reason;

 if not(`is_element/simplicial_maps_alt`(n,m)(f)) then
  reason := [convert(procname,string),"not a simplicial map from [n] to [m]",reason];
  return false;
 fi;

 return `is_mono/simplicial_maps_alt`(n,m)(f);
end:

`is_leq/simplicial_mono_alt` :=  (n::nonnegint,m::nonnegint) -> proc(f,g)
 `is_leq/simplicial_maps_alt`(n,m)(f,g);
end:

`random_element/simplicial_mono_alt` := (n::nonnegint,m::nonnegint) -> proc()
 local p,q,i,j,k,l,f;

 if n > m then return FAIL; fi;

 q := combinat[randcomb](m+1,n+1) -~ 1;
 q := sort([op(q)]);
 return q;
end:

`list_elements/simplicial_mono_alt` := proc(n::nonnegint,m::nonnegint)
 map(q -> q -~ 1,combinat[choose](m+1,n+1));
end:

`count_elements/simplicial_mono_alt` := (n::nonnegint,m::nonnegint) -> binomial(m+1,n+1);


text:
UpdateArchive := proc(lib_::string:="ppaml.mla")
  local libdir, drs, olddir, do_read, do_save, prev, fn, lib := lib_;
  drs := kernelopts(dirsep);

  libdir := eval(`if`(nargs=2,'args[2]','currentdir()'));
  if nargs=2 then
    libdir := args[2];
  else
    libdir := LibraryTools:-FindLibrary(Hakaru);
    if libdir <> NULL then
      libdir := FileTools:-ParentDirectory(libdir);
    else
      libdir := currentdir();
    end if;
  end if;
  olddir := currentdir(libdir):

  if not FileTools:-Exists(lib) then
    LibraryTools:-Create(lib):
  else
    try
      FileTools:-Remove(lib)
    catch:
      map(curry(march,'delete',lib)@curry(op,1),march('list', lib));
      if nops( march('list', lib) ) <> 0 then
        WARNING("Failed to delete all contents of %1",lib);
      end if;
    end try;
  end if:

  do_read := proc(x,$) read(cat(".",drs,x,".mpl")); NULL; end proc;
  do_save := proc(x,$) LibraryTools:-Save(x,lib); NULL; end proc;

  map(do_read,
   ["BindingTools", "Hakaru", "KB", "Partition", "Loop"
   ,"Domain", "NewSLO", "Summary"]);

  prev := kernelopts(opaquemodules=false):
  unprotect(Hakaru:-UpdateArchive);
  Hakaru:-UpdateArchive := copy(procname);
  protect(Hakaru:-UpdateArchive);

  map(do_save,
   [`gensym`, 'BindingTools', `depends/lam`, `depends/Branch`, `depends/Bind`, `depends/ary`, `depends/Plate`, `eval/lam`, `eval/Branch`, `eval/Bind`, `eval/ary`, `eval/Plate`, 'Hakaru', 'KB', `depends/forall`, `depends/Ints`, `depends/Sums`, `depends/ints`, `depends/sums`, `eval/forall`, `eval/Ints`, `eval/Sums`, `eval/ints`, `eval/sums`, `eval/Int`, `eval/Sum`, `eval/Product`, `eval/int`, `eval/sum`, `eval/product`, 'Loop', `depends/Integrand`, `depends/LO`, `eval/Integrand`, `eval/LO`, 'Domain', 'NewSLO', 'Partition', `depends/Bucket`, `depends/Index`, `eval/Bucket`, `eval/Index`, 'Summary', 'Bind', seq(parse(cat("`print/",q,"`")),q=[ints,sums,Ints,Sums]) ]);

  Domain:-Improve:-ModuleLoad():
  kernelopts(opaquemodules=prev):
  do_save('Domain:-Improve'):
  currentdir(olddir);
  NULL;
end proc:

UpdateArchive():

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
TITLE

    DakotaModel_deterministic;
    ! Objective value 4165.0
    ! line 23 has a syntax error

INDEX

    PRODUCTS := (Desk, Table, Chair);
    RESOURCES := (Lumber, Finishing, Carpentry);

DATA

    Price[PRODUCTS] := (60, 40, 10);
    Cost[RESOURCES] := (2, 4, 5.2);
    Consumption[RESOURCES, PRODUCTS] := (8, 6, 1
					4, 2, 1.5
					2, 1.5, .5);
    Demand[PRODUCTS] := (150, 125, 300);

INTEGER VARIABLES

   ProductionPRODUCTS];

VARIABLES

   Purchase[RESOURCES];

MAX

    Profit = sum(PRODUCTS: Price * Production) - sum(RESOURCES: Cost * Purchase);

SUBJECT TO

    ResourceConstraint[RESOURCES]: sum(PRODUCTS: Consumption * Production) <= Purchase;
    DemandConstraint[PRODUCTS]: Production <= Demand;

END

max_stars_repo_name: pashtetGP/optconvert

max_issues_repo_name: pashtetGP/optconvert

max_forks_repo_name: pashtetGP/optconvert


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 *)

kappa := 0.8040:
mu    := 0.2195149727645171:
m     := 100:

gamm  := m*mu/kappa:
Cx    := kappa/m:

f0 := s -> 1 + add(Cx * (gamm*s^2/(1 + gamm*s^2))^i, i=1..m):

f  := x -> f0(X2S*x):


text:
cdir := "/home/scott/Documents/research/software/hqca/hqca/maple/test_cases/":
with(LinearAlgebra): with(ArrayTools): with(QuantumChemistry):
loaddata := readdata(cat(cdir,"tmp7bh4a4ao.rdm"), 8):
Flatten := proc(x) local n, a, i, j, k, l; `local`(a, n, i, j, k, l); `description`("convert chemists to numpy and flatten an array to form a matrix"); n := Size(x); a := Matrix(1 .. n[1]*n[2], 1 .. n[3]*n[4], datatype = float[8]); for i to round(n[1]) do for j to round(n[2]) do for k to round(n[3]) do for l to round(n[4]) do a[(i - 1)*n[1] + j, (k - 1)*n[3] + l] := x[i, j, k, l]; end do; end do; end do; end do; return a; end proc:
New := Array(1 .. 4, 1 .. 4, 1 .. 4, 1 .. 4, datatype = float[8]):
for i in loaddata[3 .. ()] do
    New[round(i[1]), round(i[2]), round(i[3]), round(i[4])] := i[5]:
end do:
Transpose(Re(Eigenvalues(Flatten(New))));
pure := Purify2RDM(New, spin_free = false, electron_number=4, conv_tol = 0.00000000001):
ExportMatrix(cat(cdir, "_temp_purified.csv"), Flatten(pure[rdm2])):

max_stars_repo_name: damazz/HQCA

max_issues_repo_name: damazz/HQCA

max_forks_repo_name: damazz/HQCA


text:
# This is a row reduction algorithm for matrices over the p-local integers.
# It returns a triple [L,P,x].  Here P is an invertible matrix with P.M = L,
# so the row span of L is the same as that of M.  All rows of zeros in L
# appear at the bottom.  Every nonzero row contains an entry called a pivot.
# All pivot entries are powers of p.  If a pivot entry is p^v, then
# + all entries to the left are divisible by p^(v+1)
# + all entries to the right are divisible by p^v
# + all entries above or below are in the range {0,...,p^v - 1}
# These conditions force many entries to be zero.  I think that
# they characterise the result up to reordering the rows.  We use the
# following ordering:
# + If two rows have different valuations, then the one with lower
#    valuation comes first
# + If two rows have the same valuation, then the one with the pivot
#    further to the left comes first.
# The component x in the return value gives information abot the pivots.
# Explicitly, if the i'th entry in x is [j,v] then there is a pivot entry
# of p^v in position (i,j).

Zpl_reduce := proc(M,p)
 local L,L1,P,P1,r,c,i,j,k,h,m,J,K,v,v_min,v_min_pos,pivots,cmp,x;
 L := Copy(M);
 r,c := Dimension(M);
 P := Matrix(r,r);
 for i from 1 to r do P[i,i] := 1; od;
 v := Matrix(r,c);
 i := 0;
 J := {seq(j,j=1..r)};
 K := {seq(k,k=1..c)};
 pivots := NULL;
 v_min_pos := [0,0];
 while v_min_pos <> NULL do
  v_min := infinity;
  v_min_pos := NULL;
  for k in K do
   for j in J do
    v[j,k] := padic_val(L[j,k],p);
    if v[j,k] < v_min then
     v_min := v[j,k];
     v_min_pos := [j,k];
    fi;
   od;
  od;
  if v_min_pos = NULL then
   break;
  fi;
  j,k := op(v_min_pos);
  pivots := pivots,[j,k,v_min];
  m := p^v_min/L[j,k];
  RowOperation(L,j,m,inplace=true);
  RowOperation(P,j,m,inplace=true);
  for h from 1 to r do
   if h <> j then
    m := (L[h,k] - modp(L[h,k],p^v_min))/p^v_min;
    RowOperation(L,[h,j],-m,inplace=true);
    RowOperation(P,[h,j],-m,inplace=true);
   fi;
  od;
  J := sort(J minus {j});
  K := sort(K minus {k});
 od;
 cmp := proc(x,y)
  if x[3] < y[3] then return true; fi;
  if x[3] > y[3] then return false; fi;
  return evalb(x[2] < y[2]);
 end:
 pivots := sort([pivots],cmp);
 L1 := convert(L,listlist);
 P1 := convert(P,listlist);
 L1 := Matrix([seq(L1[x[1]],x in pivots),seq(L1[j],j in J)]);
 P1 := Matrix([seq(P1[x[1]],x in pivots),seq(P1[j],j in J)]);
 return [L1,P1,[seq([x[2],x[3]],x in pivots)]];
end:

######################################################################
# This is an implementation of the p-local integers that can be
# used with LinearAlgebra[Generic][SmithForm].  However, Smith
# forms are not unique, and it does not seem easy to control which
# version we get if we use LinearAlgebra[Generic][SmithForm].
# The above function Zpl_reduce does essentially the same job
# but with a more precisely specified result.

Zpl[`0`] := 0;
Zpl[`1`] := 1;
Zpl[`+`] := `+`;
Zpl[`-`] := `-`;
Zpl[`*`] := `*`;
Zpl[`=`] := `=`;

Zpl[Quo] := proc(a,b,r)
 local b0,q0,r0;
 b0 := p^padic_val(b,p);
 r0 := modp(a,b0);
 q0 := (a - r0)/b;
 if nargs >= 3 then r := r0; fi;
 return q0;
end:

Zpl[Rem] := proc(a,b,q)
 local b0,q0,r0;
 b0 := p^padic_val(b,p);
 r0 := modp(a,b0);
 q0 := (a - r0)/b;
 if nargs >= 3 then q := q0; fi;
 return r0;
end:

Zpl[EuclideanNorm] := ((m) -> padic_val(m,p));
Zpl[UnitPart] := ((m) -> m/p^padic_val(m,p));
Zpl[Gcdex] := proc(a,b,s,t)
 local u,v;
 u := padic_val(a,p);
 v := padic_val(b,p);
 if u <= v then
  if nargs >= 4 then
   s := p^u/a; t := 0;
  fi;
  return p^u;
 else
  if nargs >= 4 then
   s := 0; t := p^v/b;
  fi;
  return p^v;
 fi;
end:


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

# Real numbers mod integers

`is_element/RZ` := (x) -> type(x,realcons):

`is_equal/RZ` := (x,y) -> type(simplify(y-x),integer):

`is_leq/RZ` := NULL:

`normalize/RZ` := (x) -> simplify(x - floor(x)):

`random_element/RZ` := proc() rand(0..719)()/(720); end:

`list_elements/RZ` := NULL;
`count_elements/RZ` := NULL;

`d/RZ` := (x,y) -> abs(x - y - round(x - y));

`eta/RZ` := (x) -> exp(2*Pi*I*x);

`is_cyclic/RZ` := proc(x::list)
 local n,y,z,i;
 n := nops(x);
 if n <= 2 then return true; fi;
 y := simplify(x -~ x[1]);
 y := simplify(y -~ map(floor,y));
 z := [seq(y[i+1] - y[i],i=1..n-1)];
 return evalb(simplify(min(op(z))) >= 0);
end:

`is_strictly_cyclic/RZ` := proc(x::list)
 local n,y,z,i;
 n := nops(x);
 if n <= 2 then return true; fi;
 y := simplify(x -~ x[1]);
 y := simplify(y -~ map(floor,y));
 z := [seq(y[i+1] - y[i],i=1..n-1)];
 return evalb(simplify(min(op(z))) >= 0);
end:


text:
program gcd3;
procedure sort_pair_descending(pair: array[2] of integer);
begin
    if pair[0] < pair[1] then
    begin
        var tmp: integer;
        tmp := pair[0];
        pair[0] := pair[1];
        pair[1] := pair[0];
    end;
end;

procedure gcd(pair: array[2] of integer, var result: integer);
begin
    sort_pair_descending(pair);

    var new: integer;
    new := pair[0] % pair[1];
    pair[0] := pair[1];
    pair[1] := new;

    if pair[1] > 0 then
        return gcd(pair, result) {* semicolons are only used between statements in a block *}
    else
        result := pair[0]
    ;
end;

begin
    var pair: array[2] of integer;
    read(pair[0]);
    read(pair[1]);

    var result: integer;
    gcd(pair, result);
    writeln(result);
end.

max_stars_repo_name: sebschmi/codegen-course

max_issues_repo_name: sebschmi/codegen-course

max_forks_repo_name: sebschmi/codegen-course


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 *)
(* prefix:
  gga_k_mpbe_params *params;

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

$include "gga_x_mpbe.mpl"

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


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

lines := StringTools:-Split(input):
matchlookup := table( [ "]"="[", "}"="{", ">"="<", ")"="(",
    "["="]", "{"="}", "<"=">", "("=")" ] );
scores :=  table( [ "]"=57, "}"=1197, ">"=25137, ")"=3 ] ):
scores2 := table( ["("=1, "["=2, "{"=3, "<"=4] ):

points := 0:
badlines := NULL:
for j to nops(lines) do

    bracestack := DEQueue();
    for i to length(lines[j]) do
        if lines[j][i] in { "[", "{", "<", "(" } then
            push_back(bracestack, lines[j][i]);
        elif lines[j][i] in { "]", "}", ">", ")" } then
            c := pop_back(bracestack);
            if not c = matchlookup[lines[j][i]] then
                printf("Expected %a, but found %a instead.\n",
                    matchlookup[c], lines[j][i] );
                points += scores[lines[j][i]];
                badlines := badlines, j;
                break;
            end if;
        end if;
    end do;

    acscore[j] := 0;
    while not empty(bracestack) do
        c := pop_back(bracestack);
        acscore[j] := 5*acscore[j] + scores2[c];
    end do;
end do:

answer1 := points;
answer2 := round(Statistics:-Median([seq( acscore[i],
    i in ({ seq(1..nops(lines)) } minus {badlines}))] ));


text:
var $extern_var extern i32 public

func $test_extern_var ( var %i i32 ) i32 {
  return (dread i32 $extern_var) }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
######################################################################
# SW(N)(A) is the sphere associated to the vector space W(N)(A).
# In the main LaTeX document this is called S(NWA).
#
# Elements are represented as nonzero elements of the space W(N)(A),
# with an appropriate equivalence relation built into the definition
# the function `is_equal/SW`.

`is_element/SW` := (N::posint) -> (A::set) -> proc(x)
 local x0,a;
 global reason;

 if not(`is_element/W`(N)(A)(x)) then
  reason := [convert(procname,string),"x is not in W(N)(A)",x,N,A,reason];
  return false;
 fi;

 x0 := `normalise/W`(N)(A)(x);

 for a in A do
  if x0[a] <> [0$N] then return true; fi;
 od;

 reason := [convert(procname,string),"x is zero",x];
 return false;
end;

`normalise/SW` := (N::posint) -> (A::set) -> proc(x)
 local x0,u,r,a;

 x0 := `normalise/W`(N)(A)(x);
 r := max(0,seq(op(map(abs,x0[a])),a in A));
 for a in A do
  x0[a] := x0[a] /~ r;
 od;

 return(eval(x0));
end;

`bottom_normalise/SW` := (N::posint) -> (A::set) -> proc(x)
 local x0,u,r,a;

 x0 := `bottom_normalise/W`(N)(A)(x);
 r := max(0,seq(op(map(abs,x0[a])),a in A));
 for a in A do
  x0[a] := x0[a] /~ r;
 od;

 return(eval(x0));
end;

`is_equal/SW` := (N::posint) -> (A::set) -> proc(x,y)
 local x0,y0,a;

 x0 := `normalise/SW`(N)(A)(x);
 y0 := `normalise/SW`(N)(A)(y);

 for a in A do
  if not(`is_equal/R`(N)(x0[a],y0[a])) then
   return false;
  fi;
 od;

 return true;
end;

`dist/SW` := (N::posint) -> (A::set) -> proc(x,y)
 local x0,y0,a,d;

 x0 := `normalise/SW`(N)(A)(x);
 y0 := `normalise/SW`(N)(A)(y);

 d := 0;
 for a in A do
  d := d + `d_infinity/R`(N)(x0[a],y0[a]);
 od;

 return d;
end;

`normalise_2/SW` := (N::posint) -> (A::set) -> proc(x)
 local x0,u,r2,r,a;

 x0 := table();
 u := [0$N];

 for a in A do
  x0[a] := x[a];
  u := u +~ x[a];
 od;
 u := u /~ nops(A);
 r2 := 0;
 for a in A do
  x0[a] := x0[a] -~ u;
  r2 := r2 + `norm_2/R`(N)(x0[a])^2;
 od;
 r := sqrt(r2);

 for a in A do
  x0[a] := x0[a] /~ r;
 od;

 return(eval(x0));
end;

`dist_2/SW` := (N::posint) -> (A::set) -> proc(x,y)
 local x0,y0,a,d;

 x0 := `normalise_2/SW`(N)(A)(x);
 y0 := `normalise_2/SW`(N)(A)(y);

 d := 0;
 for a in A do
  d := d + `d_2/R`(N)(x0[a],y0[a])^2;
 od;

 d := combine(simplify(expand(d)));
 d := combine(simplify(expand(sqrt(d))));
 return d;
end;

`is_leq/SW` := NULL;

`random_element/SW` := (N::posint) -> (A::set) -> proc()
 local x,ok;

 ok := false;

 while not(ok) do
  x := `random_element/W`(N)(A)();
  ok := `is_element/SW`(N)(A)(x);
 od:

 return eval(x);
end;

`list_elements/SW` := NULL;
`count_elements/SW` := NULL;


text:
make_example_surfaces := proc()
 global example_surface;
 local P,T,R,dp,a,u,i,j,k;

 dp := (u,v) -> add(u[i] * v[i],i=1..3):
 R := (x) -> [-x[1]/2+sqrt(3)*x[2]/2,-sqrt(3)*x[1]/2-x[2]/2,x[3]]:

 P := table():
 P[0] := `map/pants`(x -> x,`base/pants`):
 for i from 1 to 5 do P[i] := `sprout/pants`(P[i-1],modp(i-1,3)+1); od:
 T := table():
 T[0] := `pants_joiner/pants_tube`(P[0],2,P[1],3):
 T[1] := `pants_joiner/pants_tube`(P[2],1,P[3],2):
 T[2] := `pants_joiner/pants_tube`(P[4],3,P[5],1):
 example_surface[0] := display(
  seq(`rough_plot/pants`(op(X)),X in entries(P)),
  seq(`plot/pants_tube`(op(X)),X in entries(T))
 ):

 P := table():
 T := table():
 P[0] := eval(`base/pants`):
 P[1] := `sprout/pants`(P[0],1):
 T[0] := `pants_joiner/pants_tube`(P[0],2,P[1],3):
 T[1] := `pants_joiner/pants_tube`(P[0],3,P[1],2):
 example_surface[1] := display(
  seq(`rough_plot/pants`(op(X)),X in entries(P)),
  seq(`plot/pants_tube`(op(X)),X in entries(T))
 ):

 P := table(): T := table():
 P[0] := eval(`base/pants`):
 P[1] := `sprout/pants`(P[0],1):
 P[2] := `map/pants`((x) -> x +~ [0,3,0],P[0]):
 P[3] := `map/pants`((x) -> x +~ [0,3,0],P[1]):
 T[0] := `pants_joiner/pants_tube`(P[0],2,P[2],3):
 T[1] := `pants_joiner/pants_tube`(P[1],3,P[3],2):
 T[3] := `pants_joiner/pants_tube`(P[0],3,P[1],2):
 T[4] := `pants_joiner/pants_tube`(P[2],2,P[3],3):
 example_surface[2] := display(
  seq(`rough_plot/pants`(op(X)),X in entries(P)),
  seq(`plot/pants_tube`(op(X)),X in entries(T))
 ):

 P := table(): T := table():
 P[0] := eval(`base/pants`):
 P[1] := `sprout/pants`(P[0],1):
 T[0] := `pants_joiner/pants_tube`(P[0],2,P[0],3):
 T[1] := `pants_joiner/pants_tube`(P[1],2,P[1],3):
 example_surface[3] := display(
  seq(`rough_plot/pants`(op(X)),X in entries(P)),
  seq(`plot/pants_tube`(op(X)),X in entries(T))
 ):

 a := 2:
 u := table():
 P := table():
 T := table():
 u[0] := simplex_embedding([1,0,0,0]):
 u[1] := simplex_embedding([0,1,0,0]):
 u[2] := simplex_embedding([0,0,1,0]):
 u[3] := simplex_embedding([0,0,0,1]):
 for i from 0 to 3 do
  P[i] := table([
   "centre" = a *~ u[i],
   "circle_centres" = [seq(simplify(a *~ u[i] +~ (u[j] -~ dp(u[j],u[i]) *~ u[i]) *~ (sqrt(27/32))),j in {0,1,2,3} minus {i})],
   "top" = (a + 1/2) *~ u[i]
  ]);
 od:
 T[0,1] := `pants_joiner/pants_tube`(P[0],1,P[1],1):
 T[0,2] := `pants_joiner/pants_tube`(P[0],2,P[2],1):
 T[0,3] := `pants_joiner/pants_tube`(P[0],3,P[3],1):
 T[1,2] := `pants_joiner/pants_tube`(P[1],2,P[2],2):
 T[1,3] := `pants_joiner/pants_tube`(P[1],3,P[3],2):
 T[2,3] := `pants_joiner/pants_tube`(P[2],3,P[3],3):
 example_surface[4] := display(
  seq(`rough_plot/pants`(op(X)),X in entries(P)),
  seq(`plot/pants_tube`(op(X)),X in entries(T)),
  axes=none,scaling=constrained
 ):

 P := table():
 T := table():
 a := 2:
 for i in [-1,1] do
  for j in [-1,1] do
   for k in [-1,1] do
    P[i,j,k] := table([
     "centre" = a *~ [i,j,k],
     "circle_centres" = evalf(map(x -> (a *~ [i,j,k] +~ x /~ sqrt(8)),
			    [[-2*i,j,k],[i,-2*j,k],[i,j,-2*k]])),
     "top" = evalf((a + 1/(2*sqrt(3))) *~ [i,j,k])
    ]);
   od:
  od:
 od:
 T[ 0, 1, 1] := `pants_joiner/pants_tube`(P[ 1, 1, 1],1,P[-1, 1, 1],1):
 T[ 0, 1,-1] := `pants_joiner/pants_tube`(P[ 1, 1,-1],1,P[-1, 1,-1],1):
 T[ 0,-1, 1] := `pants_joiner/pants_tube`(P[ 1,-1, 1],1,P[-1,-1, 1],1):
 T[ 0,-1,-1] := `pants_joiner/pants_tube`(P[ 1,-1,-1],1,P[-1,-1,-1],1):
 T[ 1, 0, 1] := `pants_joiner/pants_tube`(P[ 1, 1, 1],2,P[ 1,-1, 1],2):
 T[ 1, 0,-1] := `pants_joiner/pants_tube`(P[ 1, 1,-1],2,P[ 1,-1,-1],2):
 T[-1, 0, 1] := `pants_joiner/pants_tube`(P[-1, 1, 1],2,P[-1,-1, 1],2):
 T[-1, 0,-1] := `pants_joiner/pants_tube`(P[-1, 1,-1],2,P[-1,-1,-1],2):
 T[ 1, 1, 0] := `pants_joiner/pants_tube`(P[ 1, 1, 1],3,P[ 1, 1,-1],3):
 T[ 1,-1, 0] := `pants_joiner/pants_tube`(P[ 1,-1, 1],3,P[ 1,-1,-1],3):
 T[-1, 1, 0] := `pants_joiner/pants_tube`(P[-1, 1, 1],3,P[-1, 1,-1],3):
 T[-1,-1, 0] := `pants_joiner/pants_tube`(P[-1,-1, 1],3,P[-1,-1,-1],3):
 example_surface[5] := display(
  seq(`rough_plot/pants`(op(X)),X in entries(P)),
  seq(`plot/pants_tube`(op(X)),X in entries(T)),
  axes=none,scaling=constrained
 ):

 return eval(example_surface);
end:


text:
# Teach Maple (through depends and eval) about our new binding forms.
# forall bind from 1st arg to 2nd arg.
# Ints,Sums,ints,sums bind from 2nd arg to 1st arg, and also from each element
#   of the 4th arg to the other elements on the left and to the 3rd arg.

`depends/forall` := proc(bvar, pred, x, $)
  depends(pred, x minus convert(bvar, 'set'))
end proc:

`depends/Ints` := proc(body, bvar, rng, loops, x, y, $)
  local xx, i;
  if nargs = 5 then
    xx := x; # y is ununsed and x lists the variables to look for
  else
    xx := y; # x is some KB and y lists the variables to look for
    if depends(x, y) then return true end if;
  end if;
  # don't remove bvar from xx!
  if depends(body, xx minus {bvar}) then return true end if;
  for i from nops(loops) to 1 by -1 do
    if depends(op([i,2],loops), xx) then return true end if;
    xx := xx minus {op([i,1],loops)};
  end do;
  depends(rng, xx)
end proc:
`depends/Sums` := eval(`depends/Ints`):
`depends/ints` := eval(`depends/Ints`):
`depends/sums` := eval(`depends/Ints`):

`eval/forall` := proc(e, eqs, $)
  local bvar, pred;
  bvar, pred := op(e);
  eval(op(0,e), eqs)(BindingTools:-generic_evalat(bvar, pred, eqs))
end proc:

`eval/Ints` := proc(e, eqs, $)
  local body, bvar, rng, loops, n, i;
  body, bvar, rng, loops := op(1..4, e);
  bvar, body := BindingTools:-generic_evalat(bvar, body, eqs);
  eval(op(0,e), eqs)(body, bvar,
                     BindingTools:-generic_evalatstar(rng, loops, eqs),
                     op(5..-1,e))
end proc:
`eval/Sums` := eval(`eval/Ints`):
`eval/ints` := eval(`eval/Ints`):
`eval/sums` := eval(`eval/Ints`):

`eval/Int` := proc(e, eqs, $)
  local body, bound, bvar;
  body, bound := op(1..2, e);
  if bound :: name then
    bound, body := BindingTools:-generic_evalat(bound, body, eqs);
  elif bound :: `=` then
    bvar := lhs(bound);
    bvar, body := BindingTools:-generic_evalat(bvar, body, eqs);
    bound := bvar = eval(rhs(bound), eqs);
  else
    body, bound := BindingTools:-generic_evalatstar(body, bound, eqs);
  end if;
  eval(op(0,e), eqs)(body, bound, op(eval([op(3..-1,e)], eqs)))
end proc:
`eval/Sum`     := eval(`eval/Int`):
`eval/Product` := eval(`eval/Int`):
`eval/int`     := eval(`eval/Int`):
`eval/sum`     := eval(`eval/Int`):
`eval/product` := eval(`eval/Int`):

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

Loop := module ()
  option package;
  local intssums, enter_piecewise, wrap,
        t_peel, do_peel, do_split, do_graft, do_rebase_lower, do_rebase_upper,
        Binder, Stmt, t_binder, t_stmt, t_exp,
        ModuleLoad;
  export
     # These first few are smart constructors (for themselves):
         ints, sums,
     # while these are "proper functions"
         mk_HArray, genLoop, unproducts, unproduct,
         peel, split, graft, rebase_lower, rebase_upper, Print;
  # these names are not assigned (and should not be).  But they are
  # used as global names, so document that here.
  global Ints, Sums, csgn, sum;
  uses Hakaru, KB;

  t_binder := 'Binder(identical(product, Product, sum, Sum), t_kb)';
  t_stmt   := 'Stmt(anything, list, list)';
  t_exp    := '{Stmt(identical(exp), [], []),
                Stmt(identical(`^`), [anything], [])}';

  ints := proc() intssums('ints', 'int', _passed) end proc;
  sums := proc() intssums('sums', 'sum', _passed) end proc;

  Print := proc(kind::identical(Ints,Sums,ints,sums),e,v::name,r0::range,rs::list(name=range),$)
    local k;
    k := StringTools[Capitalize](StringTools[Chop](convert(kind,string)));
    k := parse(cat("`print/",eval(k),"`"));
    eval(k)(e,v[op(rs)]=r0);
  end proc;

  intssums := proc(makes::name, make::name,
                   e::anything, x::name, rr::range, ll::list(name=range),
                   kb::t_kb:=empty, $)
    local t, r, l, kb2, w0, pp, i;
    t := `if`(make=int, HReal(open_bounds(rr)), HInt(closed_bounds(rr)));
    r, l, kb2 := genLoop(rr, ll, kb, 'Integrand'(x,e));
    w0, pp := unproducts(e, x, l, kb2);
    if depends(w0, x) then
      'makes(e, x, rr, ll)'
    else
      pp := 'make'(pp, x=r);
      for i in l do pp := 'product'(pp, i) end do;
      NewSLO:-simplify_factor_assuming(w0 * pp, kb);
    end if
  end proc;

  mk_HArray := proc(t::t_type, loops::list(name=range), $)
    local res, i;
    res := t;
    for i in loops do res := HArray(res) end do;
    res
  end proc;

  genLoop := proc(e, loops::list(name=range), kb::t_kb)
    local kb1, rng, ind, do_subst, i;
    kb1 := kb;
    rng := table();
    ind := table();
    do_subst := e -> foldl(((e,eq) -> eval(e, op([lhs(eq),1],loops)=rhs(eq))),
                           e, entries(ind, 'pairs'));
    for i from nops(loops) to 1 by -1 do
      rng[i] := do_subst(op([i,2],loops));
      ind[i], kb1 := genType(op([i,1],loops), HInt(closed_bounds(rng[i])),
                               kb1, _rest);
    end do;
    do_subst(e), zip(`=`, [entries(ind, 'nolist')],
                          [entries(rng, 'nolist')]), kb1
  end proc;

  unproducts := proc(w, x::name, loops::list(name=range), kb::t_kb, $)
    local w0, pp, j, w1, w2, loop;
    w0 := 1;
    pp := w;
    for j from nops(loops) to 1 by -1 do
      w1, pp := op(unproduct(pp, x, op(j,loops), [], `*`, kb, kb));
      # separate out each of the factors in w1, as they might have different
      # variable dependencies, which can be exploited by other routines
      w2 := convert(w1, 'list', '`*`');
      for loop in [op(j+1..-1, loops)] do
        w2 := map((w -> product(eval(w, x=idx(x,lhs(loop))), loop)), w2);
      end do;
      w0 := w0 * `*`(op(w2));
      # w0 := w0 * foldl(product, w1, op(j+1..-1, loops));
    end do;
    # Rewrite ... * piecewise(i<=kk-1, xs^2, 1)
    #      to ... * xs ^ (2*piecewise(i<=kk-1, 1))
    # because the latter is easier to integrate and recognize with respect to xs
    pp := for_poly(pp, # TODO: move this rewrite into eval_factor (then remove for_poly?)?
      proc (p)
        local n, s, r;
        if p :: 'And(specfunc(piecewise), anyfunc(freeof(x), anything, 1))' then
          n := indets(op(1,p),name);
          if op(2,p) :: {'exp(anything)', '`^`'('freeof'(n),'anything')} then
            s, r := selectremove(depends, convert(op([2,-1],p),'list','`*`'), n);
            subsop(-1 = `*`(piecewise(op(1,p), `*`(op(s))), op(r)), op(2,p))
          else
            s, r := selectremove(depends, convert(op(2,p),'list','`*`'), n);
            `*`(op(r)) ^ piecewise(op(1,p), 1)
              * `if`(nops(s) > 0, piecewise(op(1,p), `*`(op(s)), 1), 1)
          end if
        else
          p
        end if
      end proc);
    w0, pp
  end proc;

  # Find [w1,pp] so that
  #   wrap(heap,w,mode,kb1,kb0)
  #   = w1*product(eval(pp,var=idx(var,lhs(loop))),loop)
  # making w1 depend on var as little as possible.
  # The flag "mode" should be `+` if "heap" contains an entry of the form
  # t_exp, or `*` otherwise.
  unproduct := proc(w, var::name, loop::(name=range),
                    heap::list, mode::identical(`*`,`+`),
                    kb1::t_kb, kb0::t_kb, $)
    local ind, res, dummy, kb, kbThen, i, w1, pp, s, r, x;
    if not depends(w, var) then
      return [wrap(heap, w, mode, kb1, kb0), 1]
    end if;
    # The `..' (indices) of all occurences of `idx(var,..)'
    ind := map2(op, 2, indets(w, idx(identical(var), anything)));
    if nops(ind) = 1 then
      ind := op(ind);
      # Make sure ind contains no bound variables before lifting it!
      # So, check that "extract using indets" and "rename using eval" commute.
      s := indets(ind, 'name');
      s := map(proc(x,$) local y; `if`(depends(ind,x), x=y, NULL) end proc, s);
      if indets(eval(w, s), idx(identical(var), anything))
                = {idx(var, eval(ind, s))} then
        # use kb as a local context, and 'solve' for the innermost bound var
        #   used in ind.
        kb  := assert(lhs(loop)=ind, kb1); # BUG! bijectivity assumed!
        # if the assertion makes the KB false, it isn't needed(?)
        ASSERT(type(kb, t_kb),
          sprintf("assert(%a,%a) produced a false KB!\n", lhs(loop)=ind,kb1)):
        res := subs(idx(var,ind) = idx(var,lhs(loop)), w);
        res := wrap(heap, res, mode, kb, kb0);
        res := subs(idx(var,lhs(loop))=dummy, res);
        if not depends(res, var) then
          return [1, subs(dummy=var, res)]
        end if
      end if
    end if;
    # distribute the unproduct over each part
    if w :: mode then
      res := map(unproduct, `if`(mode=`*`, list_of_mul(w), [op(w)]),
                 var, loop, heap, mode, kb1, kb0);
      return [`*`(op(map2(op,1,res))), `*`(op(map2(op,2,res)))]
    end if;
    # for piecewise, just map right in (using KB for context tracking)
    if w :: 'specfunc(piecewise)' then
      kb := kb1;
      for i from 1 to nops(w) do
        if i :: even then
          kbThen := assert(    op(i-1,w) , kb);
          # accumulate the rest, regardless
          kb     := assert(Not(op(i-1,w)), kb);

          # if kbThen is invalid, deal appropriately
          if not type(kbThen,t_kb) then next end if;

          w1[i], pp[i] := op(unproduct(op(i,w),var,loop,heap,mode,kbThen,kb0));
          # if kb is inconsistent, we're done
          if not type(kb,t_kb) then break end if;
        elif i = nops(w) then
          if not type(kb,t_kb) then next end if;
          w1[i], pp[i] := op(unproduct(op(i,w),var,loop,heap,mode,kb    ,kb0))
        end if
      end do;
      return [`*`(entries(w1,'nolist')), `*`(entries(pp,'nolist'))]
    end if;
    if mode = `*` then
      if w :: (anything^freeof(var)) then
        return unproduct(op(1,w), var, loop,
                         [op(heap), Stmt(`^`, [], [op(2,w)])], `*`, kb1, kb0)
      elif w :: exp(anything) then
        return unproduct(op(1,w), var, loop,
                         [op(heap), Stmt(exp, [], [])], `+`, kb1, kb0)
      elif w :: (freeof(var)^anything) then
        return unproduct(op(2,w), var, loop,
                         [op(heap), Stmt(`^`, [op(1,w)], [])], `+`, kb1, kb0)
      end if
    end if;
    if mode = `+` and w :: `*` then
      s, r := selectremove(depends, w, var);
      if s :: `*` then
        # Nonlinear %1 (time to reread kiselyov-lifted?)
      else
        return unproduct(s, var, loop,
                         [op(heap), Stmt(`*`, [], [r])], `+`, kb1, kb0)
      end if
    end if;
    if w :: And(specfunc(`if`(mode = `*`, {product, Product}, {sum, Sum})),
                anyfunc(anything, name=range(freeof(var)))) then
      x, kb := genType(op([2,1],w), HInt(closed_bounds(op([2,2],w))), kb1,
                       w, var, loop, heap);
      return unproduct(eval(op(1,w), op([2,1],w)=x), var, loop,
                       [op(heap), Binder(op(0,w), kb1)], mode, kb, kb0)
    end if;
    if mode = `*` and w :: '{`+`, specfunc({sum, Sum})}' then
      # Maybe this w is one of those big sums involving products that are
      # always equal to 1, left behind by the density of Categorical
      w1 := graft(split(peel(w)));
      w1 := combine(rebase_upper(w1));
      w1 := combine(rebase_lower(w1));
      return [wrap(heap, w1, mode, kb1, kb0), 1]
    end if;
    return [wrap(heap, w, mode, kb1, kb0), 1]
  end proc;

  enter_piecewise := proc(ee, kb0::t_kb, mode::identical(`*`,`+`), $)
    local e, kb, mo;
    e  := ee;
    kb := kb0;
    mo := mode();
    while e :: 'specfunc(piecewise)' and nops(e) = 3 do
      if op(3,e) = mo then
        kb := assert(op(1,e), kb);
        e := op(2,e);
      elif op(2,e) = mo then
        kb := assert(Not(op(1,e)), kb);
        e := op(3,e);
      else
        break;
      end if;
    end do;
    ASSERT(type(kb,t_kb), "enter_piecewise: KB of piecewise conditions is not valid, so the input piecewise must not be valid.");
    e, kb
  end proc;

  # heap is a list of t_binder and t_stmt
  # one invariant to maintain:
  #  wrap(heap, mode(a,b)) = wrap(heap,a) * wrap(heap,b)

  # Also, a t_binder is a Product/product/Sum/sum.
  # And, a t_stmt is a 1-hole context of a (multiplicative) AST,
  #   used for [a ^ hole, hole ^ b, exp(hole), c * hole].

  wrap := proc(heap::list, e1, mode1::identical(`*`,`+`),
               kb1::t_kb, kb0::t_kb, $)
    local e, kb, mode, i, entry, rest, var, new_rng, make, logmake,
       dom_spec, w, wPow, wExp, arrrgs;
    e    := e1;
    kb   := kb1;
    mode := mode1;
    for i from nops(heap) to 1 by -1 do
      entry := op(i,heap);
      if entry :: t_binder then
        if not (op(1,entry) in `if`(mode=`+`, {sum,Sum},
                                              {product,Product})) then
          print("Warning: heap mode inconsistency", heap, mode1)
        end if;
        e, kb := enter_piecewise(e, kb, mode);
        rest := kb_subtract(kb, op(2,entry));
        new_rng, rest := selectremove(type, rest,
          {[identical(genType), name, specfunc(HInt)],
           [identical(genLet), name, anything]});
        if not (new_rng :: [list]) then
          error "kb_subtract should return exactly one gen*"
        end if;
        new_rng := op(new_rng);
        var     := op(2,new_rng);
        if op(1,new_rng) = genType then
          make    := op(1,entry);
          new_rng := range_of_HInt(op(3,new_rng));
        else # op(1,new_rng) = genLet
          make    := eval;
          new_rng := op(3,new_rng);
        end if;
        dom_spec, rest := selectremove(depends,
          map(proc(a::[identical(assert),anything],$) op(2,a) end proc, rest),
          var);

        # Like e := make(piecewise(And(op(dom_spec)), e, mode()), var=new_rng);
        # but try to simplify by pushing the make and the piecewise into e
        (e, w) := selectremove(depends, convert(e, 'list', `*`), var);
        w := `*`(op(w));
        if mode = `*` then
          (wPow, e) := selectremove(type, e, '`^`'('freeof'(var), 'anything'));
          (wExp, e) := selectremove(type, e, 'exp(anything)');
          wExp := `+`(op(map2(op, 1, wExp)));
          logmake := t -> `if`(make=eval,eval,Sum)
                              (piecewise_And(dom_spec,t,0), var=new_rng);
          w := `*`( `if`(w = 1, 1, w ^ logmake(1))
                  , op(map((f -> op(1,f) ^ logmake(op(2,f))), wPow))
                  , `if`(wExp = 0, 1, exp(logmake(wExp))) );
        end if;
        e := w * make(piecewise_And(dom_spec, `*`(op(e)), mode()), var=new_rng);

        # `rest' originally comes from `entry', and this reinserts it after some
        # processing. The original KB was valid, so this one should also be valid.
        kb := foldr(assert, op(2,entry), op(rest));
        ASSERT(type(kb,t_kb), "Loop/wrap{t_binder}: rebuilding kb produced contradiction");

      elif entry :: t_stmt then
        # We evaluate arrrgs first, in case op(1,stmt) is an operation (such as
        # piecewise) that looks ahead at how many arguments it is passed before
        # they are evaluated.
        arrrgs := op(op(2,entry)), e, op(op(3,entry));
        e      := op(1,entry)(arrrgs);
        if entry :: t_exp then
          if mode <> `+` then
            print("Warning: heap mode inconsistency?", heap, mode1)
          end if;
          mode := `*`;
        end if
      else error "Malformed heap entry %1", entry end if
    end do;
    if mode <> `*` then
      print("Warning: heap mode inconsistency??", heap, mode1)
    end if;
    e, kb := enter_piecewise(e, kb, mode);
    rest := kb_subtract(kb, kb0);
    rest := map(proc(a::[identical(assert),anything],$) op(2,a) end proc, rest);
    piecewise_And(rest,e,mode())
  end proc;

  # Rewrite product(...piecewise(i+42=lo+42,th,el)...,i=lo..hi)
  # to eval(...th...,i=lo)*product(...el...,i=lo+1..hi)
  t_peel := 'And(specfunc({sum,Sum,product,Product}),
                 anyfunc(anything, name=range))';
  peel := proc(e, $)
    subsindets(e, t_peel, do_peel);
  end proc;
  do_peel := proc(e, $)
    local make, body, x, r, cond, line, test, x0, r0, cond0, here, rest;
    if not (e :: t_peel) then return e end if;
    if   op(0,e) in '{sum    ,Sum    }' then make := `+`;
    elif op(0,e) in '{product,Product}' then make := `*`; end if;
    body, r := op(e);
    x, r    := op(r);
    for cond in indets(body, '{`=`(algebraic),`<`,`<=`}') do
      line := `-`(op(cond));
      if 1 <> degree(line, x) then next end if;
      if cond :: `=` then
        cond0 := true;
        if Testzero(eval(line, x=lhs(r))) then
          x0 := lhs(r); r0 := applyop(`+`, 1, r, 1);
        elif Testzero(eval(line, x=rhs(r))) then
          x0 := rhs(r); r0 := applyop(`-`, 2, r, 1);
        else next end if;
      else
        test := proc(x0, x1, $)
          local res;
          res := evalb(not op(0,cond)(eval(line, x=x1), 0));
          `if`(res in '{true,false}'
                 and evalb(op(0,cond)(eval(line, x=x0), 0)) = res,
               res,
               FAIL)
        end proc;
        cond0 := test(lhs(r), lhs(r)+1);
        if cond0 <> FAIL then
          x0 := lhs(r); r0 := applyop(`+`, 1, r, 1);
        else
          cond0 := test(rhs(r), rhs(r)-1);
          if cond0 <> FAIL then
            x0 := rhs(r); r0 := applyop(`-`, 2, r, 1);
          else next end if;
        end if;
      end if;
      here := subs(cond =           cond0 , body);
      rest := subs(cond = evalb(not cond0), body);
      return make(eval(here, x=x0), do_peel(eval(subsop(1=rest, [2,2]=r0, e))));
    end do;
    return e;
  end proc;

  # Expand sum(a*(b-c),q) to sum(a*b,q)-sum(a*c,q)
  split := proc(e, $)
    subsindets(e, 'And(specfunc({sum,Sum}),
                       anyfunc(And(`*`,Not(`*`(Not(`+`)))),name=anything))',
                  do_split);
  end proc;
  do_split := proc(e, $)
    local terms, x;
    terms := convert(expand(op(1,e), op(indets(op(1,e), function))),
                     'list', `+`);
    x := op([2,1],e);
    `+`(op(map(proc(term, $)
      local s, r;
      s, r := selectremove(depends, convert(term, 'list', `*`), x);
      `*`(op(r), subsop(1=`*`(op(s)),e))
    end proc, terms)))
  end proc;

  # Simplify f(lo-1)*product(f(i),i=lo..hi) to product(f(i),i=lo-1..hi)
  graft := proc(e, $)
    subsindets(e, 'Or(And(`*`,Not(`*`(Not(specfunc({product,Product}))))),
                      And(`+`,Not(`+`(Not(specfunc({sum    ,Sum    }))))))',
                  do_graft);
  end proc;
  do_graft := proc(e, $)
    local produce, factors, i, j;
    produce := `if`(e::`*`, '{product,Product}',
                            '{sum    ,Sum    }');
    factors := sort(convert(e,'list'),
                    key = (factor -> -numboccur(factor,produce)));
    for i from nops(factors) to 2 by -1 do
      for j from i-1 to 1 by -1 do
        if op(j,factors) :: 'And'('specfunc'(produce),
                                  'anyfunc(anything,name=range)') then
          if Testzero(op(i,factors) - eval(op([j,1],factors),
               op([j,2,1],factors)=op([j,2,2,1],factors)-1)) then
            factors := subsop(i=NULL,applyop(`-`,[j,2,2,1],factors,1));
            break
          elif Testzero(op(i,factors) - eval(op([j,1],factors),
                 op([j,2,1],factors)=op([j,2,2,2],factors)+1)) then
            factors := subsop(i=NULL,applyop(`+`,[j,2,2,2],factors,1));
            break
          end if
        end if
      end do
    end do;
    op(0,e)(op(factors))
  end proc;

  # Normalize sum(f(i),i=2..hi) to sum(f(i+2),i=0..hi-2)
  rebase_lower := proc(e, $)
    subsindets(e, 'And(specfunc({sum,Sum,product,Product}),
                       anyfunc(anything,
                         name=Not({0,SymbolicInfinity,undefined})..anything))',
                  do_rebase_lower);
  end proc;
  do_rebase_lower := proc(e, $)
    subsop([2,2,1]=0,
           applyop(`-`,
                   [2,2,2],
                   applyop(eval, 1, e, op([2,1],e)=op([2,1],e)+op([2,2,1],e)),
                   op([2,2,1],e)))
  end proc;

  # Normalize sum(f(i),i=lo..2) to sum(f(i+2),i=lo-2..0)
  rebase_upper := proc(e, $)
    subsindets(e, 'And(specfunc({sum,Sum,product,Product}),
                       anyfunc(anything,
                         name=anything..Not({0,SymbolicInfinity,undefined})))',
                  do_rebase_upper);
  end proc;
  do_rebase_upper := proc(e, $)
    subsop([2,2,2]=0,
           applyop(`-`,
                   [2,2,1],
                   applyop(eval, 1, e, op([2,1],e)=op([2,1],e)+op([2,2,2],e)),
                   op([2,2,2],e)))
  end proc;

  ModuleLoad := proc($)
    # Override csgn to work a little bit harder on piecewise and sum
    # (to get rid of csgn(1/2+1/2*sum(piecewise(...,1,0),...))
    #  produced by int on a Gaussian mixture model)
    unprotect(csgn);
    csgn := overload([
      # Handle if the csgn of a piecewise doesn't depend on which branch
      proc(a :: specfunc(piecewise), $)
        option overload;
        local r, i;
        r := {seq(`if`(i::even or i=nops(a), csgn(op(i,a)), NULL),
                  i=1..nops(a))};
        if nops(r)=1 then return op(r) end if;
        if not assigned(_Envsignum0) then
          r := r minus {0};
          if nops(r)=1 then return op(r) end if;
        end if;
        error "invalid input: cannot csgn %1", a;
      end proc,
      # Handle if the csgn of a sum doesn't depend on the bound variable
      proc(a :: And(specfunc({sum, Sum}), anyfunc(anything, name=range)), $)
        option overload;
        local r;
        r := csgn(op(1,a));
        if not depends(r,op([2,1],a)) then
          return signum(op([2,2,2],a)+1-op([2,2,1],a)) * r
        end if;
        error "invalid input: cannot csgn %1", a;
      end proc,
      csgn]);
    protect(csgn);

    # Override sum to fail faster
    unprotect(sum);
    sum := overload([
      proc(f :: Not({`+`,`+`^integer}),
           k :: name=And(range,Not(range(rational))), $)
        option overload;
        if not (`-`(op(rhs(k))) :: 'rational')
           and not (subsindets(f, 'anything^integer', f->op(1,f))
                    :: '{`+`,And(`*`,Not(`*`(Not(`+`))))}')
           and depends(indets(f, '{specfunc(idx), specindex(idx)}'), op(1,k))
        then
          'procname(_passed)'
        else
          error "invalid input: cannot fast-fail sum(%1, %2)", _passed
        end if
      end proc,
      sum]);
    protect(sum);

    :-`print/Ints` := curry(thismodule:-Print,Ints);
    :-`print/Sums` := curry(thismodule:-Print,Sums);
    :-`print/ints` := curry(thismodule:-Print,ints);
    :-`print/sums` := curry(thismodule:-Print,sums);

  end proc;
  ModuleLoad():
end module; # Loop


text:
# validate new test case (Figure8_7.java)
with(simplex);

Constraints := [

# conservation of units at each node
e13+e14+e41+e15 = 400, # CHI
e21+e23+e24+e25 = 400, # DC

e13+e23     = 300, # HOU
e14+e41+e24 = 200, # ATL
e15+e25     = 300, # BOS

# maximum flow on individual edges
0 <= e13, e13 <= 200,
0 <= e14, e14 <= 350,
0 <= e41, e41 <= 200,
0 <= e21, e21 <= 200,
0 <= e23, e23 <= 280,
0 <= e15, e15 <= 200,
0 <= e24, e24 <= 200,
0 <= e25, e25 <= 350
];

Cost := 3*e41 + 7*e13 + 3*e14 + 4*e21 + 4*e23 + 7*e24 + 6*e15 + 6*e25;

# Invoke linear programming to solve problem
minimize (Cost, Constraints, NONNEGATIVE);


text:
program Write;	{sample24}
begin;
    writeln('It's OK?')
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing


text:
(*
 Copyright (C) 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: lda_exc *)
(* prefix:
  lda_c_chachiyo_mod_params *params;

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

(* Functional is based on Chachiyo correlation *)
$include "lda_c_chachiyo.mpl"
(* .. but with a different scaling function *)
g := z -> (opz_pow_n(z,2/3) + opz_pow_n(-z,2/3))/2:
g_zeta := zeta -> 2*(1 - g(zeta)^3):

f_chachiyo := (rs, zeta) -> e0(rs) + (e1(rs) - e0(rs))*g_zeta(zeta):
f := (rs, zeta) -> f_chachiyo(rs, zeta):