text:
# 交互式求解函数

$ifndef _INTERACTION_
$define _INTERACTION_

# 添加新的RepSol进行求解
resolveRep:=proc(rep,solInd,{newIsol:=[],newRep:=[],nocheck::boolean:=false})
    return resolveSol(rep:-osol[solInd],_options);
end proc:

# 添加新的InvSol进行求解
resolveSol:=proc(sol::InvSol,{newIsol:=[],newRep:=[],nocheck::boolean:=false})
    local s,nreps;
    if (newIsol=[]) and (newRep=[]) then
        error "至少设置newIsol和newRep之一的值";
    end if;
    s:=Object(sol);
    if (newIsol<>[]) then
        setIsol(s,newIsol,_options['nocheck']);
    end if;
    if (newRep<>[]) then
        setRep(s,newRep,_options['nocheck']);
    end if;
    resolve(s);
    nreps:=addReps(getNewSols());
    printf("新产生代表元:");
    printReps(nreps);
    return nreps;
end proc:

# 获取新的代表元
fetchNewRep:=proc(rep,solInd,con)
    return fetchSolRep(rep:-osol[solInd],addcon=con);
end proc:

# 检查两个代表元之间能否相互转化
canTransform:=proc(r1,i1,r2,i2)
    local s1,s2,eqs;
    s1:=r1:-osol[i1];
    s2:=r2:-osol[i2];
    eqs:=table();
    eqs[1]:=testTransform(s1,s2,s1);
    if (eqs[1]=true) then
        return true;
    end if;
    eqs[2]:=testTransform(s1,s2,s2);
    if (eqs[2]=true) then
        return true;
    end if;
    eqs[3]:=testTransform(s2,s1,s1);
    if (eqs[3]=true) then
        return true;
    end if;
    eqs[4]:=testTransform(s2,s1,s2);
    if (eqs[4]=true) then
        return true;
    end if;
    printf("无法求解变换方程");
    print~([entries(eqs,nolist)]);
    return false;
end proc:

# 检查两个解之间能否相互转化
# TODO　这里没有考虑转化是否有成立的条件
testTransform:=proc(s1,s2,base)
    local eq,sol,con;
    eq,sol,con:=solveTeq(s1:-rvec,s2:-rvec,base);
    if (sol<>[]) then
        print(s2:-rep);
        printf("可以在");
        print(base:-isol);
        printf("下，通过");
        print(sol);
        printf("转化为");
        print(s1:-rep);
        return true;
    else
        return eq;
    end if;
end proc:

# 输出所有InvSol对象
printSols:=proc(sols)
    local n,i;
    n:=numelems(sols);
    for i from 1 to n do
        printf("---------------------------------------------------------");
        printf("sols[%d]",i);
        printSol(sols[i]);
    end do;
    return sols;
end proc:

# 输出所有RepSol对象
printReps:=proc(reps)
    local i,n;
    n:=numelems(reps);
    for i from 1 to n do
        printf("代表元 [%d]---------------------------",i);
        printRep(reps[i]);
    end do;
end proc:

# 简要输出代表元及其成立条件
printRepCon:=proc()
    map(x->print([x:-rep,getCon(x)[x:-sid]]),getReps()):
    return;
end proc:

# 简要输出代表元及其成立条件以及不变量方程和变换方程的解
summaryReps:=proc()
    local i,n,r,_reps,id;
    _reps:=getReps();
    n:=numelems(_reps);
    for i from 1 to n do
        r:=_reps[i];
        printf("代表元 [%d]",i);
        print(r:-rep);
        printf("具有条件:");
        for id in r:-sid do
            print(getCon(r)[id]);
            print(r:-isol[id]);
            print(r:-tsol[id]);
            printf("-------------------------------------");
        end do;
    end do;
    return;
end proc:

$endif


text:
with(DifferentialThomas):
with(Tools):
Ranking([t], [[S, A, E, C, Ninv, J, I], [y2, y] ]):
sys := [
diff(S(t), t) - (-b*S(t)*Ninv(t)*A(t)*q - b*S(t)*Ninv(t)*II(t) - b*S(t)*Ninv(t)*J(t)),
diff(A(t), t) - (-E(t)*k*r + E(t)*k - A(t)*g1),
diff(E(t), t) - (b*S(t)*Ninv(t)*A(t)*q + b*S(t)*Ninv(t)*II(t) + b*S(t)*Ninv(t)*J(t) - E(t)*k),
diff(C(t), t) - (alpha*II(t)),
diff(Ninv(t), t) - (0),
diff(J(t), t) - (alpha*II(t) - g2*J(t)),
diff(II(t), t) - (-alpha*II(t) + E(t)*k*r - g1*II(t)),
y2(t) - (Ninv(t)),
y(t) - (C(t))
];
res := CodeTools[CPUTime](ThomasDecomposition(sys));

ext: mpl


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

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

# Minimum and maximum functions that I can differentiate
m_min := (x1, x2) -> x1 + (x2 - x1)*Heaviside(x1 - x2):
m_max := (x1, x2) -> x1 + (x2 - x1)*Heaviside(x2 - x1):
m_abs := (x) -> m_max(x, -x):

# a series of useful definitions

M_C         := 137.0359996287515: (* speed of light *)

RS_FACTOR   := (3/(4*Pi))^(1/3):
X2S         := 1/(2*(6*Pi^2)^(1/3)):
X2S_2D      := 1/(2*(4*Pi)^(1/2)):

X_FACTOR_C    := 3/8*(3/Pi)^(1/3)*4^(2/3):
X_FACTOR_2D_C := 8/(3*sqrt(Pi)):
K_FACTOR_C    := 3/10*(6*Pi^2)^(2/3):

MU_GE       := 10/81:
MU_PBE      := 0.06672455060314922*(Pi^2)/3:
KAPPA_PBE   := 0.8040:

# generic conversion functions
n_total    := rs -> (RS_FACTOR/rs)^3:
n_spin     := (rs, z) -> (1 + z)*n_total(rs)/2:
sigma_spin := (rs, z, xs) -> xs^2*n_spin(rs, z)^(8/3):
t_total    := (z, ts0, ts1) ->
  (ts0*((1 + z)/2)^(5/3) + ts1*((1 - z)/2)^(5/3)):
u_total    := (z, us0, us1) -> t_total(z, us0, us1):

# useful formulas that enter several functionals follow

# von Weizsäcker term
t_vw := (z, xt, us0, us1) -> (xt^2 - u_total(z, us0, us1))/8:

# See Eq. (9) of Perdew1992_13244
f_zeta    := z -> ((1 + z)^(4/3) + (1 - z)^(4/3) - 2)/(2^(4/3) - 2):
f_zeta_2d := z -> 1/2*((1 + z)^(3/2) + (1 - z)^(3/2)):

# used in several correlation functionals
mphi := z -> ((1 + z)^(2/3) + (1 - z)^(2/3))/2:
tt   := (rs, z, xt) -> xt/(4*2^(1/3)*mphi(z)*sqrt(rs)):

# in the paper it is beta_a = 0.066725
beta_Hu_Langreth := rs -> 0.066724550603149220*(1 + 0.1*rs)/(1 + 0.1778*rs):

# This is the Stoll decomposition in our language
lda_stoll_par  := (lda_func, rs, z, spin) ->
  lda_func(rs*(2/(1 + z))^(1/3), spin)*(1 + z)/2:

lda_stoll_perp := (lda_func, rs, z) ->
  + lda_func(rs, z)
  - lda_stoll_par(lda_func, rs,  z,  1)
  - lda_stoll_par(lda_func, rs, -z, -1):

gga_stoll_par  := (gga_func, rs, z, xs, spin) ->
 + gga_func(rs*(2/(1 + z))^(1/3), spin, xs, xs*(1 + spin)/2, xs*(1 - spin)/2)
 * (1 + z)/2:

# Curvature of the Fermi hole (without the current term)
Fermi_D := (xs, ts) -> 1 - xs^2/(8*ts):

# correction to Fermi_D in JCP 127, 214103 (2007); doi: http://dx.doi.org/10.1063/1.2800011
Fermi_D_a := 1e-4:
Fermi_D_corrected := (xs, ts) -> (1 - xs^2/(8*ts)) * (1 - exp(-ts^2/a^2)):

# Becke function used in several correlation functionals
b88_R_F := (f_x, rs, z, xs) ->
  1/(2*X_FACTOR_C*n_spin(rs, z)^(1/3)*f_x(xs)):

b88_zss := (css, f_x, rs, z, xs) ->
  2*css*b88_R_F(f_x, rs, z, xs):
b88_zab := (cab, f_x, rs, z, xs0, xs1) ->
  cab*(b88_R_F(f_x, rs, z, xs0) + b88_R_F(f_x, rs, -z, xs1)):

# Power series often used in mggas
mgga_w := t -> (K_FACTOR_C - t)/(K_FACTOR_C + t):
mgga_series_w := (a, n, t) -> add(a[i]*mgga_w(t)^(i-1), i=1..n):

# Used in screened functionals
kF := (rs, z) -> (3*Pi^2*(1 + z))^(1/3) * RS_FACTOR/rs:
nu := (rs, z) -> params_a_omega/kF(rs, z):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
# A striped preorder is one for which comparability is an equivalence relation

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

`is_element/striped_preord` := (A::set) -> proc(R)
 local E;
 global reason;

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

 E := R union `op/autorel`(A)(R);

 if not(`is_transitive/autorel`(A)(E)) then
  reason := [convert(procname,string),"R u R^op is not transitive",R];
  return false;
 fi;

 return true;
end;

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

`is_equal/striped_preord`     := eval(`is_equal/autorel`);
`is_leq/striped_preord`       := eval(`is_leq/preord`);
`is_separated/striped_preord` := eval(`is_separated/autorel`);

`op/striped_preord` := eval(`op/autorel`);

######################################################################
# A075729

`count_elements/striped_preord` := proc(A::set)
 local n,x,egf;
 n := nops(A);
 egf := convert(series(exp(1/(2-exp(x))-1),x=0,n+1),polynom,x):
 return n! * coeff(egf,x,n);
end:

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

`list_elements/striped_preord` := proc(A::set)
 local RR,PI,pi,LL,QQ,B,Q,L;

 RR := [];

 PI := `list_elements/partitions`(A);
 for pi in PI do
  LL := [{}];
  for B in pi do
   QQ := `list_elements/total_preord`(B);
   LL := [seq(seq(L union Q,Q in QQ),L in LL)];
  od;
  RR := [op(RR),op(LL)];
 od:

 return RR;
end:

######################################################################
# Note that the function below is not actually a rank function, but
# it is an ingredient in certain rank functions to be defined
# elsewhere.

`rank/preord` := (A::set) -> proc(R)
 nops(`block_partition/preord`(A)(R));
end:

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

`random_element/striped_preord` := (A::set) -> proc()
 local pi,Q,B;

 pi := `random_element/partitions`(A)();
 Q := {};
 for B in pi do
  Q := Q union `random_element/total_preord`(B)();
 od;
 return Q;
end:

######################################################################
# Returns a striped preorder R0 > R with
# (or FAIL if R is already maximal).

`bump/striped_preord` := (A::set) -> proc(R)
 local pi,B,a,b,B0,R0;

 pi := `block_partition/preord`(A)(R);

 for B in pi do
  if nops(B) > 1 then
   a := B[1];
   B0 := B minus {a};
   R0 := R minus {seq([b,a],b in B0)};
   return R0;
  fi;
 od;

 return FAIL;
end:

######################################################################
# Returns a striped preorder R1 with R < R0 <= S
# (or FAIL if this is impossible)

`relative_bump/striped_preord` := (A::set) -> proc(R,S)
 local pi,B,BB,SB,C,a,b,c,B0,R0;

 pi := `block_partition/preord`(A)(R);

 for B in pi do
  BB := `top/autorel`(B);
  SB := S intersect BB;
  if BB minus S <> {} then
   for a in B do
    C := select(b -> member([b,a],SB) and not(member([a,b],SB)),B);
    if C = {} then
     B0 := B minus {a};
     R0 := R minus {seq([b,a],b in B0)};
     return R0;
    fi;
   od;
   return FAIL; # should not happen
  fi;
 od;

 return FAIL;
end:

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

`describe/striped_preord` := (A::set) -> proc(R)
 local s0,s1,s2,pi,B,C,r,m,j,k,c;

 s0 := "";
 pi := `block_partition/preord`(A)(R union `op/autorel`(A)(R));
 for B in pi do
  r := `rank_table/preord`(B)(R intersect `top/autorel`(B));
  m := max(map(b -> r[b],B));
  s1 := "";
  for j from 0 to m do
   C := select(b -> r[b]=j,B);
   s2 := "";
   for k from 1 to nops(C) do
    s2 := cat(s2,`if`(k>1,"=",""),sprintf("%A",C[k]));
   od:
   s1 := cat(s1,`if`(j>0,"<",""),s2);
  od;
  s0 := cat(s0,`if`(s0="","",","),s1);
 od:

 return s0;
end:


text:
######################################################################
# trees(A) is the set of trees on A
# (Section sec-trees)

`is_element/trees` := (A::set) -> proc(TT)
 local n,i,j,T,U,is_min,M;
 global reason;

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

 n := nops(TT);
 is_min := table;
 for i from 1 to n do
  is_min[i] := true;
  T := TT[i];
  if not type(T,set) and T <> {} and T minus A = {} then
   reason := [convert(procname,string),"T is not a nonempty subset of A",T,A];
   return false;
  fi;
 od;

 # Now check nesting, and along the way, check which sets are minimal.
 for i from 1 to n do
  T := TT[i];
  for j from 1 to i-1 do
   U := TT[j];
   if T intersect U <> {} then
    if T minus U = {} then
     is_min[j] := false;
    else
     if U minus T = {} then
      is_min[i] := false;
     else
      reason := [convert(procname,string),"T and U are neither disjoint nor nested",T,U];
      return false;
     fi;
    fi;
   fi;
  od;
 od;

 # Check that the union of the minimal sets is A.
 M := NULL;
 for i from 1 to n do
  if is_min[i] then M := M,op(TT[i]); fi;
 od;
 if A minus {M} <> {} then
  reason := [convert(procname,string),"The minimal sets do not cover A",M,A];
  return false;
 fi;

 return true;
end;

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

`is_equal/trees` := (A::set) -> (TT,UU) -> evalb(TT = UU);

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

`is_leq/trees` := (A::set) -> (TT,UU) -> evalb(TT minus UU = {});

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

`is_separated/trees` := (A::set) -> proc(TT)
 local a;
 return evalb({seq({a},a in A)} minus TT = {});
end:

######################################################################
# Suppose we have a set A, a tree TT on A, a set T in TT, and a point a that
# does not lie in A.  We can then construct a new tree on A u {a} by attaching
# a to the node T.  There are two slightly different ways to do this.
# We can either attach a to a new node just above T, or we can attach a to
# the node T itself.  (However, the second version is invalid if |T|=1).

`attach_to_tree` := (A::set) -> proc(TT,T,a,above_)
 local above,CT,SCT,NCT;

 above := `if`(nargs > 3, above_, false);

 if not(`is_element/trees`(A)(TT) and
        member(T,TT) and
        not(member(a,A))) then
  return FAIL;
 fi;

 if nops(T)=1 and not(above) then
  return FAIL;
 fi;

 # Sets that contain (or are equal to) T
 CT := select(U -> (T minus U = {}),TT);

 SCT := CT minus {T};

 # Sets that do not contain T
 NCT := TT minus CT;

 if above then
  return {{a},op(NCT),T,op(map(U -> U union {a},CT))};
 else
  return {{a},op(NCT),op(map(U -> U union {a},CT))};
 fi;
end;

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

`random_element/trees`:= (A::set) -> proc()
 local n,a,A1,T0,T1,T,above;

 # Degenerate cases for small A
 if A={} then return {}; fi;
 if nops(A)=1 then return {A}; fi;

 # If A is not too small, we write it as {a} union A1.  Then we recursively
 # choose a random tree on A1, and attach a to it.

 n := nops(A);
 a := A[rand(1..n)()];
 A1 := A minus {a};
 if n = 2 then
  above := true;
 else
  above := evalb(rand(0..1)() = 1);
 fi;

 T0 := `random_element/trees`(A1)();
 T1 := select(U -> nops(U) > 1,T0);

 if above then
  T := T0[rand(1..nops(T0))()];
 else
  T := T1[rand(1..nops(T1))()];
 fi;

 T := `attach_to_tree`(A1)(T0,T,a,above);

 return T;
end;

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

`big_sets/trees` := (TT) -> select(T -> nops(T) > 1,TT);

######################################################################
# Given a tree TT, this returns a table P indexed by the non-maximal
# elements of TT, such that P[T] is the parent of T in TT.

`parent_map` := (A::set) -> proc(TT)
 local TS,n,i,j,T,U,parent;

 # if not(`is_element/trees`(A)(TT)) then return FAIL; fi;

 TS := sort([op(TT)],(a,b) -> (nops(a) < nops(b)));
 n := nops(TS);
 parent := table();

 for i from 1 to n do
  T := TS[i];
  parent[T] := FAIL;
  for j from i+1 to n do
   U := TS[j];
   if T minus U = {} then
    parent[T] := U;
    break;
   fi;
  od;
 od;

 return eval(parent);
end;

######################################################################
# Given a tree TT, this returns a table C indexed by the
# elements of TT, such that C[T] is the set of children of T in
# TT.  In particular, C[T] = {} if T is minimal.

`children_map` := (A::set) -> proc(TT)
 local parent,children,TT0,T,U;

 parent := eval(`parent_map`(A)(TT));

 children := table();
 for T in TT do
  children[T] := {};
 od:

 TT0 := TT minus {A};

 for T in TT0 do
  U := parent[T];
  if U <> FAIL then
   children[U] := {op(children[U]),T};
  fi;
 od;

 return eval(children);
end;

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

corolla := proc(A::set)
 local a;
 return {A,seq({a},a in A)};
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 *)
(* prefix:
  mgga_c_vsxc_params *params;

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

$define lda_c_pw_params
$define lda_c_pw_modified_params
$include "lda_c_pw.mpl"

$include "gvt4.mpl"

vsxc_comp := (rs, z, spin, xs, ts) ->
  + lda_stoll_par(f_pw, rs,  z,  1)
  * gtv4(params_a_alpha_ss, params_a_dss, xs, 2*(ts - K_FACTOR_C))
  * Fermi_D(xs, ts):

(* The parallel and perpendicular components of the energy *)
vsxc_fpar  := (rs, z, xs0, xs1, ts0, ts1) ->
  + vsxc_comp(rs,  z,  1, xs0, ts0)
  + vsxc_comp(rs, -z, -1, xs1, ts1):

vsxc_fperp := (rs, z, xs0, xs1, ts0, ts1) ->
  + lda_stoll_perp(f_pw, rs,  z)
  * gtv4(params_a_alpha_ab, params_a_dab, sqrt(xs0^2 + xs1^2), 2*(ts0 + ts1 - 2*K_FACTOR_C)):

vsxc_f := (rs, z, xs0, xs1, ts0, ts1) ->
  + vsxc_fpar (rs, z, xs0, xs1, ts0, ts1)
  + vsxc_fperp(rs, z, xs0, xs1, ts0, ts1):

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


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

`eta/cubes` := (k::posint) -> (A::set) ->
 `if`(nops(A) = 1,table(op(A) = `id/single_cubes`(k)),FAIL);

`gamma/cubes` := (k::posint) -> (A::set,B::set) -> (p) -> proc(f,g)
 local a,h;

 h := table();
 for a in A do
  h[a] := `compose/single_cubes`(k)(f[p[a]],g[p[a]][a]);
 od;
 return eval(h);
end;


text:
kernelopts(assertlevel=2): # be strict on all assertions while testing
kernelopts(opaquemodules=false): # allow testing of internal routines
if not (NewSLO :: `module`) then
  WARNING("loading NewSLO failed");
  `quit`(3);
end if;

with(Hakaru):
with(NewSLO):
with(Partition):

#####################################################################
#
# disintegration tests
#
#####################################################################
d1 := Bind(Lebesgue(-infinity,infinity), x, Ret(Pair(-5*x,3/x))):
d1r := {Weight(1/5,Ret(-15/t))}:

# should try d2 with 1/c as well
d2 := Bind(Lebesgue(-infinity,infinity), x, Ret(Pair((-1/7)*x-1,3))):
d2r := {Weight(7, Ret(3))}:

#The next two tests are simplified versions of the Borel-Kolmogorov paradox.
#https://en.wikipedia.org/wiki/Borel-Kolmogorov_paradox
d3 := Bind(Uniform(0,1), x, Bind(Uniform(0,1), y, Ret(Pair(x-y,f(x,y))))):
d3r := {
  PARTITION([Piece(Or(t <= -1,1 < t)
                  ,Msum())
            ,Piece(And(t <= 0, -1 < t)
                  ,Weight(t+1,Bind(Uniform(-t,1),`y`,Ret(f(t+y,y)))))
            ,Piece(And(t <= 1, 0 < t)
                  ,Weight(1-t,Bind(Uniform(0,1-t),`y`,Ret(f(t+y,y)))))
            ]) ,
  PARTITION([Piece(Or(t <= -1,1 < t)
                  ,Msum())
            ,Piece(And(t <= 0, -1 < t)
                  ,Weight(t+1,Bind(Uniform(0,t+1),`x`,Ret(f(x,-t+x)))))
            ,Piece(And(t <= 1, 0 < t)
                  ,Weight(1-t,Bind(Uniform(t,1),`x`,Ret(f(x,-t+x)))))
            ]),
PARTITION([Piece(Or(t <= -1,1 <= t),Msum())
          ,Piece(And(-1 < t,t <= 0),Weight(t+1,Bind(Uniform(0,t+1),x,Ret(f(x,-t+x)))))
          ,Piece(And(0 < t,t < 1),Weight(1-t,Bind(Uniform(t,1),x,Ret(f(x,-t+x)))))
          ])
}:

BUniform := proc(x,b,$) Bind(Uniform(0,1), x, b) end proc:
d3_3 := BUniform(x, BUniform(y, BUniform(z, Ret(Pair(x+y+z,f(x,y,z)))))):
d3_3_r := {}:

d4 := Bind(Uniform(0,1), x, Bind(Uniform(0,1), y, Ret(Pair(x/y,x)))):
d4r := {
  Weight(1/abs(t)^2,
    Bind(Uniform(0,1),x,
         piecewise(x < t,Weight(x,Ret(x)),Msum()))),
  piecewise(0 < t,
    Bind(Uniform(0,1),y,
         piecewise(t < 1/y,
           Weight(y,Ret(t*y)),
           Msum())),
    Msum()),
  PARTITION([Piece(t <= 0,Msum())
           , Piece(And(t <= 1, 0 < t),Weight(1/t,Bind(Uniform(0,t),x,Weight(x,Ret(x)))))
           , Piece(1 < t,Weight(1/2/t^2,BetaD(2,1)))]),
  PARTITION([Piece(t <= 0,Msum()), Piece(0 < t and t < 1,Weight(1/t,Bind(Uniform(0,t),x,Weight(x,Ret(x))))), Piece(1 <= t,Weight(1/2/t^2,BetaD(2,1)))])
}:

# like d3 but positive, and the entire parametric family.
d3posfam := Bind(Uniform(0, 1), x, Bind(Uniform(0, 1), y, Ret(Pair(x+y+K_0, f(x, y))))):
d3posfam_r := {
  PARTITION([ Piece(t <= K_0, Msum())
            , Piece(And(K_0 < t, t <= 1+K_0), Weight(t-K_0, Bind(Uniform(0, t-K_0), x, Ret(f(x, t-x-K_0)))))
            , Piece(And(t < 2+K_0, 1+K_0 < t), Weight(2-t+K_0, Bind(Uniform(-1+t-K_0, 1), x, Ret(f(x, t-x-K_0)))))
            , Piece(2+K_0 <= t, Msum())]) }:
d3posfam_ctx := [K_0::real]:

d5 := Bind(Gaussian(0,1), x, Bind(Gaussian(x,1), y, Ret(Pair(y,x)))):
d5r := {Weight((1/2)*exp(-(1/4)*t^2)/Pi^(1/2), Gaussian((1/2)*t, (1/2)*2^(1/2)))}:

d6 := Bind(Gaussian(0,1), x, Bind(Gaussian(x,1), y, Ret(Pair(x,y)))):
d6r := {Weight(1/2*2^(1/2)/Pi^(1/2)*exp(-1/2*t^2),Gaussian(t,1))}:

# note (y+y), which gives trouble for a syntactic approach
normalFB1 :=
  Bind(Gaussian(0,1), x,
  Bind(Gaussian(x,1), y,
  Ret(Pair((y+y)+x, _Unit)))):

normalFB1r := {Weight(1/26*exp(-1/26*t^2)/Pi^(1/2)*13^(1/2)*2^(1/2),Ret(_Unit))}:

# tests taken from haskell/Tests/Disintegrate.hs
# use same names, to be clearer
norm0a :=
  Bind(Gaussian(0,1), x,
  Bind(Gaussian(x,1), y,
  Ret(Pair(y, x)))):
  # note that the answer below is much nicer than the one expected in Haskell
norm0r := {Weight(1/2*exp(-1/4*t^2)/Pi^(1/2),Gaussian(t/2,sqrt(2)/2))}:

norm1a :=
  Bind(Gaussian(3,2), x,Ret(piecewise(x<0, Pair(-x, _Unit), Pair(x, _Unit)))):
norm1b :=
  Bind(Gaussian(3,2), x,piecewise(x<0, Ret(Pair(-x, _Unit)), Ret(Pair(x, _Unit)))):

norm1r_w := (1/4)*sqrt(2)*exp(-9/8)*((exp(t))^(3/4)/(exp(t^2))^(1/8)+1/((exp(t^2))^(1/8)*(exp(t))^(3/4)))/sqrt(Pi):
norm1r := {
  Weight( piecewise(t < 0  , 0
                   ,0 <= t , norm1r_w
                   )
        , Ret(_Unit)
        ) ,
  Weight( PARTITION([Piece(t < 0  , 0)
                    ,Piece(0 <= t , norm1r_w)
                   ])
        , Ret(_Unit)
        )
}:

assume(s::real, noiseT >= 3, noiseT <= 8, noiseE >= 1, noiseE <= 8);
easyRoad:= [
  Bind(Uniform(3, 8), noiseT,
  Bind(Uniform(1, 4), noiseE,
  Bind(Gaussian(0, noiseT), x1,
  Bind(Gaussian(x1, noiseE), m1,
  Bind(Gaussian(x1, noiseT), x2,
  Bind(Gaussian(x2, noiseE), m2,
  Ret(Pair(Pair(m1,m2), Pair(noiseT,noiseE)))
  )))))),
  Pair(s,t)
]:
#The first expression below comes from the actual output of disint, hand-
#simplified 1) to bring factors into the innnermost integral, 2) to combine
#products of exps, and 3) to express the polynomial arg of exp in a logical way
#by sub-factoring.
easyRoadr:= {
  Weight(                #Weight 1
    Pi/8,
    Bind(                #Bind 1
      Uniform(3, 8), noiseT,
      Weight(            #Weight 2
        1/noiseT^2,
        Bind(            #Bind 2
          Uniform(1, 4), noiseE,
          Weight(        #Weight 3
            int(         #Int 1
              int(       #Int 2
                exp(
                  -(x2^2/2 - x1*x2 + x1^2)/noiseT^2 -
                  ((t-x2)^2 + (s-x1)^2)/noiseE^2
                )*2/Pi/noiseE,
                x2= -infinity..infinity
              ),         #-Int 2
              x1= -infinity..infinity
            ),           #-Int 1
            Ret(Pair(noiseT, noiseE))
          )              #-Weight 3
        )                #-Bind 2
      )                  #-Weight 2
    )                    #-Bind 1
  ),                     #-Weight 1

  #Hopefully, that's equivalent to...
  Bind(Uniform(3, 8), noiseT,
  Bind(Uniform(1, 4), noiseE,
  Bind(Gaussian(0, noiseT), x1,
  Bind(Weight(density[Gaussian](x1, noiseE)(s), Ret(_Unit)), _,
  Bind(Gaussian(x1, noiseT), x2,
  Weight(density[Gaussian](x2, noiseE)(t), Ret(Pair(noiseT, noiseE)))
  )))))
}:
helloWorld:=
  Bind(Gaussian(0,1), mu,
  Bind(Plate(n, k, Gaussian(mu, 1)), nu,
  Ret(Pair(nu, mu))
  )):
helloWorldr:= {
  Bind(Gaussian(0,1), mu,
  Plate(n, i, Weight(density[Gaussian](mu, 1)(idx(t,i)), Ret(mu)))
  )
}:

pair_x_x := Bind(Uniform(0, 1), x, Ret(Pair(x, x))):
pair_x_x_r := {
  PARTITION([Piece(t < 0, Msum()), Piece(And(0 <= t, t <= 1), Ret(t)),
             Piece(1 < t, Msum())]) }:

#This first block of tests is to test the basic functionality of disint, and,
#to some extent, the system as a whole. These tests may be meaningless to you,
#the statistician and end user of this Hakaru product; they aren't meant to
#have any statistical meaning.--Carl 2016Oct04

TestDisint(
     [Ret(Pair(sqrt(Pi), x)), t &M Ret(7)],
     {Msum()},
     label= "(d0_2) `Dirac` test 1"
);

TestDisint(
     [Ret(Pair(sqrt(Pi), x^2)), t &M Ret(sqrt(Pi))],
     {Ret(x^2)},
     label= "(d0_3) `Dirac` test 2"
);

TestDisint(
     [Bind(Lebesgue((-1,1)*~infinity), x, Ret(Pair(sqrt(Pi), x^2))),
      t &M Ret(sqrt(Pi))
     ],
     {Bind(Lebesgue((-1,1)*~infinity), x1, Ret(x1^2))},
     label= "(d0_4) `Dirac` test with `Bind`"
);

TestDisint(d1, d1r, label = "(d1) Disintegrate linear function");
TestDisint(d2, d2r, label = "(d2) Disintegrate linear function II");
TestDisint(d5, d5r, label = "(d5) Disintegrate N(0,1)*N(x,1), over y");
TestDisint(d6, d6r, label = "(d6) Disintegrate N(0,1)*N(x,1), over x");
TestDisint(norm0a, norm0r,
     label = "(norm0a) U(0,1) >>= \x -> U(x,1) >>= \y -> Ret(y,x)"
);

## This one is kind of cosmetic; it would be 'fixed' properly if the
## disintegration process did not use 'improve' to do "domain information
## discovery", but rather had a specific function (and then improve could
## indeed do this integral).
# should work now
TestDisint( normalFB1, normalFB1r,
     label = "(d7_normalFB1) Disintegrate N(0,1)*N(x,1), over (y+y)+x"
            );

TestDisint(d3, d3r, label = "(d3) Disintegrate U(0,1) twice, over x-y");

######################################################################
#
# These tests fail, and are expected to.  Move them up when they
# start passing (and are expected to).
#
# They are, however, roughly in order of what we'd like to have work.
#

# change of variables
TestDisint(d4, d4r, label = "(d4) Disintegrate U(0,1) twice, over x/y");
# funky piecewise
TestDisint(norm1a, norm1r,
     label = "(norm1a) U(0,1) into Ret of pw"
);
TestDisint(norm1b, norm1r,
     label = "(norm1b) U(0,1) into pw of Ret"
);
#In this one the function in the inequality, x+x^3, is injective but nonlinear.
TestDisint(
     Bind(Gaussian(0,1), x, Ret(Pair(x+x^3, f(x)))),
     {}, #I don't know what to expect.
     label= "(d0_5) Injective nonlinear inequality"
);

TestDisint(pair_x_x, pair_x_x_r, label="(pair_x_x) Disintegrate U(0,1) over Ret(x,x)");
TestDisint(d3_3, d3_3_r, [], 180
          , label = "(d3_3) Disintegrate U(0,1) thrice, over x+y+z");
TestDisint(d3posfam, d3posfam_r, d3posfam_ctx
          , label = "(d3posfam) Disintegrate U(0,1) twice, over x+y+K");

#This one is a basic test of the Counting wrt-var type.
#This one gives the Weight(-1, ...) error
TestDisint(
     [Bind(PoissonD(2), n, Ret(Pair(3,n))), n_wrt &M Counting((-1,1)*~infinity)],
     {},  #I don't know what to expect.
     label= "(d0_1) `Counting` test; `Weight` bug (currently failing)"
);

TestDisint(
     helloWorld, helloWorldr,
     [n::integer, n > 0],
     label= "(helloWorld) Plate of Normals"
);

# tests which take too long
TestDisint(
     easyRoad, easyRoadr, [],
     120, #takes 6 - 8 minutes to `improve` on an Intel i7
     label= "(easyRoad) Combo of Normals with distinct Uniform noises"
);

# Example from Chad (April 9, 2017)
prog :=
 Bind(Gaussian(70,30), temp,
 Bind(Gaussian(temp, 10), withErr,
 Ret(Pair(min(withErr, 100), temp)))):
TestDisint( prog, {}, label= "censored temperature example");

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
# h := discrim(a1*x^4 + x^5 + a2*x^3 + a3*x^2 + a4*x + a5, x):

# # (1)-(47) are from Wang Dongming - So1ving po1ynomia1 equations: Characteristic sets and triangu1ar systems
# # 共16个例子

wangzero :=
[
    [   #(2) wang_ex9
        [x^8 + 4*x^6*y^2 + 6*x^4*y^4 + 4*x^2*y^6 + y^8 - 2*x^6*y - 6*x^4*y^3 - 6*x^2*y^5 - 2*y^7 - 5*x^6 - 11*x^4*y^2 - 7*x^2*y^4 - y^6 + 10*x^4*y + 12*x^2*y^3 + 2*y^5 + 4*x^4 + 4*x^2*y^2 - 8*x^2*y,
        8*x^7 + 24*x^5*y^2 + 24*x^3*y^4 + 8*x*y^6 - 12*x^5*y - 24*x^3*y^3 - 12*x*y^5 - 30*x^5 - 44*x^3*y^2 - 14*x*y^4 + 40*x^3*y + 24*x*y^3 + 16*x^3 + 8*x*y^2 - 16*x*y,
        8*x^6*y + 24*x^4*y^3 + 24*x^2*y^5 + 8*y^7 - 2*x^6 - 18*x^4*y^2 - 30*x^2*y^4 - 14*y^6 - 22*x^4*y - 28*x^2*y^3 - 6*y^5 + 10*x^4 + 36*x^2*y^2 + 10*y^4 + 8*x^2*y - 8*x^2
        ],
        [x, y]
    ],
    [   #(7) wang_ex15
        [8*x^2-2*x*y-6*x*z+3*x+3*y^2-7*y*z+10*y+10*z^2-8*z-4,
         10*x^2-2*x*y+6*x*z-6*x+9*y^2-y*z-4*y-2*z^2+5*z-9,
         5*x^2+8*x*y+4*x*z+8*x+9*y^2-6*y*z+2*y-z^2-7*z+5
        ],
        [x, y, z]
    ],
    [
        #(8) wang_ex16
        [2 + 10*p^2 - 2*p*q - 2*q + c^2*(q - 1)^2 + 2*c*d*(1 - q)*(q - p) + 4*p*q*d*(d - c) + 4*d^2*(1 - 2*p) + 4*d^2*(p - q) + 4*d*c*(p - 1) + 4*p*q*(c + 1) - 8*p + d^2*(p - q)^2,
        d*(2*p + 1)*(q - p) + c*(p + 2)*(1 - q) - 6*p*d + 2*c*(-p*q + p + q - 1) + 6*p^2*d,
        -4*(p - 1)^2 + 2*p*(p - q) - 2*q + 2,
        4 - 4*q^2 + 4*p + 3*c^2*(q - 1)^2 - 3*d^2*(p - q)^2 + 12*d^2*(p - 1)^2 + 4*p*(p - 2) + 12*d*c*(p + q + p*q - 1)
        ],
        [q, c, p, d]
    ],
    [
        #(10) wang_ex21
        [x + y + z + t + u,
         x*y + y*z + z*t + t*u + u*x,
         x*y*z + y*z*t + z*t*u + t*u*x + u*x*y,
         x*y*z*t + y*z*t*u + z*t*u*x + t*u*x*y + u*x*y*z,
         x*y*z*t*u + 1
        ],
        [x, y, z, t, u]
    ],
    [   #(13) wang_ex26
        [x^2 + y^2 + z^2,
         z*y*x,
         y^2*x^2 + z^2*x^2 + z^2*y^2,
         u^2 + 1/3*t^2,
         u^3 - t^2*u,
         2*x^2*u - y^2*u - z^2*u + y^2*t - z^2*t,
         -y^2*x^2*u - z^2*x^2*u + 2*z^2*y^2*u - y^2*x^2*t + z^2*x^2*t,
         2*x^2*u^2 - y^2*u^2 - z^2*u^2 - 2*y^2*t*u + 2*z^2*t*u - 2/3*x^2*t^2 + 1/3*y^2*t^2 + 1/3*z^2*t^2,
         - y^2*x^2*u^2 - z^2*x^2*u^2 + 2*z^2*y^2*u^2 + 2*y^2*x^2*t*u - 2*z^2*x^2*t*u + 1/3*y^2*x^2*t^2 + 1/3*z^2*x^2*t^2 - 2/3*z^2*y^2*t^2,
         - 3*y^2*x^4*t*u^2 + 3*z^2*x^4*t*u^2 + 3*y^4*x^2*t*u^2 - 3*z^4*x^2*t*u^2 - 3*z^2*y^4*t*u^2 + 3*z^4*y^2*t*u^2
         + 1/3*y^2*x^4*t^3 - 1/3*z^2*x^4*t^3 - 1/3*y^4*x^2*t^3 + 1/3*z^4*x^2*t^3 + 1/3*z^2*y^4*t^3 - 1/3*z^4*y^2*t^3
        ],
        [x, y, z, t, u]
    ],
    [   #(14) wang_ex27
        [y + u + v - 1,
         z + t + 2*u - 3,
         y + t + 2*v - 1,
         x - y - z - t - u - v,
         t*u*x^2 - 1569/31250*y*z^3,
         z*v - 587/15625*y*t
        ],
        [x, y, z, t, u, v]
    ],
    [
        #(15) wang_ex28
        [x^2 - x + 2*y^2 + 2*z^2 + 2*t^2,
         2*x*y + 2*y*z + 2*z*t - y,
         2*x*z + y^2 + 2*y*t - z,
         x + 2*y + 2*z + 2*t - 1
        ],
        [x, y, z, t]
    ],
    [   #(21) wang_ex34
        [B2*C2 + B3*C3 + B4*C4 + B5*C5 - 1/2,
         B2*C2^2 + B3*C3^2 + B4*C4^2 + B5*C5^2 - 1/3,
         B3*A32*C2 + B4*A42*C2 + B4*A43*C3 + B5*A52*C2 + B5*A53*C3 + B5*A54*C4 - 1/6,
         B2*C2^3 + B3*C3^3 + B4*C4^3 + B5*C5^3 - 1/4,
         B3*C3*A32*C2 + B4*C4*A42*C2 + B4*C4*A43*C3 + B5*C5*A52*C2 + B5*C5*A53*C3 + B5*C5*A54*C4 - 1/8,
         B3*A32*C2^2 + B4*A42*C2^2 + B4*A43*C3^2 + B5*A52*C2^2 + B5*A53*C3^2 + B5*A54*C4^2 - 1/12,
         B4*A43*A32*C2 + B5*A53*A32*C2 + B5*A54*A42*C2 + B5*A54*A43*C3-1/24,
         B2*C2^4 + B3*C3^4 + B4*C4^4 + B5*C5^4 - 1/5,
         B3*C3^2*A32*C2 + B4*C4^2*A42*C2 + B4*C4^2*A43*C3 + B5*C5^2*A52*C2 + B5*C5^2*A53*C3 + B5*C5^2*A54*C4 - 1/10,
         B3*C2^2*A32*C3 + B4*C2^2*A42*C4 + B4*C3^2*A43*C4 + B5*C2^2*A52*C2 + B5*C3^2*A53*C5 + B5*C4^2*A54*C5 - 1/15,
         B4*C4*A43*A32*C2 + B5*C5*A53*A32*C2 + B5*C5*A54*A42*C2
         + B5*C5*A54*A43*C3 - 1/30,
         B3*A32^2*C2^2 + B4*A42^2*C2^2 + 2*B4*A42*C2*A43*C3 + B4*A43^2*C3^2 + B5*A52^2*C2^2 + B5*A53^2*C3^2 + B5*A54^2*C4^2 + 2*B5*A52*C2*A53*C3 + 2*B5*A52*C2*A54*C4 + 2*B5*A53*C3*A54*C4 - 1/20,
         B3*A32*C2^3 + B4*A42*C2^3 + B4*A43*C3^3 + B5*A52*C2^3 + B5*A53*C3^3 + B5*A54*C4^3 - 1/20,
         B4*A43*C3*A32*C2 + B5*A53*C3*A32*C2 + B5*A54*C4*A42*C2 + B5*A54*C4*A43*C3- 1/40,
         B4*A43*A32*C2^2 + B5*A53*A32*C2^2 + B5*A54*A42*C2^2 + B5*A54*A43*C3^2- 1/60,
         B5*A54*A43*A32*C2-1/120
        ],
        [B2, B3, B4, B5, A32, A42, A43, A52, A53, A54, C4, C2, C3, C5]
    ],
    [   #(23) wang_ex36
        [U0^2 - U0 + 2*U1^2,
         U0 + 2*U1 - 1
        ],
        [U0, U1]
    ],
    [   #(24) wang_ex37
        [U0^2 - U0 + 2*U1^2 + 2*U2^2,
         2*U0*U1 + 2*U1*U2 - U1,
         U0 + 2*U1 + 2*U2 - 1
        ],
        [U0, U1, U2]
    ],
    [   #(25) wang_ex38
        [U0^2 - U0 + 2*U1^2 + 2*U2^2 + 2*U3^2,
         2*U0*U1 + 2*U1*U2 + 2*U2*U3 - U1,
         2*U0*U2 + U1^2 + 2*U1*U3 - U2,
         U0 + 2*U1 + 2*U2 + 2*U3 - 1
        ],
        [U0, U1, U2, U3]
    ],
    [   #(26) wang_ex39
        [
         U0^2 - U0 + 2*U1^2 + 2*U2^2 + 2*U3^2 + 2*U4^2,
         2*U0*U1 + 2*U1*U2 + 2*U2*U3 + 2*U3*U4 - U1,
         2*U0*U2 + U1^2 + 2*U1*U3 + 2*U2*U4 - U2,
         2*U0*U3 + 2*U1*U2 + 2*U1*U4 - U3,
         U0 + 2*U1 + 2*U2 + 2*U3 + 2*U4 - 1
        ],
        [U4, U3, U2, U1, U0]
    ],
    [   #(27) wang_ex40
        [U0^2 - U0 + 2*U1^2 + 2*U2^2 + 2*U3^2 + 2*U4^2 + 2*U5^2,
         2*U0*U1 + 2*U1*U2 + 2*U2*U3 + 2*U3*U4 + 2*U4*U5 - U1,
         2*U0*U2 + U1^2 + 2*U1*U3 + 2*U2*U4 + 2*U3*U5 - U2,
         2*U0*U3 + 2*U1*U2 + 2*U1*U4 + 2*U2*U5 - U3,
         2*U0*U4 + 2*U1*U3 + 2*U1*U5 + U2^2 - U4,
         U0 + 2*U1 + 2*U2 + 2*U3 + 2*U4 + 2*U5 - 1
        ],
        [U5, U4, U3, U2, U1, U0]
    ],
    [   #(28) wang_ex43
        [U4^4 - 20/7*A46^2,
         A46^2*U3^4 + 7/10*A46*U3^4 + 7/48*U3^4 - 50/27*A46^2 - 35/27*A46 - 49/216,
         A46^5*U4^3 + 7/5*A46^4*U4^3 + 609/1000*A46^3*U4^3 + 49/1250*A46^2*U4^3 - 27391/800000*A46*U4^3 - 1029/160000*U4^3 + 3/7*A46^5*U3*U4^2 + 3/5*A46^6*U3*U4^2 + 63/200*A46^3*U3*U4^2 + 147/2000*A46^2*U3*U4^2 + 4137/800000*A46*U3*U4^2 - 7/20*A46^4*U3^2*U4 - 77/125*A46^3*U3^2*U4 - 23863/60000*A46^2*U3^2*U4 - 1078/9375*A46*U3^2*U4 - 24353/1920000*U3^2*U4 - 3/20*A46^4*U3^3 - 21/100*A46^3*U3^3 - 91/800*A46^2*U3^3 - 5887/200000*A46*U3^3 - 343/128000*U3^3
        ],
        [U3, U4, A46]
    ],
    [   #(29) wang_ex45
        [45*P + 35*S - 165*B - 36,
         35*P + 40*Z + 25*T - 27*S,
         15*W + 25*P*S + 30*Z - 18*T - 165*B^2,
         -9*W + 15*P*T + 20*Z*S,
         W*P + 2*Z*T - 11*B^3,
         99*W - 11*S*B + 3*B^2
        ],
        [W, P, Z, T, S, B]
    ],
    [   #(30) wang_ex46
        [45*P + 35*S - 165*B - 36,
         35*P + 40*Z + 25*T - 27*S,
         15*W + 25*P*S + 30*Z - 18*T - 165*B^2,
         -9*W + 15*P*T + 20*Z*S,
         W*P + 2*Z*T - 11*B^3,
         99*W - 11*S*B + 3*B^2,
         B^2 + 33/50*B + 2673/10000
        ],
        [W, P, Z, T, S, B]
    ]
]:

max_stars_repo_name: lihaokun/StrongSfTriDec

max_issues_repo_name: lihaokun/StrongSfTriDec

max_forks_repo_name: lihaokun/StrongSfTriDec


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

`is_element/WK` := (A::set) -> proc(TTxl)
 local TT,x,l,TT1,T,n,u,children;
 global reason;

 if not(type(TTxl,list) and nops(TTxl) = 3) then
  reason := [convert(procname,string),"TTxl cannot be split as [TT,x,l]",TTxl];
  return false;
 fi;

 TT,x,l := op(TTxl);

 if not(`is_element/trees`(A)(TT)) then
  reason := [convert(procname,string),"TT is not a tree on A",TT,A];
  return false;
 fi;

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

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

 if map(op,{indices(x)}) <> TT1 then
  reason := [convert(procname,string),"x is not indexed by TT1",x,TT1];
  return false;
 fi;

 children := `children_map`(A)(TT);
 for T in TT1 do
  if not(`is_element/stasheff_star`(children[T])(x[T])) then
   reason := [convert(procname,string),"x[T] has invalid type",T,x[T],reason];
   return false;
  fi;
 od;

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

 if map(op,{indices(l)}) <> TT then
  reason := [convert(procname,string),"l is not indexed by TT",l,TT];
  return false;
 fi;

 n := nops(A);

 for T in TT do
  if not(`is_element/RR`(l[T]) and l[T] >= 0 and l[T] <= 1) then
   reason := [convert(procname,string),"l[T] is not in the unit interval",T,l[T]];
   return false;
  fi;

  if nops(T) = 1 and l[T] <> 1 then
   reason := [convert(procname,string),"l[{a}] <> 1",op(T),l[T]];
   return false;
  fi;

  if nops(T) = n and l[T] <> 1 then
   reason := [convert(procname,string),"l[A] <> 1",T,l[T]];
   return false;
  fi;
 od;

 return true;
end;

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

`is_interior/WK` := (A::set) -> proc(TTxl)
 local TT,x,l,n,T;
 global reason;

 TT,x,l := op(TTxl);
 n := nops(A);

 for T in TT do
  if nops(T) > 1 and nops(T) < n and l[T] = 1 then
   reason := [convert(procname,string),"not interior, l[T] = 1",T];
   return false;
  fi;
 od:

 return true;
end;

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

`is_reduced/WK` := (A::set) -> proc(TTxl)
 local TT,x,l,T;
 global reason;

 TT,x,l := op(TTxl);

 for T in TT do
  if l[T] = 0 then
   reason := [convert(procname,string),"not reduced, l[T] = 0",T];
   return false;
  fi;
 od:

 return true;
end;

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

`reduce_once/WK` := (A::set) -> proc(TTxl,T)
 local TT,x,l,parent,children,n,TT1,x1,l1,U,B,C,X,p,y,z;

 TT,x,l := op(TTxl);

 parent := `parent_map`(A)(TT);
 children := `children_map`(A)(TT);

 TT1 := TT minus {T};

 U := parent[T];
 C := children[U];
 B := C minus {T} union children[T];

 p := table();
 for X in C minus {T} do p[X] := X; od;
 for X in children[T] do p[X] := T; od;

 y := table();
 for X in C minus {T} do y[X] := `eta/stasheff_star`({X}); od;
 y[T] := x[T];

 z := `gamma/stasheff_star`(B,C)(p)(x[U],y);

 x1 := table();
 l1 := table();

 for X in TT1 do
  l1[X] := l[X];
  if nops(X) > 1 then
   x1[X] := `if`(X = U,eval(z),eval(x[X]));
  fi;
 od;

 return([TT1,eval(x1),eval(l1)]);
end;

`reduce/WK` := (A::set) -> proc(TTxl)
 local is_reduced,X,TT,x,l,T;

 is_reduced := false;
 X := TTxl;

 while not(is_reduced) do
  TT,x,l := op(X);
  is_reduced := true;

  for T in TT do
   if l[T] = 0 then
    is_reduced := false;
    X := `reduce_once/WK`(A)(X,T);
    break;
   fi;
  od:
 od;

 return X;
end;

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

`is_equal/WK` := (A::set) -> proc(X0,X1)
 local Y0,Y1,TT0,TT1,x0,x1,l0,l1,children,T;
 global reason;

 Y0 := `reduce/WK`(A)(X0);
 Y1 := `reduce/WK`(A)(X1);
 TT0,x0,l0 := op(Y0);
 TT1,x1,l1 := op(Y1);

 if TT0 <> TT1 then
  reason := [convert(procname,string),"TT0 <> TT1",TT0,TT1];
  return false;
 fi;

 children := `children_map`(A)(TT0);

 for T in TT0 do
  if l0[T] <> l1[T] then
   reason := [convert(procname,string),"l0[T] <> l1[T]",T,l0[T],l1[T]];
   return false;
  fi;

  if nops(T) > 1 then
   if not(`is_equal/stasheff_star`(children[T])(x0[T],x1[T])) then
    reason := [convert(procname,string),"x0[T] <> x1[T]",T,x0[T],x1[T],reason];
    return false;
   fi;
  fi;
 od;

 return true;
end;

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

`is_leq/WK` := NULL;

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

`random_element/WK` := (A::set) -> proc()
 local TT,x,l,TT1,T,n,children,d;

 d := 3;
 TT := `random_element/full_trees`(A)();

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

 x := table();
 l := table();

 children := `children_map`(A)(TT);
 for T in TT1 do
  x[T] := `random_element/stasheff_star`(children[T])();
 od;

 n := nops(A);

 for T in TT do
  if nops(T) = 1 or nops(T) = n then
   l[T] := 1;
  else
   l[T] := rand(0..d)()/d;
  fi;
 od;

 return [TT,eval(x),eval(l)];
end;

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

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

`eta/WK` := proc(A::set)
 if nops(A) <> 1 then return FAIL; fi;
 return [[op(A)],table([])];
end;

`gamma/WK` := (A::set,B::set) -> (p) -> proc(M,N)
 local F,TT,SS,UU,x,y,z,l,m,n,pi,T,T1,C,C1,U,b,L,children;

 F := fibres(A,B)(p);

 TT,x,l := op(M);
 children := children_map(B)(TT);
 pi := table();
 for T in TT do
  pi[T] := map(op,map(b -> F[b],T));
 od;

 SS := table();
 y := table();
 m := table();

 for b in B do
  SS[b],y[b],m[b] := op(N[b]);
 od;

 UU := {seq(op(SS[b]),b in B)};
 z := table();
 n := table();

 for T in TT do
  T1 := pi[T];
  if nops(T1) > 0 then
   UU := {op(UU),T1};
   n[T1] := l[T];
   if nops(T) = 1 then
    if nops(T1) > 1 then
     z[T1] := eval(y[op(T)][T1]);
    fi;
   else
    C := children[T];
    C1 := map(U -> pi[U],C);
    z[T1] :=  `act/stasheff_star`(C,C1)(pi)(x[T]);
   fi;
  fi;
 od;

 for b in B do
  for U in SS[b] do
   n[U] := m[b][U];
   if nops(U) > 1 then
    z[U] := y[b][U];
   fi;
  od;
 od;

 L := [UU,eval(z),eval(n)];
 L := `reduce/WK`(A)(L);
 return L;
end;

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

`phi/WK/K` := (A::set) -> proc(TTxl)
 local SS,TT,x,l,TT1,s,t,u,RR,R,T,f,parent,children,J,J1,JJ,i,j,k,n,m;

 TT,x,l := op(TTxl);

 TT1 := select(U -> nops(U) > 1,TT);
 parent := `parent_map`(A)(TT);
 children := `children_map`(A)(TT);

 s := table():
 RR := table():

 for T in TT1 do
  RR[T],s[T] := op(x[T]);
 od:

 R := `phi/WK/K/aux`(A,RR);
 u := table():
 n := nops(A);
 JJ := {seq(seq({seq(R[k],k=i..j)},j=i..n),i=1..n)};

 for J in JJ do
  if member(J,TT) then
   u[J] := (1+l[J])/2;
  else
   SS := select(T -> (J minus T = {}),TT);
   m := min(op(map(nops,SS)));
   T := select(T -> (nops(T) = m),SS)[1];
   J1 := select(U -> (U intersect J <> {}),children[T]);
   if {op(map(op,J1))} = J then
    u[J] := s[T][J1]/2;
   else
    u[J] := 0;
   fi;
  fi;
 od;

 return [R,eval(u)];
end;

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

`phi/WK/K/aux` := proc(T,UU)
  if nops(T) = 1 then
   return [op(T)];
  else
   map(op,map(U -> `phi/WK/K/aux`(U,UU),UU[T]));
  fi;
 end;

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

`theta/K/WK` := (A::set) -> proc(Rt)
 local R,t,r,T,TT,RR,s,l,x,C,m,J,J1,i,j,k,children;

 R,t := op(Rt);

 r := table();
 for i from 1 to nops(R) do
  r[R[i]] := i;
 od;

 TT := select(J -> t[J] > 1/2,map(op,{indices(t)}));
 RR := table():
 s := table():
 l := table():
 x := table():

 children := `children_map`(A)(TT);
 for T in TT do
  l[T] := 2*t[T] - 1;
  C := [op(children[T])];
  C := sort(C,(U0,U1) -> (r[U0[1]] < r[U1[1]]));
  RR[T] := C;
  m := nops(C);
  if m > 1 then
   s[T] := table():
   for i from 1 to m do
    for j from i to m do
     J1 := {seq(C[k],k=i..j)};
     J := map(op,J1);
     s[T][J1] := min(1,2*t[J]);
    od;
   od;
   x[T] := [RR[T],eval(s[T])];
  fi;
 od:

 return [TT,eval(x),eval(l)];
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/.
*)

$define lda_c_pw_params
$define lda_c_pw_modified_params
$include "lda_c_pw.mpl"

ux    := (mgamma, x) -> mgamma*x^2/(1 + mgamma*x^2):

(* article uses t = 2 tau convention *)
wx_ss := (t, dummy) -> (K_FACTOR_C - t)/(K_FACTOR_C + t):
wx_os := (ts0, ts1) -> (K_FACTOR_C*(ts0 + ts1) - 2*ts0*ts1)/(K_FACTOR_C*(ts0 + ts1) + 2*ts0*ts1):

b97mv_g := (mgamma, wx, cc, n, xs, ts0, ts1) ->
  add(cc[i][1]*wx(ts0, ts1)^cc[i][2]*ux(mgamma, xs)^cc[i][3], i=1..n):

b97mv_fpar  := (lda_func, rs, z, xs0, xs1, ts0, ts1) ->
  + lda_stoll_par(lda_func, rs,  z,  1) * b97mv_g(gamma_x,  wx_ss, par_x,  par_n, xs0, ts0, 0)
  + lda_stoll_par(f_pw    , rs,  z,  1) * b97mv_g(gamma_ss, wx_ss, par_ss, par_n, xs0, ts0, 0)
  + lda_stoll_par(lda_func, rs, -z, -1) * b97mv_g(gamma_x,  wx_ss, par_x,  par_n, xs1, ts1, 0)
  + lda_stoll_par(f_pw    , rs, -z, -1) * b97mv_g(gamma_ss, wx_ss, par_ss, par_n, xs1, ts1, 0):

b97mv_fos := (rs, z, xs0, xs1, ts0, ts1) ->
  lda_stoll_perp(f_pw, rs, z)
  * b97mv_g(gamma_os, wx_os, par_os, par_n, sqrt(xs0^2 + xs1^2)/sqrt(2), ts0, ts1):

f_b97mv :=  (lda_func, rs, z, xt, xs0, xs1, ts0, ts1) ->
  + b97mv_fpar(lda_func, rs, z, xs0, xs1, ts0, ts1)
  + b97mv_fos(rs, z, xs0, xs1, ts0, ts1):


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

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

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

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

(* the extra factor of 1/2 comes from the spin sum rule *)
rlda_a1 := (5/4) * 3*Pi * params_a_prefactor/X_FACTOR_C:

(* the functional is inherently unstable but that's how it is *)
rlda_f := (x, u, t) -> rlda_a1/(2*t - u/4):

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


text:
inputfile := "AOC4-input.txt":
input := StringTools:-Split(FileTools:-Text:-ReadFile(inputfile), "\n"):
numbers := [parse(input[1])]:
input := map(StringTools:-Trim@StringTools:-Squeeze,input[3..-1]):
numboards := ceil( nops(input) / 6 );
theboards := [seq(1..numboards)]:

# create board and mark arrays for each bingo board
for i to numboards do
    tmp := map(s->map(parse, StringTools:-Split(s)),input[1..5]);
    boards[i] := Array(tmp,datatype=integer);
    marks[i] := Array(1..5, 1..5, 0, datatype=integer[1]);
    if i <> numboards then
        input := input[7..-1];
    else
        input := [];
    end if;
end do:

# find the location of n on each board and mark it
drawnumber := proc(n)
global boards, marks, theboards;
local i, loc;
    for i in theboards do
        if member(n, boards[i], 'loc') then
            marks[i][loc] := 1;
        end if;
    end do;
    return NULL;
end proc:

# check if board i has a bingo marked
checkbingo :=proc(i)
global marks;
    if max(ArrayTools:-AddAlongDimension(marks[i], 1)) = 5 then
        return true;
    elif max(ArrayTools:-AddAlongDimension(marks[i], 2)) = 5 then
        return true;
    else
        return false;
    end if;
end proc:

# first problem - find the board that will win FIRST
answer1 := FAIL:
for j to nops(numbers) while answer1 = FAIL do
    n := numbers[j];
    drawnumber(n);
    for i to numboards do
        if checkbingo(i) then
            answer1 := add(boards[i]-boards[i]*marks[i])*n;
            break;
        end if;
    end do;
end do:
answer1;

# second problem - find the board that will win LAST
answer2 := FAIL:
numbingos := 0:
for j to nops(numbers) while answer2 = FAIL do
    n := numbers[j];
    drawnumber(n);
    for i in theboards do
        if checkbingo(i) then
            if nops(theboards) = 1 then
                answer2 := add(boards[i]-boards[i]*marks[i])*n;
            end if;
            theboards := subs(i=NULL, theboards);
        end if;
    end do;
end do:
answer2;


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

exp4_a1 := 199.81:
exp4_a2 := 4.3476:
exp4_c1 := 0.8524:
exp4_c2 := 1.2264:

# This is Eq. (40) of the paper.
exp4_f0 := s -> exp4_c1*(1 - exp(-exp4_a1*s^2)) + exp4_c2*(1 - exp(-exp4_a2*s^4)):
exp4_f  := x -> exp4_f0(X2S*x):

f := (rs, zeta, xt, xs0, xs1) -> gga_kinetic(exp4_f, rs, zeta, xs0, xs1):


text:
`is_element/autorel` := (A::set) -> proc(R)
 `is_element/rel`(A,A)(R);
end;

`is_equal/autorel` := (A::set) -> (R,S) -> `is_equal/rel`(A,A)(R,S);
`is_leq/autorel`   := (A::set) -> (R,S) -> `is_leq/rel`(A,A)(R,S);
`bot/autorel` := (A::set) -> `bot/rel`(A,A);
`top/autorel` := (A::set) -> `top/rel`(A,A);

`is_a_function/autorel` := (A::set) -> (R) -> `is_a_function/rel`(A,A)(R);
`is_total/autorel` := (A::set) -> (R) -> `is_total/rel`(A,A)(R);

`op/autorel` := (A::set) -> (R) -> `op/rel`(A,A)(R);

`hash/autorel` := (A::set) -> (R) -> `hash/rel`(A,A)(R);
`unhash/autorel` := (A::set) -> (r) -> `unhash/rel`(A,A)(r);

`id/autorel` := (A::set) -> `id/rel`(A);

`o/autorel` := (A::set) -> (S,R) -> `o/rel`(A,A,A)(S,R);

`is_reflexive/autorel` := (A::set) -> proc(R)
 return evalb(`id/autorel`(A) minus R = {});
end;

`is_irreflexive/autorel` := (A::set) -> proc(R)
 return evalb(`id/autorel`(A) intersect R = {});
end;

`is_symmetric/autorel` := (A::set) -> proc(R)
 return evalb(R = `op/autorel`(A)(R));
end;

`is_antisymmetric/autorel` := (A::set) -> proc(R)
 return evalb(R intersect `op/autorel`(A)(R) = {});
end;

`is_transitive/autorel` := (A::set) -> proc(R)
 local RR;
 RR := `o/autorel`(A)(R,R);
 return evalb(`is_leq/autorel`(A)(RR,R));
end;

`is_separated/autorel` := (A::set) -> proc(R)
 return evalb((R intersect `op/autorel`(A)(R)) minus `id/autorel`(A) = {});
end;

`list_elements/autorel` := (A::set) -> `list_elements/rel`(A,A);
`count_elements/autorel` := (A::set) -> 2^(nops(A)^2);

`random_element/autorel` := (A::set) -> (p_) -> `random_element/rel`(A,A)(args);

`transitive_closure/autorel` := (A::set) -> proc(R)
 local S,n0,n1;

 S := `id/autorel`(A) union R;
 n0 := 0;
 n1 := nops(S);
 while n1 > n0 do
  n0 := n1;
  S := `o/autorel`(A)(S,S);
  n1 := nops(S);
 od;

 return S;
end;

`cofunctor/autorel` := (A::set,B::set) -> (f) -> proc(R)
 local AA;

 AA := `top/autorel`(A);
 return select(aa -> member([f[aa[1]],f[aa[2]]],R),AA);
end;

`reindex_standard/autorel` := (A::set) -> proc(R)
 local n,i,f,N,R0;

 n := nops(A);
 f := table();
 for i from 1 to n do f[A[i]] := i; od;
 N := {seq(i,i=1..n)};
 R0 := map(aa -> [f[aa[1]],f[aa[2]]],R);
 return [N,R0];
end;


text:
# int foo(void) {
#   int n = 1;
#   int i;
#   for (i=1; i<10; i++)
#     n = n * i;
# }

func $foo (var %i i32) i32 {
  while (
    lt i32 i32 ( dread i32 %i, constval i32 10) ) {
      dassign %i (
	mul i32 (dread i32 %i, dread i32 %i))
      dassign %i (
        add i32 (dread i32 %i, constval i32 1) )}
  return ( constval i32 1 )
}
 # 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:
# https://matplotlib.org/users/customizing.html

#figure.titlesize : small        # size of the figure title (Figure.suptitle())
figure.figsize   : 8, 5         # figure size in inches
#figure.dpi       : 200          # figure dots per inch
figure.facecolor : white
figure.edgecolor : white
figure.autolayout : True        # Automatically apply 'plt.tight_layout'

lines.linewidth : 3
lines.markersize : 10
legend.fontsize : 15
legend.handlelength : 1
axes.titlesize : 17
axes.labelsize : 15
axes.linewidth : 1
axes.grid : True
grid.linestyle : :
xtick.labelsize : 14 # font size of tick labels
xtick.direction : in
ytick.direction : in
ytick.labelsize : 14
xtick.major.size : 10      # major tick size in points
xtick.minor.size : 5     # minor tick size in points
ytick.major.size : 10
ytick.minor.size : 5

#xtick.major.size : 20
#xtick.major.width : 4
#xtick.minor.size : 10
#xtick.minor.width : 2

xtick.minor.visible  : True   # turn minor ticks on by default
ytick.minor.visible  : True

# --- I like this but most people don't have it
#font.family : Optima

# --- this one is still distinct but should be more widely usable
font.family : Courier New
font.weight : bold
axes.labelweight : bold


text:
with(DifferentialAlgebra):
ring_diff := DifferentialRing(blocks = [[x2, x4, x1, x3], [y] ], derivations = [t]):
sys := [
diff(x2(t), t) - (alpha*x1(t) - beta*x2(t)),
diff(x4(t), t) - ((gama*sigma*x2(t)*x4(t) - delta*sigma*x3(t)*x4(t)) / (x3(t))),
diff(x1(t), t) - ((-b*c*x1(t) - b*x1(t)*x4(t) + 1) / (c + x4(t))),
diff(x3(t), t) - (gama*x2(t) - delta*x3(t)),
y(t) - (x1(t))
];
res := CodeTools[CPUTime](RosenfeldGroebner(sys, ring_diff, singsol=none));


text:
######################################################################
# Maps of finite sets

`is_element/maps` := (A::set,B::set) -> proc(p)
 local dom,cod;
 global reason;

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

 # Note, if p has indices 5, 6 and 7 then indices(p) will return
 # [5],[6],[7] (an unbracketed list with each individual entry in
 # square brackets).  The construct map(op,{indices(p)}) gives
 # {5,6,7} instead.
 dom := map(op,{indices(p)});
 if dom <> A then
  reason := [convert(procname,string),"wrong domain",dom,A];
  return false;
 fi;

 # Similarly, if p[5]=55 and p[6]=66 and p[7]=77 then entries(p)
 # is [55],[66],[77] but map(op,{entries(p)}) is {55,66,77}.
 cod := map(op,{entries(p)});
 if cod minus B <> {} then
  reason := [convert(procname,string),"wrong codomain",cod,B];
  return false;
 fi;

 return true;
end;

`is_equal/maps` := (A::set,B::set) -> proc(p,q)
 local a;
 global reason;

 for a in A do
  if p[a] <> q[a] then
   reason := [convert(procname,string),"different value",a,p[a],q[a]];
   return false;
  fi;
 od;

 return true;
end;

`is_leq/maps` := NULL;

`random_element/maps` := (A::set,B::set) -> proc()
 local p,n,r,a;

 p := table();
 if B = {} then
  if A = {} then
   return eval(p);
  else
   return FAIL;
  fi;
 fi;

 # nops(B) gives the number of operands in B.  Here B is supposed
 # to be a set, so the operands are just the elements.
 n := nops(B);

 # Note that the line below defines r to be a function such that each
 # invocation of r() will return a new random number in the range
 # {1,...,n}.
 r := rand(1..n);

 for a in A do p[a] := B[r()]; od;

 return eval(p);
end;

# To list the maps from A to B, we let a be the first element in A,
# and put A0 = A \ {a}.  By induction we can enumerate all the maps
# p0 : A0 -> B, and we generate all possible extensions of p0 by
# sending a to each possible element of B.

`list_elements/maps` := proc(A::set,B::set)
 local A0,a,b,M0,M,p0,p;

 # The unique map from the empty set to B is represented by a
 # table with no entries; this starts the induction.
 if A = {} then
  return [table()];
 fi;

 A0 := sort([op(A)]);
 a := A0[1];
 A0 := {op(2..-1,A0)};
 M0 := `list_elements/maps`(A0,B);
 M := [];

 for b in B do
  for p0 in M0 do
   p := copy(p0);
   p[a] := b;

   # The next line is the standard way to append an element to the
   # list M.  We use op() to strip off the square brackets, and then
   # put them back in with the new element eval(p) inside.  The
   # reason for eval(p) instead of p is technical and will not be
   # explained here, except to say that p is a table and it is
   # often good to apply eval() to tables if they are going to be
   # returned as values of functions or otherwise passed around.
   M := [op(M),eval(p)];
  od;
 od;

 return M;
end;

`count_elements/maps` := (A::set,B::set) -> nops(B) ^ nops(A);

# fibres(A,B)(p) assumes that p is a function from A to B, and
# returns a table indexed by the elements of B, whose entries
# are the fibres p^{-1}{b}.

fibres := (A::set,B::set) -> proc(p)
 local F,a,b;

 F := table();
 for b in B do F[b] := {}; od;
 for a in A do F[p[a]] := {op(F[p[a]]),a}; od;
 for b in B do F[b] := sort(F[b]); od;

 return eval(F);
end;

# `id/maps`(A) is the identity map A -> A
`id/maps` := proc(A::set) local a; table({seq(a=a,a in A)}); end:

# `compose/maps`(A,B,C)(p,q) assumes that p : A -> B and q : B -> C,
# and it returns the composite q o p : A -> C.
`compose/maps` := (A::set,B::set,C::set) -> proc(p,q) local a;
  table({seq(a = q[p[a]],a in A)}); end:


text:
with(Groebner):
with(PolynomialIdeals):

J := PolynomialIdeal({2*f1*f2*f3*f4*f5 - 9823275,
5*f1*f2*f3*f4 + 24/5*f1*f2*f3*f5 + 21/5*f1*f2*f4*f5 + 16/5*f1*f3*f4*f5 + 9/5*f2*f3*f4*f5 - 4465125,
4*f1*f2*f3 + 6*f1*f2*f4 + 18/5*f1*f2*f5 + 6*f1*f3*f4 + 24/5*f1*f3*f5 + 14/5*f1*f4*f5 + 4*f2*f3*f4 + 18/5*f2*f3*f5 + 14/5*f2*f4*f5 + 8/5*f3*f4*f5 - 441486,
3*f1*f2 + 4*f1*f3 + 4*f1*f4 + 12/5*f1*f5 + 3*f2*f3 + 4*f2*f4 + 12/5*f2*f5 + 3*f3*f4 + 12/5*f3*f5 + 7/5*f4*f5 - 15498,
2*f1 + 2*f2 + 2*f3 + 2*f4 + 6/5*f5 - 215,
f1 + 2*f2 + 3*f3 + 4*f4 + 5*f5 + 6*t}):
st := time[real]():
Basis(J, tdeg(f1,f2,f3,f4,f5,t), characteristic=2^31-1):
time[real]() - st;

J := PolynomialIdeal({2*f1*f2*f3*f4*f5*f6 - 1404728325,
6*f1*f2*f3*f4*f5 + 35/6*f1*f2*f3*f4*f6 + 16/3*f1*f2*f3*f5*f6 + 9/2*f1*f2*f4*f5*f6 + 10/3*f1*f3*f4*f5*f6 + 11/6*f2*f3*f4*f5*f6 - 648336150,
5*f1*f2*f3*f4 + 8*f1*f2*f3*f5 + 14/3*f1*f2*f3*f6 + 9*f1*f2*f4*f5 + 7*f1*f2*f4*f6 + 4*f1*f2*f5*f6 + 8*f1*f3*f4*f5 + 7*f1*f3*f4*f6 + 16/3*f1*f3*f5*f6 + 3*f1*f4*f5*f6 + 5*f2*f3*f4*f5 + 14/3*f2*f3*f4*f6 + 4*f2*f3*f5*f6 + 3*f2*f4*f5*f6 +
5/3*f3*f4*f5*f6 - 67597623,
4*f1*f2*f3 + 6*f1*f2*f4 + 6*f1*f2*f5 + 7/2*f1*f2*f6 + 6*f1*f3*f4 + 8*f1*f3*f5 + 14/3*f1*f3*f6 + 6*f1*f4*f5 + 14/3*f1*f4*f6 + 8/3*f1*f5*f6 + 4*f2*f3*f4 + 6*f2*f3*f5 + 7/2*f2*f3*f6 + 6*f2*f4*f5 + 14/3*f2*f4*f6 + 8/3*f2*f5*f6 + 4*f3*f4*f5 + 7/2*f3*f4*f6 + 8/3*f3*f5*f6 +
3/2*f4*f5*f6 - 2657700,
3*f1*f2 + 4*f1*f3 + 4*f1*f4 + 4*f1*f5 + 7/3*f1*f6 + 3*f2*f3 + 4*f2*f4 + 4*f2*f5 + 7/3*f2*f6 + 3*f3*f4 + 4*f3*f5 + 7/3*f3*f6 + 3*f4*f5 + 7/3*f4*f6 + 4/3*f5*f6 - 46243,
2*f1 + 2*f2 + 2*f3 + 2*f4 + 2*f5 + 7/6*f6 - 358}):
st := time[real]():
Basis(J, tdeg(f1,f2,f3,f4,f5,f6), characteristic=2^31-1):
time[real]() - st;

J := PolynomialIdeal({2*f1*f2*f3*f4*f5*f6*f7 - 273922023375, 7*f1*f2*f3*f4*f5*f6 + 48/7*f1*f2*f3*f4*f5*f7 + 45/7*f1*f2*f3*f4*f6*f7 + 40/7*f1*f2*f3*f5*f6*f7 + 33/7*f1*f2*f4*f5*f6*f7 + 24/7*f1*f3*f4*f5*f6*f7 + 13/7*f2*f3*f4*f5*f6*f7 - 127830277575,
6*f1*f2*f3*f4*f5 + 10*f1*f2*f3*f4*f6 + 40/7*f1*f2*f3*f4*f7 + 12*f1*f2*f3*f5*f6 + 64/7*f1*f2*f3*f5*f7 + 36/7*f1*f2*f3*f6*f7 + 12*f1*f2*f4*f5*f6 + 72/7*f1*f2*f4*f5*f7 + 54/7*f1*f2*f4*f6*f7 + 30/7*f1*f2*f5*f6*f7 + 10*f1*f3*f4*f5*f6 + 64/7*f1*f3*f4*f5*f7 + 54/7*f1*f3*f4*f6*f7 + 40/7*f1*f3*f5*f6*f7 + 22/7f1*f4*f5*f6*f7 + 6*f2*f3*f4*f5*f6 + 40/7*f2*f3*f4*f5*f7 + 36/7*f2*f3*f4*f6*f7 + 30/7*f2*f3*f5*f6*f7 + 22/7*f2*f4*f5*f6*f7 + 12/7*f3*f4*f5*f6*f7 - 13829872635,
5*f1*f2*f3*f4 + 8*f1*f2*f3*f5 + 8*f1*f2*f3*f6 + 32/7*f1*f2*f3*f7 + 9*f1*f2*f4*f5 + 12*f1*f2*f4*f6 + 48/7*f1*f2*f4*f7 + 9*f1*f2*f5*f6 + 48/7*f1*f2*f5*f7 + 27/7*f1*f2*f6*f7 + 8*f1*f3*f4*f5 + 12*f1*f3*f4*f6 +
48/7*f1*f3*f4*f7 + 12*f1*f3*f5*f6 + 64/7*f1*f3*f5*f7 + 36/7*f1*f3*f6*f7 + 8*f1*f4*f5*f6 + 48/7*f1*f4*f5*f7 + 36/7*f1*f4*f6*f7 + 20/7*f1*f5*f6*f7 + 5*f2*f3*f4*f5 + 8*f2*f3*f4*f6 + 32/7*f2*f3*f4*f7 + 9*f2*f3*f5*f6 + 48/7*f2*f3*f5*f7 + 27/7*f2*f3*f6*f7 + 8*f2*f4*f5*f6 + 48/7*f2*f4*f5*f7 + 36/7*f2*f4*f6*f7 + 20/7*f2*f5*f6*f7 + 5*f3*f4*f5*f6 + 32/7*f3*f4*f5*f7 + 27/7*f3*f4*f6*f7 + 20/7*f3*f5*f6*f7 + 11/7*f4*f5*f6*f7 - 585849123,
4*f1*f2*f3 + 6*f1*f2*f4 + 6*f1*f2*f5 + 6*f1*f2*f6 + 24/7*f1*f2*f7 + 6*f1*f3*f4 + 8*f1*f3*f5 + 8*f1*f3*f6 + 32/7*f1*f3*f7 + 6*f1*f4*f5 + 8*f1*f4*f6 + 32/7*f1*f4*f7 + 6*f1*f5*f6 + 32/7*f1*f5*f7 + 18/7*f1*f6*f7 + 4*f2*f3*f4 + 6*f2*f3*f5 + 6*f2*f3*f6 + 24/7*f2*f3*f7 + 6*f2*f4*f5 + 8*f2*f4*f6 + 32/7*f2*f4*f7 + 6*f2*f5*f6 + 32/7*f2*f5*f7 + 18/7*f2*f6*f7 + 4*f3*f4*f5 + 6*f3*f4*f6 + 24/7*f3*f4*f7 + 6*f3*f5*f6 + 32/7*f3*f5*f7 + 18/7*f3*f6*f7 + 4*f4*f5*f6 + 24/7*f4*f5*f7 + 18/7*f4*f6*f7 + 10/7*f5*f6*f7 - 11675085,
3*f1*f2 + 4*f1*f3 + 4*f1*f4 + 4*f1*f5 + 4*f1*f6 + 16/7*f1*f7 + 3*f2*f3 + 4*f2*f4 + 4*f2*f5 + 4*f2*f6 + 16/7*f2*f7 + 3*f3*f4 + 4*f3*f5 + 4*f3*f6 + 16/7*f3*f7 + 3*f4*f5 + 4*f4*f6 + 16/7*f4*f7 + 3*f5*f6 + 16/7*f5*f7 + 9/7*f6*f7 - 116053,
2*f1 + 2*f2 + 2*f3 + 2*f4 + 2*f5 + 2*f6 + 8/7*f7 - 553}):
st := time[real]():
Basis(J, tdeg(f1,f2,f3,f4,f5,f6,f7), characteristic=2^31-1):
time[real]() - st;


text:
program sample13; { square root }
var x, high, low, mid, can : integer;
begin
	writeln('Input x for calculating root x');
	readln(x);
	if x < 0 then begin
		writeln('can not calculate a root of negative number')
	end
	else begin
		low := 0;
		high := 181;
		while (high - low) >= 2 do begin
			mid := (high + low) div 2;
			can := mid * mid;
			if x < can then high := mid
			else if can < x then low := mid
			else begin
				high := mid;
				low := mid
			end
		end;
		if high = low then writeln('root ', x, ' = ', low)
		else if (high * high - x) > (x - low * low) then
			writeln('root ', x, ' = ', low)
		else writeln('root ', x, ' = ', high)
	end
end.

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing


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

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

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

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

cap_f0 := s -> 1 - params_a_alphaoAx*s*log(1 + s)/(1 + params_a_c*log(1 + s)):
cap_f  := x -> cap_f0(X2S*x):

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


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_pw86_params *params;

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

$ifdef gga_x_rpw86_params
params_a_aa := 15*0.1234:
params_a_bb := 17.33:
params_a_cc := 0.163:
$endif

pw86_f0 := s -> (1 + params_a_aa*s^2 + params_a_bb*s^4 + params_a_cc*s^6)^(1/15):
pw86_f  := x -> pw86_f0(X2S*x):

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


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

sys := [
diff(x1(t), t) = a*x1(t) + b*x1(t) - x2(t)*c*x1(t),
diff(x2(t), t) = -a*b*x2(t) + d*x2(t)*x1(t),
y1(t) = x1(t)
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


text:
program sample29;
	/* Calculator program	*/
	/*  c n : clear & set n	*/
	/*  + n : add n			*/
	/*  - n : subtract n	*/
	/*  * n : multiply n	*/
	/*  / n : divide by n		*/
	/*  o : off				*/
var unused1: integer;
    UnusedArrayForTest: array[200] of char;
procedure gcmlcm(m, n, gc, lc : integer);
	/* gc := GCM(m,n), lc := LCM(m,n) */
	/* m and n are not changed        */
var a,b,r : integer;
begin
a := m;
b := n;
while b <> 0 do begin
r := a - (a div b) * b;
a := b;
b := r
end;
gc := a;
lc := (m div gc) * n
end;
procedure abs(a, b : integer);
begin
if a < 0 then b := -a  else b := a
end;
procedure gcm(a, b, gc : integer);
var lc, aa, bb : integer;
begin
if (a = 0) or (b = 0) then gc := 1
else begin
call abs(a,aa);  call abs(b,bb);
call gcmlcm(aa, bb, gc, lc)
end
end;
procedure lcm(a, b, lc : integer);
var gc, aa, bb : integer;
begin
if (a = 0) or (b = 0) then lc := 1
else begin
call abs(a,aa);  call abs(b,bb);
call gcmlcm(aa, bb, gc, lc)
end
end;
var unusedchar : char;
procedure reduce(a1, a2 : integer);
  var gc:integer;
begin
if a1 = 0 then begin a2 := 1; return end;
if a2 = 0 then begin a1 := 1; return end;
if a2 < 0 then begin a1 := -a1;      a2 := - a2 end;
 call gcm(a1, a2, gc);
a1 := a1 div gc;
a2 := a2 div gc
end;

procedure sum(x1, x2, y1, y2 : integer);
var lc, y11 : integer;
begin
call lcm(x2, y2, lc);
x1 := x1 * (lc div x2);
y11 := y1 * (lc div y2);
x1 := x1 + y11;
x2 := lc;
call reduce(x1, x2)
end;
procedure sub(x1, x2, y1, y2 : integer);
var lc, y11 : integer;
begin
call sum(x1, x2, -y1, y2)
end;
procedure mult(x1, x2, y1, y2 : integer);
var gc,y22,y11 : integer;
begin
call gcm(x1, y2, gc);
x1 := x1 div gc;
y22 := y2 div gc;
call gcm(x2, y1, gc);
x2 := x2 div gc;
y11 := y1 div gc;
x1 := x1 * y11;
x2 := x2 * y22;
call reduce(x1, x2)
end;
procedure divide(x1, x2, y1, y2 : integer);
begin
call mult(x1, x2, y2, y1)
end;
var unusedarray : array[100] of char;
procedure printfinal(a, b : integer);
begin
if a = 0 then writeln('Final Result =', a)
else if b = 1 then writeln('Final Result =', a)
else writeln('Final Result =', a, '/', b)
end;
procedure printtemp(a, b : integer);
begin
if a = 0 then writeln('Temporary Result =', a)
else if b = 1 then writeln('Temporary Result =', a)
else writeln('Temporary Result =', a, '/', b)
end;

var	x1, x2, y1,y2 : integer;
var	com :char;
endflag : boolean;
begin
writeln('   *** Calculator -- h for help ***');
x1 := 0;  x2 :=1;
endflag := false;
while not endflag do begin
writeln(' Please input command :');
readln(com, y1);  y2 := 1;
if (com = 'c')  or (com = 'C') then begin x1 := y1; x2 := y2 end
else if com = '+' then call sum(x1, x2, y1, y2)
else if com = '-' then call sub(x1, x2, y1, y2)
else if com = '*' then call mult(x1, x2, y1, y2)
else if com = '/' then call divide(x1, x2, y1, y2)
else if (com = 'o') or (com = 'O') then endflag := true
else begin
writeln;
writeln('Calculator Usage:');
writeln('  c number : clear & set it');
writeln('  + number : add it');
writeln('  - number : subtract it');
writeln('  * number : multiply it');
writeln('  / number : divide by it');
writeln('  o        : off(terminate execution)');
writeln
end;
if endflag then call printfinal(x1, x2)
else call printtemp(x1, x2)
end
end.


text:
# Definitions for Parallel Robot Dynamics Code Generation
# Einleitung
# Erstelle Definitionen für die Maple-Skripte zur Berechnung von Kinematik und Dynamik des parallelen Roboters
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# para -> Berechnung für eine parallelen Roboter
# definitions -> Definitionen
# 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
# Init
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
read "../helper/proc_MatlabExport":
read "../helper/proc_convert_t_s":
# Lese Umgebungsvariable für Codegenerierung.
read "../robot_codegen_definitions/robot_env_par":
printf("Generiere Parameter für %s\n",robot_name):
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", leg_name):
# Link-Index, für den die Jacobi-Matrix aufgestellt wird. Hier wird angenommen, dass der Endeffektor das letzte Segment (=Link) ist. Die Jacobi-Matrix kann hier aber für beliebige Segmente aufgestellt werden. (0=Basis)
LIJAC:=NL-1:
# Ergebnisse der analytischen Jacobi-Matrix (Translatorisch)
px, py, pz := 0, 0, 0:
alphaxs_base, betays_base, gammazs_base := 0, 0, 0:
read sprintf("../codeexport/%s/tmp/jacobia_transl_%d_maple.m", leg_name, LIJAC):
b_transl := b_transl:
# Lese Umgebungsvariable für Codegenerierung.
read "../robot_codegen_definitions/robot_env_par":
# Additional Definitions
# Parameter definieren, wenn nicht vorher schon geschehen
if not assigned(J_SP) then
   J_SP := Matrix(3,3,[XX,XY,XZ,XY,YY,YZ,XZ,YZ,ZZ]):
end if:
if not assigned(J_P_P) then
   J_P_P := Matrix(3,3,[XXFP,XYFP,XZFP,XYFP,YYFP,YZFP,XZFP,YZFP,ZZFP]):
end if:
if not assigned(xE_s) then
   xE_s:=<x_all[1];x_all[2];x_all[3];x_all[4];x_all[5];x_all[6]>:
end if:

# Gravity vector in world frame
if not assigned(g_world) then
  g_world := Matrix(3, 1):
  g_world(1 .. 3, 1) := <g1, g2, g3>:
end if:

# Parallel Robotics Definitions
# Erstelle für Gelenkkoordinaten, -geschwindigkeiten und -beschleunigungen für jedes Bein
J := simplify(b_transl):
vars := indets(J):
counter := 0:
thetaList := convert(vars,list):
dList := convert(d,list):
for i to NQJ do
  var := parse(sprintf("qJ%ds", i)):
  if has(thetaList, var) then
    counter := counter + 1:
  end if:
  if has(dList, var) then
    counter := counter + 1:
  end if:
end do:
NQJ_parallel := counter:
counter := 0:
for i to ColumnDimension(J) do
  if not(Equal(J(..,i),Vector(3,[0,0,0]))) then
    counter := counter + 1:
  end if:
end do:
NQJ_parallel := counter:

qJ_i_s := Matrix(NQJ_parallel,N_LEGS):
qJD_i_s := Matrix(NQJ_parallel,N_LEGS):
qJDD_i_s := Matrix(NQJ_parallel,N_LEGS):
for i from 1 to NQJ_parallel do
  for j from 1 to N_LEGS do
     k := j-1:
     qJ_i_s(i,j):=parse(sprintf("qJ%ds", i+k*NQJ_parallel)):
     qJD_i_s(i,j):=parse(sprintf("qJD%ds", i+k*NQJ_parallel)):
     qJDD_i_s(i,j):=parse(sprintf("qJDD%ds", i+k*NQJ_parallel)):
  end do:
end do:
# Erstelle Vektor der Basisdrehungen jedes Beines
frame_A_i := Matrix(3,N_LEGS):
for i to N_LEGS do
   for j to 3 do
      if has(indets(leg_frame(4..6)),leg_frame(j+3)) then
         frame_A_i(j,i) := leg_frame(j+3)[i]:
      else
         frame_A_i(j,i) := leg_frame(j+3):
      end if:
   end do:
end do:
# Plattformschwerpunkt
r_P_sP := Vector(3,symbol=r_sP):
s_P_sP := Vector(3,symbol=s_sP):
angleConvLeg := leg_frame(7):
angleConv := xE_s(7):
# Erstelle EE-Koordinaten: xE_t, xED_t, xEDD_t -> Variablen mit Zeitabhängigkeit
#                                          xE_s, xED_s, xEDD_s -> Variablen ohne Zeitabhängigkeit
xE_t:=<xE(t);yE(t);zE(t);PsiE(t);ThetaE(t);PhiE(t)>:
xED_t:=diff~(xE_t,t):
xEDD_t:=diff~(xED_t,t):
# Zähle Freiheitsgrade des Roboters und setze nicht benötigte zu null.
xE_s := Matrix(xE_s(1..6,1)):
xED_s := copy(xE_s):
xEDD_s := copy(xE_s):
FHG_trans := 0:FHG_rot := 0:
counter := 6:
for i to 6 do
   if not(xE_s(i) = 0) then
   	 Tmp := xE_s(i):
      xED_s(i) := D||Tmp:
      xEDD_s(i) := DD||Tmp:
   else
   	 counter := counter - 1:
      xED_s(i) := 0:
      xEDD_s(i) := 0:
      xE_t(i) := 0:
      xED_t(i) := 0:
      xEDD_t(i) := 0:
   end if:
end do:
NX := counter:
counter := 0:
for i from 4 to 6 do
	if not(leg_frame(i) = 0) then
		counter := counter + 1:
	end if:
end do:
legAngles := counter:
# Export
# Maple-Export
save g_world, J_P_P, s_P_sP, NX, NQJ_parallel, angleConvLeg, angleConv, r_P_sP, frame_A_i, qJ_i_s, qJD_i_s, qJDD_i_s, xE_t, xED_t, xEDD_t, xE_s, xED_s, xEDD_s, J_SP, sprintf("../codeexport/%s/tmp/para_definitions", robot_name):
varScript := <NX;NQJ_parallel;legAngles>:
MatlabExport(varScript, sprintf("../codeexport/%s/tmp/var_parallel.m", robot_name), 1);

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


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

`is_element/stasheff_trees` := (A::set) -> proc(RTT)
 local R,TT,r,rf,TT1,T,T1,i;
 global reason;

 if not(type(RTT,list) and nops(RTT) = 2) then
  reason := [convert(procname,string),"RTT cannot be split as [R,TT]",RTT];
  return false;
 fi;

 R,TT := op(RTT);

 if not(`is_element/ord`(A)(R)) then
  reason := [convert(procname,string),"R is not an ordering",R,reason];
  return false;
 fi;

 if not(`is_element/full_trees`(A)(TT)) then
  reason := [convert(procname,string),"TT is not a full tree",TT,reason];
  return false;
 fi;

 r := table();
 for i from 1 to nops(R) do r[R[i]] := i; od;
 rf := (a) -> r[a];

 for T in TT do
  T1 := map(rf,T);
  if nops(T1) <> (max(T1) - min(T1) + 1) then
   reason := [convert(procname,string),"T is not an R-interval",T,R,T1];
   return false;
  fi;
 od;

 return true;
end;

`is_equal/stasheff_trees` := (A::set) -> proc(RTT1,RTT2)
 local R1,TT1,R2,TT2;

 R1,TT1 := op(RTT1);
 R2,TT2 := op(RTT2);
 return `is_equal/ord`(A)(R1,R2) and `is_equal/trees`(A)(TT1,TT2);
end:

`is_leq/stasheff_trees` := (A::set) -> proc(RTT1,RTT2)
 local R1,TT1,R2,TT2;

 R1,TT1 := op(RTT1);
 R2,TT2 := op(RTT2);
 return `is_equal/ord`(A)(R1,R2) and `is_leq/trees`(A)(TT1,TT2);
end:

`random_element/stasheff_trees` := (A::set) -> proc()
 local n,R,TT;

 n := nops(A);

 if n = 0 then return FAIL; fi;

 R := `random_element/ord`(A)();
 TT := `random_element/standard_stasheff_trees`(n)();
 TT := map(U -> map(u -> R[u],U),TT);

 return([R,TT]);
end;

`list_elements/stasheff_trees` := proc(A::set)
 local n,RR,R,TTT0,TT0,TT,XX;

 n := nops(A);

 RR := `list_elements/ord`(A);
 TTT0 := `list_elements/standard_stasheff_trees`(n);

 XX := NULL;
 for R in RR do
  for TT0 in TTT0 do
   TT := map(U -> map(u -> R[u],U),TT0);
   XX := XX,[R,TT];
  od;
 od;
 return([XX]);
end;

`count_elements/stasheff_trees` := proc(A::set)
 local n,m;
 n := nops(A);
 m := `count_elements/standard_stasheff_trees`(n);
 return n! * m;
end;


text:
# A (p,q)-Steiner system on a set A is a family of subsets of A (called blocks),
# each of size q, such that each subset of size p in A is contained ina unique
# block.

steiner_number_test := proc(p::posint,q::posint,n::posint)
 local i,m;
 for i from 0 to p-1 do
  m := binomial(n-i,p-i)/binomial(q-i,p-i);
  if not(type(m,integer)) then
   return false;
  fi;
 od:
 return true;
end:

`is_element/steiner_systems` := (A::set) -> (p::posint,q::posint) -> proc(S)
 global reason;

 local P,B,CC;

 if not(type(S,set(set))) then
  reason := ["is_element/steiner_systems","S is not a set of sets",S];
  return false;
 fi;

 if map(op,S) minus A <> {} then
  reason := ["is_element/steiner_systems","S is not a set of subsets of A",S,A];
  return false;
 fi;

 if map(nops,S) <> {q} then
  reason := ["is_element/steiner_systems","S is not a set of sets of size q",S,q];
  return false;
 fi;

 P := combinat[choose](A,p);

 for B in P do
  CC := select(C -> B minus C = {},S);
  if nops(CC) <> 1 then
   reason := ["is_element/steiner_systems","The set B is not contained in a unique block",B,CC];
   return false;
  fi;
 od:

 return true;
end:

`res/steiner_system` := (A::set,B::set) -> proc(S)
 local T;

 if B minus A <> {} then
  error("B is not a subset of A");
 fi;

 if nops(S[1]) + nops(B) < nops(A) then
  error("B is too small");
 fi;

 T := select(C -> C union B = A,S);
 T := map(C -> C intersect B,T);

 return T;
end:

# A (1,q)-Steiner system on A is just a partition of A into disjoint blocks of
# size q.  There is an obvious way to make such a partition when A = {1,..,m*q}

`steiner_system/congruence` := proc(q,m)
 local i,j;
 [{seq(i,i=1..q*m)},
  {seq({seq(j,j=i*q+1..(i+1)*q)},i=0..m-1)},
  [1,q,m*q]];
end:

# For any prime p, the lines in (Z/p)^2 give a (2,p)-Steiner system.

`steiner_system/affine` := proc(p)
 local V,S,i,j,m,c;

 if not(isprime(p)) then
  return FAIL;
 fi;

 V := {seq(seq([i,j],j=0..p-1),i=0..p-1)};
 S := {seq({seq([i,j],j=0..p-1)},i=0..p-1),
       seq(seq({seq([i,modp(m*i+c,p)],i=0..p-1)},c=0..p-1),m=0..p-1)};
 return [V,S,[2,p,p^2]];
end:


text:
`is_element/shellings` := (A::set) -> (R) -> proc(S)
 global reason;
 local S0,S1,d,n,i,j,k,c,N,K;

 if A = {} then
  reason := ["is_element/shellings","A is empty"];
  return false;
 fi;

 if not type(S,list) then
  reason := ["is_element/shellings","S is not a list",S];
  return false;
 fi;

 S0 := `maximal_chains/partord`(A)(R);

 if {op(S)} <> S0 then
  reason := ["is_element/shellings","S is not an enumeration of the maximal simplices",S,S0];
  return false;
 fi;

 if nops(map(nops,S0)) <> 1 then
  reason := ["is_element/shellings","The maximal simplices are not all of the same dimension",S0];
  return false;
 fi;

 d := nops(S0[1]) - 1;
 n := nops(S);
 S1 := map(c -> {op(c)},S);
 N := table();

 for i from 2 to n do
  N[i] := select(j -> nops(S1[j] intersect S1[i]) = d,{seq(k,k=1..i-1)});
  if N[i] = {} then
   reason := ["is_element/shellings","The initial shelling condition fails for i",S,i];
   return false;
  fi;
 od:

 for i from 1 to n-1 do
  for j from i+1 to n do
   c := S1[i] intersect S1[j];
   K := select(k -> c minus S1[k] = {} ,N[j]);
   if K = {} then
    reason := ["is_element/shellings","The shelling condition fails for (i,j)",S,i,j];
    return false;
   fi;
  od;
 od;

 return true;
end:

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

`is_element/edge_lex_labellings` := (A::set) -> (R) -> proc(l)
 global reason;
 local H,C,d,l0,cut,inc,cmp,R0,C0,C1,c0,c1,ab,a,b,u,i;

 if `bottom_element/partord`(A)(R) = FAIL then
  reason := ["is_element/edge_lex_labellings","There is no bottom element",A,R];
  return false;
 fi;

 if `top_element/partord`(A)(R) = FAIL then
  reason := ["is_element/edge_lex_labellings","There is no top element",A,R];
  return false;
 fi;

 H := `hasse_diagram/partord`(A)(R);

 if not(is_table_on(H)(l)) then
  reason := ["is_element/edge_lex_labellings",
             "l is not a table indexed by the edges of the Hasse diagram",l,H];
  return false;
 fi;

 C := `maximal_chains/partord`(A)(R);
 if nops(map(nops,C)) <> 1 then
  reason := ["is_element/shellings","The maximal simplices are not all of the same dimension",C];
  return false;
 fi;

 d := nops(C[1]) - 1;

 l0 := (c) -> [seq(l[[c[i],c[i+1]]],i=1..nops(c)-1)];
 cmp := proc(c,d)
  local u;

  u := select(x -> x <> 0,l0(d) -~ l0(c));
  return evalb(u <> [] and u[1] > 0);
 end:

 cut := proc(c,a,b)
  local i,c0;

  i := 1;
  while c[i] <> a do i := i+1; od;
  c0 := NULL;
  while c[i] <> b do
   c0 := c0,c[i];
   i := i+1;
  od;
  c0 := [c0,b];
  return c0;
 end:

 inc := proc(c)
  local u,v,m;
  u := l0(c);
  v := [op(u),u[-1]] -~ [u[1],op(u)];
  m := min(op(v));
  return evalb(m >= 0);
 end:

 R0 := select(ab -> ab[1] <> ab[2],R);
 for ab in R0 do
  a,b := op(ab);
  C0 := select(c -> member(a,c) and member(b,c),C);
  C0 := map(cut,C0,a,b);
  C1 := select(inc,C0);

  if nops(C1) <> 1 then
   reason :=
    ["is_element/shellings",
     "The interval [a,b] does not have a unique increasing maximal chain",
     ab,C0,C1];
   return false;
  fi;

  c1 := op(C1);
  for c0 in C0 minus C1 do
   if not cmp(c1,c0) then
    reason :=
     ["is_element/shellings",
      "In the interval [a,b], the increasing chain is not lex-smallest",
      ab,c0,c1];
    return false;
   fi;
  od;
 od:

 return true;
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: work_gga_x *)

$define gga_x_pbe_tca_params
$include "gga_x_pbe.mpl"

$define gga_x_pw91_params
$include "gga_x_pw91.mpl"

malpha :=  1:
mbeta  := 19:

fab := x -> 1/(1 + exp(-malpha*(x - mbeta))):
f   := x -> (1 - fab(x))*f_pbe(x) + fab(x)*f_pw91(x):


text:
with(Groebner):
with(PolynomialIdeals):

J := PolynomialIdeal({x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10, x1*x2 + x1*x3 + x1*x4 + x1*x5 + x1*x6 + x1*x7 + x1*x8 + x1*x9 + x1*x10 + x2*x3 + x2*x4 + x2*x5 + x2*x6 + x2*x7 + x2*x8 + x2*x9 + x2*x10 + x3*x4 + x3*x5 + x3*x6 + x3*x7 + x3*x8 + x3*x9 + x3*x10 + x4*x5 + x4*x6 + x4*x7 + x4*x8 + x4*x9 + x4*x10 + x5*x6 + x5*x7 + x5*x8 + x5*x9 + x5*x10 + x6*x7 + x6*x8 + x6*x9 + x6*x10 + x7*x8 + x7*x9 + x7*x10 + x8*x9 + x8*x10 + x9*x10, x1*x2*x3 + x1*x2*x4 + x1*x2*x5 + x1*x2*x6 + x1*x2*x7 + x1*x2*x8 + x1*x2*x9 + x1*x2*x10 + x1*x3*x4 + x1*x3*x5 + x1*x3*x6 + x1*x3*x7 + x1*x3*x8 + x1*x3*x9 + x1*x3*x10 + x1*x4*x5 + x1*x4*x6 + x1*x4*x7 +
x1*x4*x8 + x1*x4*x9 + x1*x4*x10 + x1*x5*x6 + x1*x5*x7 + x1*x5*x8 + x1*x5*x9 + x1*x5*x10 + x1*x6*x7 + x1*x6*x8 + x1*x6*x9 + x1*x6*x10 + x1*x7*x8 + x1*x7*x9 + x1*x7*x10 + x1*x8*x9 + x1*x8*x10 + x1*x9*x10 + x2*x3*x4 + x2*x3*x5 + x2*x3*x6 + x2*x3*x7 + x2*x3*x8 + x2*x3*x9 + x2*x3*x10 + x2*x4*x5 + x2*x4*x6 + x2*x4*x7 + x2*x4*x8 + x2*x4*x9 + x2*x4*x10 + x2*x5*x6 + x2*x5*x7 + x2*x5*x8 + x2*x5*x9 + x2*x5*x10 + x2*x6*x7 + x2*x6*x8 + x2*x6*x9 + x2*x6*x10 + x2*x7*x8 + x2*x7*x9 + x2*x7*x10 + x2*x8*x9 + x2*x8*x10 + x2*x9*x10 + x3*x4*x5 + x3*x4*x6 + x3*x4*x7 + x3*x4*x8 + x3*x4*x9 + x3*x4*x10 + x3*x5*x6 + x3*x5*x7 + x3*x5*x8 + x3*x5*x9 + x3*x5*x10 + x3*x6*x7 + x3*x6*x8 + x3*x6*x9 + x3*x6*x10 + x3*x7*x8 + x3*x7*x9 + x3*x7*x10 + x3*x8*x9 + x3*x8*x10 + x3*x9*x10 + x4*x5*x6 + x4*x5*x7 + x4*x5*x8 + x4*x5*x9 + x4*x5*x10 + x4*x6*x7 + x4*x6*x8 + x4*x6*x9 + x4*x6*x10 + x4*x7*x8 + x4*x7*x9 + x4*x7*x10 + x4*x8*x9 + x4*x8*x10 + x4*x9*x10 + x5*x6*x7 + x5*x6*x8 + x5*x6*x9 + x5*x6*x10 + x5*x7*x8 + x5*x7*x9 + x5*x7*x10 + x5*x8*x9 + x5*x8*x10 + x5*x9*x10 + x6*x7*x8 + x6*x7*x9 + x6*x7*x10 + x6*x8*x9 + x6*x8*x10 + x6*x9*x10 + x7*x8*x9 + x7*x8*x10 + x7*x9*x10 + x8*x9*x10, x1*x2*x3*x4 + x1*x2*x3*x5 + x1*x2*x3*x6 + x1*x2*x3*x7 + x1*x2*x3*x8 + x1*x2*x3*x9 + x1*x2*x3*x10 + x1*x2*x4*x5 + x1*x2*x4*x6 + x1*x2*x4*x7 + x1*x2*x4*x8 + x1*x2*x4*x9 + x1*x2*x4*x10 + x1*x2*x5*x6 + x1*x2*x5*x7 + x1*x2*x5*x8 + x1*x2*x5*x9 + x1*x2*x5*x10 + x1*x2*x6*x7 + x1*x2*x6*x8 + x1*x2*x6*x9 + x1*x2*x6*x10 + x1*x2*x7*x8 + x1*x2*x7*x9 + x1*x2*x7*x10 + x1*x2*x8*x9 + x1*x2*x8*x10 + x1*x2*x9*x10 +
x1*x3*x4*x5 + x1*x3*x4*x6 + x1*x3*x4*x7 + x1*x3*x4*x8 + x1*x3*x4*x9 + x1*x3*x4*x10 + x1*x3*x5*x6 + x1*x3*x5*x7 + x1*x3*x5*x8 + x1*x3*x5*x9 + x1*x3*x5*x10 + x1*x3*x6*x7 + x1*x3*x6*x8 + x1*x3*x6*x9 + x1*x3*x6*x10 + x1*x3*x7*x8 + x1*x3*x7*x9 + x1*x3*x7*x10 + x1*x3*x8*x9 + x1*x3*x8*x10 + x1*x3*x9*x10 + x1*x4*x5*x6 + x1*x4*x5*x7 + x1*x4*x5*x8 + x1*x4*x5*x9 + x1*x4*x5*x10 + x1*x4*x6*x7 + x1*x4*x6*x8 + x1*x4*x6*x9 + x1*x4*x6*x10 + x1*x4*x7*x8 + x1*x4*x7*x9 + x1*x4*x7*x10 + x1*x4*x8*x9 + x1*x4*x8*x10 + x1*x4*x9*x10 + x1*x5*x6*x7 + x1*x5*x6*x8 + x1*x5*x6*x9 + x1*x5*x6*x10 + x1*x5*x7*x8 + x1*x5*x7*x9 + x1*x5*x7*x10 + x1*x5*x8*x9 + x1*x5*x8*x10 + x1*x5*x9*x10 + x1*x6*x7*x8 + x1*x6*x7*x9 + x1*x6*x7*x10 + x1*x6*x8*x9 + x1*x6*x8*x10 + x1*x6*x9*x10 + x1*x7*x8*x9 + x1*x7*x8*x10 + x1*x7*x9*x10 + x1*x8*x9*x10 + x2*x3*x4*x5 + x2*x3*x4*x6 + x2*x3*x4*x7 + x2*x3*x4*x8 + x2*x3*x4*x9 + x2*x3*x4*x10
+ x2*x3*x5*x6 + x2*x3*x5*x7 + x2*x3*x5*x8 + x2*x3*x5*x9 + x2*x3*x5*x10 + x2*x3*x6*x7 + x2*x3*x6*x8 + x2*x3*x6*x9 + x2*x3*x6*x10 + x2*x3*x7*x8 + x2*x3*x7*x9 + x2*x3*x7*x10 + x2*x3*x8*x9 + x2*x3*x8*x10 + x2*x3*x9*x10 + x2*x4*x5*x6 + x2*x4*x5*x7 + x2*x4*x5*x8 + x2*x4*x5*x9 + x2*x4*x5*x10 + x2*x4*x6*x7 + x2*x4*x6*x8 + x2*x4*x6*x9 + x2*x4*x6*x10 + x2*x4*x7*x8 + x2*x4*x7*x9 + x2*x4*x7*x10 + x2*x4*x8*x9 + x2*x4*x8*x10 + x2*x4*x9*x10 + x2*x5*x6*x7
+ x2*x5*x6*x8 + x2*x5*x6*x9 + x2*x5*x6*x10 + x2*x5*x7*x8 + x2*x5*x7*x9 + x2*x5*x7*x10 + x2*x5*x8*x9 + x2*x5*x8*x10 + x2*x5*x9*x10 + x2*x6*x7*x8 + x2*x6*x7*x9 + x2*x6*x7*x10 + x2*x6*x8*x9 + x2*x6*x8*x10 + x2*x6*x9*x10 + x2*x7*x8*x9 + x2*x7*x8*x10 + x2*x7*x9*x10 + x2*x8*x9*x10 + x3*x4*x5*x6 + x3*x4*x5*x7 + x3*x4*x5*x8 + x3*x4*x5*x9 + x3*x4*x5*x10 + x3*x4*x6*x7 + x3*x4*x6*x8 + x3*x4*x6*x9 + x3*x4*x6*x10 + x3*x4*x7*x8 + x3*x4*x7*x9 + x3*x4*x7*x10 + x3*x4*x8*x9 + x3*x4*x8*x10 + x3*x4*x9*x10 + x3*x5*x6*x7 + x3*x5*x6*x8 + x3*x5*x6*x9 + x3*x5*x6*x10 + x3*x5*x7*x8 + x3*x5*x7*x9 + x3*x5*x7*x10 + x3*x5*x8*x9 + x3*x5*x8*x10 + x3*x5*x9*x10 + x3*x6*x7*x8 + x3*x6*x7*x9 + x3*x6*x7*x10 + x3*x6*x8*x9 + x3*x6*x8*x10 + x3*x6*x9*x10 + x3*x7*x8*x9 + x3*x7*x8*x10 + x3*x7*x9*x10 + x3*x8*x9*x10 + x4*x5*x6*x7 + x4*x5*x6*x8 + x4*x5*x6*x9 + x4*x5*x6*x10 + x4*x5*x7*x8 + x4*x5*x7*x9 + x4*x5*x7*x10 + x4*x5*x8*x9 + x4*x5*x8*x10 + x4*x5*x9*x10 + x4*x6*x7*x8 + x4*x6*x7*x9 + x4*x6*x7*x10 + x4*x6*x8*x9 + x4*x6*x8*x10 + x4*x6*x9*x10 + x4*x7*x8*x9 + x4*x7*x8*x10 + x4*x7*x9*x10 + x4*x8*x9*x10 + x5*x6*x7*x8 + x5*x6*x7*x9 + x5*x6*x7*x10 + x5*x6*x8*x9 + x5*x6*x8*x10 + x5*x6*x9*x10 + x5*x7*x8*x9 + x5*x7*x8*x10 + x5*x7*x9*x10 + x5*x8*x9*x10 + x6*x7*x8*x9 + x6*x7*x8*x10 + x6*x7*x9*x10 + x6*x8*x9*x10 + x7*x8*x9*x10, x1*x2*x3*x4*x5 + x1*x2*x3*x4*x6 + x1*x2*x3*x4*x7 + x1*x2*x3*x4*x8 + x1*x2*x3*x4*x9 + x1*x2*x3*x4*x10 + x1*x2*x3*x5*x6 + x1*x2*x3*x5*x7 + x1*x2*x3*x5*x8 + x1*x2*x3*x5*x9 + x1*x2*x3*x5*x10 + x1*x2*x3*x6*x7 + x1*x2*x3*x6*x8 + x1*x2*x3*x6*x9 + x1*x2*x3*x6*x10 + x1*x2*x3*x7*x8 + x1*x2*x3*x7*x9 + x1*x2*x3*x7*x10 + x1*x2*x3*x8*x9 + x1*x2*x3*x8*x10 + x1*x2*x3*x9*x10 + x1*x2*x4*x5*x6 + x1*x2*x4*x5*x7 + x1*x2*x4*x5*x8 + x1*x2*x4*x5*x9 + x1*x2*x4*x5*x10 + x1*x2*x4*x6*x7 + x1*x2*x4*x6*x8 + x1*x2*x4*x6*x9 + x1*x2*x4*x6*x10 + x1*x2*x4*x7*x8 + x1*x2*x4*x7*x9 + x1*x2*x4*x7*x10 + x1*x2*x4*x8*x9 + x1*x2*x4*x8*x10 + x1*x2*x4*x9*x10 + x1*x2*x5*x6*x7 + x1*x2*x5*x6*x8 + x1*x2*x5*x6*x9 + x1*x2*x5*x6*x10 + x1*x2*x5*x7*x8 + x1*x2*x5*x7*x9 + x1*x2*x5*x7*x10 + x1*x2*x5*x8*x9 + x1*x2*x5*x8*x10 + x1*x2*x5*x9*x10 + x1*x2*x6*x7*x8 + x1*x2*x6*x7*x9 + x1*x2*x6*x7*x10 + x1*x2*x6*x8*x9 + x1*x2*x6*x8*x10 + x1*x2*x6*x9*x10 + x1*x2*x7*x8*x9 + x1*x2*x7*x8*x10 + x1*x2*x7*x9*x10 + x1*x2*x8*x9*x10 + x1*x3*x4*x5*x6 + x1*x3*x4*x5*x7 + x1*x3*x4*x5*x8 + x1*x3*x4*x5*x9 + x1*x3*x4*x5*x10 + x1*x3*x4*x6*x7 +
x1*x3*x4*x6*x8 + x1*x3*x4*x6*x9 + x1*x3*x4*x6*x10 + x1*x3*x4*x7*x8 + x1*x3*x4*x7*x9 + x1*x3*x4*x7*x10 + x1*x3*x4*x8*x9 + x1*x3*x4*x8*x10 + x1*x3*x4*x9*x10 + x1*x3*x5*x6*x7 + x1*x3*x5*x6*x8 + x1*x3*x5*x6*x9 + x1*x3*x5*x6*x10 + x1*x3*x5*x7*x8 + x1*x3*x5*x7*x9 + x1*x3*x5*x7*x10 + x1*x3*x5*x8*x9 + x1*x3*x5*x8*x10 + x1*x3*x5*x9*x10 + x1*x3*x6*x7*x8 + x1*x3*x6*x7*x9 + x1*x3*x6*x7*x10 + x1*x3*x6*x8*x9 + x1*x3*x6*x8*x10 + x1*x3*x6*x9*x10 + x1*x3*x7*x8*x9 + x1*x3*x7*x8*x10 + x1*x3*x7*x9*x10 + x1*x3*x8*x9*x10 + x1*x4*x5*x6*x7 + x1*x4*x5*x6*x8 + x1*x4*x5*x6*x9 + x1*x4*x5*x6*x10 + x1*x4*x5*x7*x8 + x1*x4*x5*x7*x9 + x1*x4*x5*x7*x10 + x1*x4*x5*x8*x9 + x1*x4*x5*x8*x10 + x1*x4*x5*x9*x10 + x1*x4*x6*x7*x8 + x1*x4*x6*x7*x9 + x1*x4*x6*x7*x10 + x1*x4*x6*x8*x9 + x1*x4*x6*x8*x10 + x1*x4*x6*x9*x10 + x1*x4*x7*x8*x9 + x1*x4*x7*x8*x10 + x1*x4*x7*x9*x10 + x1*x4*x8*x9*x10 + x1*x5*x6*x7*x8 + x1*x5*x6*x7*x9 + x1*x5*x6*x7*x10 + x1*x5*x6*x8*x9 + x1*x5*x6*x8*x10 + x1*x5*x6*x9*x10 +
x1*x5*x7*x8*x9 + x1*x5*x7*x8*x10 + x1*x5*x7*x9*x10 + x1*x5*x8*x9*x10 + x1*x6*x7*x8*x9 + x1*x6*x7*x8*x10 + x1*x6*x7*x9*x10 + x1*x6*x8*x9*x10 + x1*x7*x8*x9*x10 + x2*x3*x4*x5*x6 + x2*x3*x4*x5*x7 + x2*x3*x4*x5*x8 + x2*x3*x4*x5*x9 + x2*x3*x4*x5*x10 + x2*x3*x4*x6*x7 + x2*x3*x4*x6*x8 + x2*x3*x4*x6*x9 + x2*x3*x4*x6*x10 + x2*x3*x4*x7*x8 + x2*x3*x4*x7*x9 + x2*x3*x4*x7*x10 + x2*x3*x4*x8*x9 + x2*x3*x4*x8*x10 + x2*x3*x4*x9*x10 + x2*x3*x5*x6*x7 + x2*x3*x5*x6*x8 + x2*x3*x5*x6*x9 + x2*x3*x5*x6*x10 + x2*x3*x5*x7*x8 + x2*x3*x5*x7*x9 + x2*x3*x5*x7*x10 + x2*x3*x5*x8*x9 + x2*x3*x5*x8*x10 + x2*x3*x5*x9*x10
+ x2*x3*x6*x7*x8 + x2*x3*x6*x7*x9 + x2*x3*x6*x7*x10 + x2*x3*x6*x8*x9 + x2*x3*x6*x8*x10 + x2*x3*x6*x9*x10 + x2*x3*x7*x8*x9 + x2*x3*x7*x8*x10 + x2*x3*x7*x9*x10 + x2*x3*x8*x9*x10 + x2*x4*x5*x6*x7 + x2*x4*x5*x6*x8 + x2*x4*x5*x6*x9 + x2*x4*x5*x6*x10 + x2*x4*x5*x7*x8 + x2*x4*x5*x7*x9 + x2*x4*x5*x7*x10 + x2*x4*x5*x8*x9 + x2*x4*x5*x8*x10 + x2*x4*x5*x9*x10 + x2*x4*x6*x7*x8 + x2*x4*x6*x7*x9 + x2*x4*x6*x7*x10 + x2*x4*x6*x8*x9 + x2*x4*x6*x8*x10 + x2*x4*x6*x9*x10 + x2*x4*x7*x8*x9 + x2*x4*x7*x8*x10 + x2*x4*x7*x9*x10 + x2*x4*x8*x9*x10 + x2*x5*x6*x7*x8 + x2*x5*x6*x7*x9 + x2*x5*x6*x7*x10 + x2*x5*x6*x8*x9 + x2*x5*x6*x8*x10 + x2*x5*x6*x9*x10 + x2*x5*x7*x8*x9 + x2*x5*x7*x8*x10
+ x2*x5*x7*x9*x10 + x2*x5*x8*x9*x10 + x2*x6*x7*x8*x9 + x2*x6*x7*x8*x10 + x2*x6*x7*x9*x10 + x2*x6*x8*x9*x10 + x2*x7*x8*x9*x10 + x3*x4*x5*x6*x7 + x3*x4*x5*x6*x8 + x3*x4*x5*x6*x9 + x3*x4*x5*x6*x10 + x3*x4*x5*x7*x8 + x3*x4*x5*x7*x9 + x3*x4*x5*x7*x10 + x3*x4*x5*x8*x9 + x3*x4*x5*x8*x10 + x3*x4*x5*x9*x10 + x3*x4*x6*x7*x8 + x3*x4*x6*x7*x9 + x3*x4*x6*x7*x10 + x3*x4*x6*x8*x9 +
x3*x4*x6*x8*x10 + x3*x4*x6*x9*x10 + x3*x4*x7*x8*x9 + x3*x4*x7*x8*x10 + x3*x4*x7*x9*x10 + x3*x4*x8*x9*x10 + x3*x5*x6*x7*x8 + x3*x5*x6*x7*x9 + x3*x5*x6*x7*x10 + x3*x5*x6*x8*x9 + x3*x5*x6*x8*x10 + x3*x5*x6*x9*x10 + x3*x5*x7*x8*x9 + x3*x5*x7*x8*x10 + x3*x5*x7*x9*x10 + x3*x5*x8*x9*x10 + x3*x6*x7*x8*x9 + x3*x6*x7*x8*x10 + x3*x6*x7*x9*x10 + x3*x6*x8*x9*x10 + x3*x7*x8*x9*x10
+ x4*x5*x6*x7*x8 + x4*x5*x6*x7*x9 + x4*x5*x6*x7*x10 + x4*x5*x6*x8*x9 + x4*x5*x6*x8*x10 + x4*x5*x6*x9*x10 + x4*x5*x7*x8*x9 + x4*x5*x7*x8*x10 + x4*x5*x7*x9*x10 + x4*x5*x8*x9*x10 + x4*x6*x7*x8*x9 + x4*x6*x7*x8*x10 + x4*x6*x7*x9*x10 + x4*x6*x8*x9*x10 + x4*x7*x8*x9*x10 + x5*x6*x7*x8*x9 + x5*x6*x7*x8*x10 + x5*x6*x7*x9*x10 + x5*x6*x8*x9*x10 + x5*x7*x8*x9*x10 + x6*x7*x8*x9*x10, x1*x2*x3*x4*x5*x6 + x1*x2*x3*x4*x5*x7 + x1*x2*x3*x4*x5*x8 + x1*x2*x3*x4*x5*x9 + x1*x2*x3*x4*x5*x10 + x1*x2*x3*x4*x6*x7 + x1*x2*x3*x4*x6*x8 + x1*x2*x3*x4*x6*x9 + x1*x2*x3*x4*x6*x10 + x1*x2*x3*x4*x7*x8 + x1*x2*x3*x4*x7*x9 + x1*x2*x3*x4*x7*x10 + x1*x2*x3*x4*x8*x9 + x1*x2*x3*x4*x8*x10 + x1*x2*x3*x4*x9*x10 + x1*x2*x3*x5*x6*x7 + x1*x2*x3*x5*x6*x8 + x1*x2*x3*x5*x6*x9 + x1*x2*x3*x5*x6*x10 + x1*x2*x3*x5*x7*x8 + x1*x2*x3*x5*x7*x9 + x1*x2*x3*x5*x7*x10 + x1*x2*x3*x5*x8*x9 + x1*x2*x3*x5*x8*x10 + x1*x2*x3*x5*x9*x10 + x1*x2*x3*x6*x7*x8 + x1*x2*x3*x6*x7*x9 + x1*x2*x3*x6*x7*x10 + x1*x2*x3*x6*x8*x9 + x1*x2*x3*x6*x8*x10 + x1*x2*x3*x6*x9*x10 + x1*x2*x3*x7*x8*x9 + x1*x2*x3*x7*x8*x10 + x1*x2*x3*x7*x9*x10 + x1*x2*x3*x8*x9*x10 + x1*x2*x4*x5*x6*x7 + x1*x2*x4*x5*x6*x8 + x1*x2*x4*x5*x6*x9 + x1*x2*x4*x5*x6*x10 + x1*x2*x4*x5*x7*x8 + x1*x2*x4*x5*x7*x9 + x1*x2*x4*x5*x7*x10 + x1*x2*x4*x5*x8*x9 + x1*x2*x4*x5*x8*x10 + x1*x2*x4*x5*x9*x10 + x1*x2*x4*x6*x7*x8 + x1*x2*x4*x6*x7*x9 +
x1*x2*x4*x6*x7*x10 + x1*x2*x4*x6*x8*x9 + x1*x2*x4*x6*x8*x10 + x1*x2*x4*x6*x9*x10 + x1*x2*x4*x7*x8*x9 + x1*x2*x4*x7*x8*x10 + x1*x2*x4*x7*x9*x10 + x1*x2*x4*x8*x9*x10 + x1*x2*x5*x6*x7*x8 + x1*x2*x5*x6*x7*x9 + x1*x2*x5*x6*x7*x10 + x1*x2*x5*x6*x8*x9 + x1*x2*x5*x6*x8*x10 + x1*x2*x5*x6*x9*x10 + x1*x2*x5*x7*x8*x9 + x1*x2*x5*x7*x8*x10 + x1*x2*x5*x7*x9*x10 + x1*x2*x5*x8*x9*x10
+ x1*x2*x6*x7*x8*x9 + x1*x2*x6*x7*x8*x10 + x1*x2*x6*x7*x9*x10 + x1*x2*x6*x8*x9*x10 + x1*x2*x7*x8*x9*x10 + x1*x3*x4*x5*x6*x7 + x1*x3*x4*x5*x6*x8 + x1*x3*x4*x5*x6*x9 + x1*x3*x4*x5*x6*x10 + x1*x3*x4*x5*x7*x8 + x1*x3*x4*x5*x7*x9 + x1*x3*x4*x5*x7*x10 + x1*x3*x4*x5*x8*x9 + x1*x3*x4*x5*x8*x10 + x1*x3*x4*x5*x9*x10 + x1*x3*x4*x6*x7*x8 + x1*x3*x4*x6*x7*x9 + x1*x3*x4*x6*x7*x10 + x1*x3*x4*x6*x8*x9 + x1*x3*x4*x6*x8*x10 + x1*x3*x4*x6*x9*x10 + x1*x3*x4*x7*x8*x9 + x1*x3*x4*x7*x8*x10 + x1*x3*x4*x7*x9*x10 + x1*x3*x4*x8*x9*x10 + x1*x3*x5*x6*x7*x8 + x1*x3*x5*x6*x7*x9 + x1*x3*x5*x6*x7*x10 + x1*x3*x5*x6*x8*x9 + x1*x3*x5*x6*x8*x10 + x1*x3*x5*x6*x9*x10 + x1*x3*x5*x7*x8*x9 + x1*x3*x5*x7*x8*x10 + x1*x3*x5*x7*x9*x10 + x1*x3*x5*x8*x9*x10 + x1*x3*x6*x7*x8*x9
+ x1*x3*x6*x7*x8*x10 + x1*x3*x6*x7*x9*x10 + x1*x3*x6*x8*x9*x10 + x1*x3*x7*x8*x9*x10 + x1*x4*x5*x6*x7*x8 + x1*x4*x5*x6*x7*x9 + x1*x4*x5*x6*x7*x10 + x1*x4*x5*x6*x8*x9 + x1*x4*x5*x6*x8*x10 + x1*x4*x5*x6*x9*x10 + x1*x4*x5*x7*x8*x9 + x1*x4*x5*x7*x8*x10 + x1*x4*x5*x7*x9*x10 + x1*x4*x5*x8*x9*x10 + x1*x4*x6*x7*x8*x9 + x1*x4*x6*x7*x8*x10 + x1*x4*x6*x7*x9*x10 + x1*x4*x6*x8*x9*x10 + x1*x4*x7*x8*x9*x10 + x1*x5*x6*x7*x8*x9 + x1*x5*x6*x7*x8*x10 + x1*x5*x6*x7*x9*x10 + x1*x5*x6*x8*x9*x10 + x1*x5*x7*x8*x9*x10 + x1*x6*x7*x8*x9*x10 + x2*x3*x4*x5*x6*x7 + x2*x3*x4*x5*x6*x8 + x2*x3*x4*x5*x6*x9 + x2*x3*x4*x5*x6*x10 + x2*x3*x4*x5*x7*x8 + x2*x3*x4*x5*x7*x9 + x2*x3*x4*x5*x7*x10 + x2*x3*x4*x5*x8*x9 + x2*x3*x4*x5*x8*x10 + x2*x3*x4*x5*x9*x10 + x2*x3*x4*x6*x7*x8 + x2*x3*x4*x6*x7*x9 + x2*x3*x4*x6*x7*x10 + x2*x3*x4*x6*x8*x9 + x2*x3*x4*x6*x8*x10 + x2*x3*x4*x6*x9*x10 + x2*x3*x4*x7*x8*x9 + x2*x3*x4*x7*x8*x10 + x2*x3*x4*x7*x9*x10 + x2*x3*x4*x8*x9*x10 + x2*x3*x5*x6*x7*x8 + x2*x3*x5*x6*x7*x9 + x2*x3*x5*x6*x7*x10 + x2*x3*x5*x6*x8*x9 + x2*x3*x5*x6*x8*x10 + x2*x3*x5*x6*x9*x10 + x2*x3*x5*x7*x8*x9 + x2*x3*x5*x7*x8*x10 + x2*x3*x5*x7*x9*x10 + x2*x3*x5*x8*x9*x10 + x2*x3*x6*x7*x8*x9 + x2*x3*x6*x7*x8*x10 + x2*x3*x6*x7*x9*x10 + x2*x3*x6*x8*x9*x10 + x2*x3*x7*x8*x9*x10 + x2*x4*x5*x6*x7*x8 + x2*x4*x5*x6*x7*x9 + x2*x4*x5*x6*x7*x10 + x2*x4*x5*x6*x8*x9 + x2*x4*x5*x6*x8*x10 + x2*x4*x5*x6*x9*x10 + x2*x4*x5*x7*x8*x9 + x2*x4*x5*x7*x8*x10 + x2*x4*x5*x7*x9*x10 + x2*x4*x5*x8*x9*x10 + x2*x4*x6*x7*x8*x9 + x2*x4*x6*x7*x8*x10 + x2*x4*x6*x7*x9*x10 + x2*x4*x6*x8*x9*x10 + x2*x4*x7*x8*x9*x10 + x2*x5*x6*x7*x8*x9 + x2*x5*x6*x7*x8*x10 + x2*x5*x6*x7*x9*x10 + x2*x5*x6*x8*x9*x10 + x2*x5*x7*x8*x9*x10 + x2*x6*x7*x8*x9*x10 + x3*x4*x5*x6*x7*x8 + x3*x4*x5*x6*x7*x9 + x3*x4*x5*x6*x7*x10 + x3*x4*x5*x6*x8*x9 + x3*x4*x5*x6*x8*x10
+ x3*x4*x5*x6*x9*x10 + x3*x4*x5*x7*x8*x9 + x3*x4*x5*x7*x8*x10 + x3*x4*x5*x7*x9*x10 + x3*x4*x5*x8*x9*x10 + x3*x4*x6*x7*x8*x9 + x3*x4*x6*x7*x8*x10 + x3*x4*x6*x7*x9*x10 + x3*x4*x6*x8*x9*x10 + x3*x4*x7*x8*x9*x10 + x3*x5*x6*x7*x8*x9 + x3*x5*x6*x7*x8*x10 + x3*x5*x6*x7*x9*x10 + x3*x5*x6*x8*x9*x10 + x3*x5*x7*x8*x9*x10 + x3*x6*x7*x8*x9*x10 + x4*x5*x6*x7*x8*x9 + x4*x5*x6*x7*x8*x10 + x4*x5*x6*x7*x9*x10 + x4*x5*x6*x8*x9*x10 + x4*x5*x7*x8*x9*x10 + x4*x6*x7*x8*x9*x10 + x5*x6*x7*x8*x9*x10, x1*x2*x3*x4*x5*x6*x7 + x1*x2*x3*x4*x5*x6*x8 + x1*x2*x3*x4*x5*x6*x9 + x1*x2*x3*x4*x5*x6*x10 + x1*x2*x3*x4*x5*x7*x8 + x1*x2*x3*x4*x5*x7*x9 + x1*x2*x3*x4*x5*x7*x10 + x1*x2*x3*x4*x5*x8*x9 +
x1*x2*x3*x4*x5*x8*x10 + x1*x2*x3*x4*x5*x9*x10 + x1*x2*x3*x4*x6*x7*x8 + x1*x2*x3*x4*x6*x7*x9 + x1*x2*x3*x4*x6*x7*x10 + x1*x2*x3*x4*x6*x8*x9 + x1*x2*x3*x4*x6*x8*x10 + x1*x2*x3*x4*x6*x9*x10 + x1*x2*x3*x4*x7*x8*x9 + x1*x2*x3*x4*x7*x8*x10 + x1*x2*x3*x4*x7*x9*x10 + x1*x2*x3*x4*x8*x9*x10 + x1*x2*x3*x5*x6*x7*x8 + x1*x2*x3*x5*x6*x7*x9 + x1*x2*x3*x5*x6*x7*x10 + x1*x2*x3*x5*x6*x8*x9 + x1*x2*x3*x5*x6*x8*x10 + x1*x2*x3*x5*x6*x9*x10 + x1*x2*x3*x5*x7*x8*x9 + x1*x2*x3*x5*x7*x8*x10 + x1*x2*x3*x5*x7*x9*x10 + x1*x2*x3*x5*x8*x9*x10
+ x1*x2*x3*x6*x7*x8*x9 + x1*x2*x3*x6*x7*x8*x10 + x1*x2*x3*x6*x7*x9*x10 + x1*x2*x3*x6*x8*x9*x10 + x1*x2*x3*x7*x8*x9*x10 + x1*x2*x4*x5*x6*x7*x8 + x1*x2*x4*x5*x6*x7*x9 + x1*x2*x4*x5*x6*x7*x10 + x1*x2*x4*x5*x6*x8*x9 + x1*x2*x4*x5*x6*x8*x10 + x1*x2*x4*x5*x6*x9*x10 + x1*x2*x4*x5*x7*x8*x9 + x1*x2*x4*x5*x7*x8*x10 + x1*x2*x4*x5*x7*x9*x10 + x1*x2*x4*x5*x8*x9*x10 + x1*x2*x4*x6*x7*x8*x9 + x1*x2*x4*x6*x7*x8*x10 + x1*x2*x4*x6*x7*x9*x10 + x1*x2*x4*x6*x8*x9*x10 + x1*x2*x4*x7*x8*x9*x10 + x1*x2*x5*x6*x7*x8*x9 + x1*x2*x5*x6*x7*x8*x10 + x1*x2*x5*x6*x7*x9*x10 + x1*x2*x5*x6*x8*x9*x10 + x1*x2*x5*x7*x8*x9*x10 + x1*x2*x6*x7*x8*x9*x10 + x1*x3*x4*x5*x6*x7*x8 + x1*x3*x4*x5*x6*x7*x9 + x1*x3*x4*x5*x6*x7*x10 + x1*x3*x4*x5*x6*x8*x9 + x1*x3*x4*x5*x6*x8*x10 + x1*x3*x4*x5*x6*x9*x10 + x1*x3*x4*x5*x7*x8*x9 + x1*x3*x4*x5*x7*x8*x10 + x1*x3*x4*x5*x7*x9*x10 + x1*x3*x4*x5*x8*x9*x10 + x1*x3*x4*x6*x7*x8*x9 + x1*x3*x4*x6*x7*x8*x10 + x1*x3*x4*x6*x7*x9*x10 + x1*x3*x4*x6*x8*x9*x10 + x1*x3*x4*x7*x8*x9*x10 + x1*x3*x5*x6*x7*x8*x9 + x1*x3*x5*x6*x7*x8*x10 + x1*x3*x5*x6*x7*x9*x10 + x1*x3*x5*x6*x8*x9*x10 + x1*x3*x5*x7*x8*x9*x10 + x1*x3*x6*x7*x8*x9*x10 + x1*x4*x5*x6*x7*x8*x9 + x1*x4*x5*x6*x7*x8*x10 + x1*x4*x5*x6*x7*x9*x10 + x1*x4*x5*x6*x8*x9*x10 + x1*x4*x5*x7*x8*x9*x10 + x1*x4*x6*x7*x8*x9*x10
+ x1*x5*x6*x7*x8*x9*x10 + x2*x3*x4*x5*x6*x7*x8 + x2*x3*x4*x5*x6*x7*x9 + x2*x3*x4*x5*x6*x7*x10 + x2*x3*x4*x5*x6*x8*x9 + x2*x3*x4*x5*x6*x8*x10 + x2*x3*x4*x5*x6*x9*x10 + x2*x3*x4*x5*x7*x8*x9 + x2*x3*x4*x5*x7*x8*x10 + x2*x3*x4*x5*x7*x9*x10 + x2*x3*x4*x5*x8*x9*x10 + x2*x3*x4*x6*x7*x8*x9 + x2*x3*x4*x6*x7*x8*x10 + x2*x3*x4*x6*x7*x9*x10 + x2*x3*x4*x6*x8*x9*x10 + x2*x3*x4*x7*x8*x9*x10 + x2*x3*x5*x6*x7*x8*x9 + x2*x3*x5*x6*x7*x8*x10 + x2*x3*x5*x6*x7*x9*x10 + x2*x3*x5*x6*x8*x9*x10 + x2*x3*x5*x7*x8*x9*x10 + x2*x3*x6*x7*x8*x9*x10 + x2*x4*x5*x6*x7*x8*x9 + x2*x4*x5*x6*x7*x8*x10 + x2*x4*x5*x6*x7*x9*x10 + x2*x4*x5*x6*x8*x9*x10 + x2*x4*x5*x7*x8*x9*x10 + x2*x4*x6*x7*x8*x9*x10 + x2*x5*x6*x7*x8*x9*x10 + x3*x4*x5*x6*x7*x8*x9 + x3*x4*x5*x6*x7*x8*x10 + x3*x4*x5*x6*x7*x9*x10 + x3*x4*x5*x6*x8*x9*x10 + x3*x4*x5*x7*x8*x9*x10 + x3*x4*x6*x7*x8*x9*x10 + x3*x5*x6*x7*x8*x9*x10 + x4*x5*x6*x7*x8*x9*x10, x1*x2*x3*x4*x5*x6*x7*x8 + x1*x2*x3*x4*x5*x6*x7*x9 + x1*x2*x3*x4*x5*x6*x7*x10 + x1*x2*x3*x4*x5*x6*x8*x9 + x1*x2*x3*x4*x5*x6*x8*x10 + x1*x2*x3*x4*x5*x6*x9*x10 + x1*x2*x3*x4*x5*x7*x8*x9 + x1*x2*x3*x4*x5*x7*x8*x10 + x1*x2*x3*x4*x5*x7*x9*x10 + x1*x2*x3*x4*x5*x8*x9*x10 + x1*x2*x3*x4*x6*x7*x8*x9 + x1*x2*x3*x4*x6*x7*x8*x10 + x1*x2*x3*x4*x6*x7*x9*x10 + x1*x2*x3*x4*x6*x8*x9*x10 + x1*x2*x3*x4*x7*x8*x9*x10 + x1*x2*x3*x5*x6*x7*x8*x9 + x1*x2*x3*x5*x6*x7*x8*x10
+ x1*x2*x3*x5*x6*x7*x9*x10 + x1*x2*x3*x5*x6*x8*x9*x10 + x1*x2*x3*x5*x7*x8*x9*x10 + x1*x2*x3*x6*x7*x8*x9*x10 + x1*x2*x4*x5*x6*x7*x8*x9 + x1*x2*x4*x5*x6*x7*x8*x10 + x1*x2*x4*x5*x6*x7*x9*x10 + x1*x2*x4*x5*x6*x8*x9*x10 + x1*x2*x4*x5*x7*x8*x9*x10 + x1*x2*x4*x6*x7*x8*x9*x10 + x1*x2*x5*x6*x7*x8*x9*x10
+ x1*x3*x4*x5*x6*x7*x8*x9 + x1*x3*x4*x5*x6*x7*x8*x10 + x1*x3*x4*x5*x6*x7*x9*x10 + x1*x3*x4*x5*x6*x8*x9*x10 + x1*x3*x4*x5*x7*x8*x9*x10 + x1*x3*x4*x6*x7*x8*x9*x10 + x1*x3*x5*x6*x7*x8*x9*x10 + x1*x4*x5*x6*x7*x8*x9*x10 + x2*x3*x4*x5*x6*x7*x8*x9 + x2*x3*x4*x5*x6*x7*x8*x10 + x2*x3*x4*x5*x6*x7*x9*x10 + x2*x3*x4*x5*x6*x8*x9*x10 + x2*x3*x4*x5*x7*x8*x9*x10 + x2*x3*x4*x6*x7*x8*x9*x10 + x2*x3*x5*x6*x7*x8*x9*x10 + x2*x4*x5*x6*x7*x8*x9*x10 + x3*x4*x5*x6*x7*x8*x9*x10, x1*x2*x3*x4*x5*x6*x7*x8*x9 + x1*x2*x3*x4*x5*x6*x7*x8*x10 + x1*x2*x3*x4*x5*x6*x7*x9*x10 + x1*x2*x3*x4*x5*x6*x8*x9*x10 + x1*x2*x3*x4*x5*x7*x8*x9*x10 + x1*x2*x3*x4*x6*x7*x8*x9*x10 + x1*x2*x3*x5*x6*x7*x8*x9*x10 + x1*x2*x4*x5*x6*x7*x8*x9*x10 + x1*x3*x4*x5*x6*x7*x8*x9*x10 + x2*x3*x4*x5*x6*x7*x8*x9*x10, x1*x2*x3*x4*x5*x6*x7*x8*x9*x10 + 1}):

st := time[real]():
Basis(J, tdeg(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10)):
time[real]() - st;


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

$define lda_x_params
$include "lda_x.mpl"
$include "attenuation.mpl"

a_cnst := (4/(9*Pi))^(1/3)*p_a_hyb_omega_0_/2:

lda_x_ax := -RS_FACTOR*X_FACTOR_C/2^(4/3):
lda_x_erf_spin := (rs, z) ->
  lda_x_ax*opz_pow_n(z,4/3)/rs * attenuation_erf(a_cnst*rs/opz_pow_n(z,1/3)):

f_lda_x_erf := (rs, z) -> lda_x_erf_spin(rs, z) + lda_x_erf_spin(rs, -z):

f := (rs, z) -> f_lda_x_erf(rs, z):