text:
make_ico_oct_complex := proc()
 global ico_oct_complex, ico_oct_complex1;
 local n,T,T0,T1,m,v,w,vo,vi,i,j,k,E,F,e,f,nm;

 n := 102;
 T["vertices"] := [seq(i,i=1..n)];
 T["embedding_dim"] := 3;

 m := (i,j) -> 3 + 20 * i + modp(j,20):
 v := table():

 v[1] := [0,0,1]:
 v[2] := [0,0,-1]:

 for i from 0 to 19 do
  v[m(0,i)] := evalf([cos(Pi*i/10)  ,sin(Pi*i/10)  , 2] /~ sqrt(5));
  v[m(2,i)] := evalf([cos(Pi*i/10)  ,sin(Pi*i/10)  , 0]);
  v[m(4,i)] := evalf([cos(Pi*i/10)  ,sin(Pi*i/10)  ,-2] /~ sqrt(5));
 od:

 for i from 0 to 4 do
  v[m(1,4*i)]   := evalf([cos(Pi*4*i/10)*2,sin(Pi*4*i/10)*2, 1] /~ sqrt(5));
  v[m(3,4*i+2)] := evalf([cos(Pi*(4*i+2)/10)*2,sin(Pi*(4*i+2)/10)*2, -1] /~ sqrt(5));
 od:

 nm := (x) -> x /~ sqrt(add(x[i]^2,i=1..3));

 for i from 0 to 4 do
  v[m(1,4*i+2)] := nm(v[m(1,4*i  )] +~ v[m(1,4*i+4)]);
  v[m(3,4*i  )] := nm(v[m(3,4*i-2)] +~ v[m(3,4*i+2)]);
 od:

 for i from 0 to 4 do
  v[m(1,4*i+1)] := nm(v[m(1,4*i  )] +~ v[m(1,4*i+2)]);
  v[m(1,4*i+3)] := nm(v[m(1,4*i+2)] +~ v[m(1,4*i+4)]);
  v[m(3,4*i+1)] := nm(v[m(3,4*i  )] +~ v[m(3,4*i+2)]);
  v[m(3,4*i+3)] := nm(v[m(3,4*i+2)] +~ v[m(3,4*i+4)]);
 od:

 vo := table():
 for i from 1 to n do
  vo[i] := v[i] /~ `+`(op(map(abs,v[i])));
 od:

 T["embedding"] := eval(v);
 T["sphere_embedding"] := eval(v);
 T["octahedron_embedding"] := eval(vo);
 T["normal"] := copy(v);

 E := [
  seq([ 1  , m(0,i)],i=0..19),
  seq(seq([m(j,i),m(j+1,i)],i=0..19),j=0..3),
  seq([m(4,i),2],i=0..19),
  seq(seq([m(j,i),m(j,i+1)],i=0..19),j=0..4),
  seq([m(0,i),m(1,i+1)],i=1..19,2),
  seq([m(0,i),m(1,i-1)],i=1..19,2),
  seq([m(1,i+1),m(2,i)],i=1..19,2),
  seq([m(1,i-1),m(2,i)],i=1..19,2),
  seq([m(2,i),m(3,i+1)],i=1..19,2),
  seq([m(2,i),m(3,i-1)],i=1..19,2),
  seq([m(3,i+1),m(4,i)],i=1..19,2),
  seq([m(3,i-1),m(4,i)],i=1..19,2),
  NULL
 ]:

 E := sort([op(map(e -> sort([op(e)]),E))]);

 T["edges"] := E;

 F := NULL;
 for e in T["edges"] do
  i,j := op(e);
  for k from max(i,j)+1 to n do
   if k <> i and k <> j and member([i,k],E) and member([j,k],E) then
    F := F,[i,j,k];
   fi;
  od:
 od:

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

 T["octahedron_edges"] := [
  seq([ 1  , m(0,i)],i=0..19,5),
  seq(seq([m(j,i),m(j+1,i)],i=0..19,5),j=0..3),
  seq([m(4,i),2],i=0..19,5),
  seq([m(2,i),m(2,i+1)],i=0..19),
  NULL
 ]:

 T["octahedron_edges"] := sort(map(sort,T["octahedron_edges"])):

 T["octahedron_basic_vertices"] :=
  sort([1,2,m(2,0),m(2,5),m(2,10),m(2,15)]);

 T["octahedron_basic_edges"] := [
  seq([1,m(2,5*i)],i=0..3),
  seq([m(2,5*i),m(2,5*(i+1))],i=0..3),
  seq([m(2,5*i),2],i=0..3)
 ]:

 T["octahedron_basic_edges"] := sort(map(sort,T["octahedron_basic_edges"])):

 T["octahedron_basic_faces"] := [
  seq([1,m(2,5*i),m(2,5*(i+1))],i=0..3),
  seq([2,m(2,5*i),m(2,5*(i+1))],i=0..3)
 ]:

 T["octahedron_basic_faces"] := sort(map(sort,T["octahedron_basic_faces"])):

 T["icosahedron_edges"] := [
  seq([ 1  , m(0,4*i)],i=0..4),
  seq([m(0,4*i),m(1,4*i)],i=0..4),
  seq([m(1,4*i),m(2,4*i+1)],i=0..4),
  seq([m(1,4*i),m(2,4*i-1)],i=0..4),
  seq([m(2,4*i+1),m(3,4*i+2)],i=0..4),
  seq([m(2,4*i-1),m(3,4*i-2)],i=0..4),
  seq([m(3,4*i+2),m(4,4*i+2)],i=0..4),
  seq([m(4,4*i+2),2],i=0..4),
  seq([m(1,i),m(1,i+1)],i=0..19),
  seq([m(3,i),m(3,i+1)],i=0..19),
  NULL
 ]:

 T["icosahedron_edges"] := sort(map(sort,T["icosahedron_edges"])):

 T["icosahedron_basic_vertices"] :=
  sort([1,2,seq(m(1,4*i),i=0..4),seq(m(3,4*i+2),i=0..4)]);

 T["icosahedron_basic_edges"] := [
  seq([ 1  , m(1,4*i)],i=0..4),
  seq([m(1,4*i),m(3,4*i+2)],i=0..4),
  seq([m(1,4*i),m(3,4*i-2)],i=0..4),
  seq([m(3,4*i+2),2],i=0..4),
  seq([m(1,4*i),m(1,4*(i+1))],i=0..4),
  seq([m(3,4*i+2),m(3,4*i-2)],i=0..4),
  NULL
 ]:

 T["icosahedron_basic_edges"] := sort(map(sort,T["icosahedron_basic_edges"])):

 T["icosahedron_basic_faces"] := [
  seq([1,m(1,4*i),m(1,4*(i+1))],i=0..4),
  seq([m(1,4*i),m(3,4*i+2),m(1,4*i+4)],i=0..4),
  seq([m(3,4*i-2),m(1,4*i),m(3,4*i+2)],i=0..4),
  seq([2,m(3,4*i+2),m(3,4*i-2)],i=0..4),
  NULL
 ]:

 T["icosahedron_basic_faces"] := sort(map(sort,T["icosahedron_basic_faces"])):

 w := table():

 for k from 1 to 20 do
  f := T["icosahedron_basic_faces"][k];
  w[k] := (v[f[1]] +~ v[f[2]] +~ v[f[3]]) /~ 3;
  w[k] := w[k] /~ add(w[k][i]^2,i=1..3);
 od:

 T["icosahedron_forms"] := eval(w):

 T["icosahedron_norm"] := (x) ->
  max(seq(add(x[i]*w[j][i],i=1..3),j=1..20));

 vi := table():

 for i from 1 to n do
  vi[i] := v[i] /~ T["icosahedron_norm"](v[i]);
 od:

 T["icosahedron_embedding"] := eval(vi);

 ico_oct_complex := eval(T):

 return eval(T);
end:

`subdivide/ico_oct_complex` := proc(T)
 local T0,T1,ix,n0,V0,i,x;

 T0 := `triangular_subdivision/simplicial_complex`(T):
 `normalise_embedding/simplicial_complex`(T0):
 T0["normal"] := copy(T0["embedding"]):

 T0["octahedron_edges"]  := select(e -> nops(e[1]) = 1 and member(e[2],T["octahedron_edges"] ),T0["edges"]):
 T0["icosahedron_edges"] := select(e -> nops(e[1]) = 1 and member(e[2],T["icosahedron_edges"]),T0["edges"]):
 ix := table():
 V0 := T0["vertices"]:
 n0 := nops(V0):
 for i from 1 to n0 do ix[V0[i]] := i; od:
 T1 := table():
 T1["vertices"] := [seq(i,i=1..n0)]:
 T1["edges"] := map(e -> [ix[e[1]],ix[e[2]]],T0["edges"]):
 T1["faces"] := map(f -> [ix[f[1]],ix[f[2]],ix[f[3]]],T0["faces"]):
 T1["max_simplices"] := T1["faces"]:
 T1["octahedron_edges"]  := map(e -> [ix[e[1]],ix[e[2]]],T0["octahedron_edges"]):
 T1["icosahedron_edges"] := map(e -> [ix[e[1]],ix[e[2]]],T0["icosahedron_edges"]):
 T1["octahedron_basic_vertices"]  := T["octahedron_basic_vertices"]:
 T1["octahedron_basic_edges"]     := T["octahedron_basic_edges"]:
 T1["octahedron_basic_faces"]     := T["octahedron_basic_faces"]:
 T1["icosahedron_basic_vertices"] := T["icosahedron_basic_vertices"]:
 T1["icosahedron_basic_edges"]    := T["icosahedron_basic_edges"]:
 T1["icosahedron_basic_faces"]    := T["icosahedron_basic_faces"]:

 T1["icosahedron_forms"] := copy(T["icosahedron_forms"]):
 T1["icosahedron_norm"]  := copy(T["icosahedron_norm"]):

 T1["embedding_dim"] := 3:
 T1["embedding"] := table():
 T1["normal"] := table():
 T1["sphere_embedding"]      := table():
 T1["octahedron_embedding" ] := table():
 T1["icosahedron_embedding"] := table():

 for i from 1 to n0 do
  x := T0["embedding"][V0[i]];
  x := x /~ evalf(sqrt(add(x[i]^2,i=1..3)));
  T1["embedding"][i]        := x;
  T1["normal"][i]           := x;
  T1["sphere_embedding"][i] := x;
  x := x /~ add(abs(x[i]),i=1..3);
  T1["octahedron_embedding"][i] := x;
  x := x /~ T1["icosahedron_norm"](x);
  T1["icosahedron_embedding"][i] := x;
 od:

 return eval(T1):
end:

`find_paths/ico_oct_complex` := proc(T)
 local k,E,P,E0,e,e0;

 for k in ["octahedron","icosahedron"] do
  E := {op(T[cat(k,"_edges")])};
  P := NULL:

  while E <> {} do
   e := E[1];
   E := E minus {e};
   E0 := select(e0 -> member(e[-1],e0),E);
   while E0 <> {} do
    e0 := E0[1];
    if e0[1] = e[-1] then
     e := [op(e),e0[2]];
    else
     e := [op(e),e0[1]];
    fi;
    E := E minus {e0};
    E0 := select(e0 -> member(e[-1],e0),E);
   od;
   P := P,e;
  od:
  T[cat(k,"_paths")] := [P];
 od:
end:

`make_plots/ico_oct_complex` := proc(T)
 local v;

 T["embedding"] := copy(T["octahedron_embedding"]);
 `plot/simplicial_complex`(T):
 `surface_plot/simplicial_complex`(T):
 T["octahedron_plot"] := T["plot"];
 T["octahedron_surface_plot"] := T["surface_plot"];

 T["embedding"] := copy(T["icosahedron_embedding"]);
 `plot/simplicial_complex`(T):
 `surface_plot/simplicial_complex`(T):
 T["icosahedron_plot"] := T["plot"];
 T["icosahedron_surface_plot"] := T["surface_plot"];

 T["embedding"] := copy(T["sphere_embedding"]);
 `plot/simplicial_complex`(T):
 `surface_plot/simplicial_complex`(T):
 T["sphere_plot"] := T["plot"];
 T["sphere_surface_plot"] := T["surface_plot"];

 v := eval(T["embedding"]);

 T["octahedron_arc_plot"] :=
  display(
   seq(line(v[e[1]],v[e[2]],colour=blue,thickness=4),e in T1["octahedron_edges"]),
   scaling=constrained,axes=none);

 T["icosahedron_arc_plot"] :=
  display(
   seq(line(v[e[1]],v[e[2]],colour=green,thickness=4),e in T1["icosahedron_edges"]),
   scaling=constrained,axes=none);

 T["octahedron_line_plot"] :=
  display(
   seq(line(v[e[1]],v[e[2]],colour=blue,thickness=4),e in T1["octahedron_basic_edges"]),
   scaling=constrained,axes=none);

 T["icosahedron_line_plot"] :=
  display(
   seq(line(v[e[1]],v[e[2]],colour=green,thickness=4),e in T1["icosahedron_basic_edges"]),
   scaling=constrained,axes=none);

end:

`make_javascript/ico_oct_complex` := proc(T)
 local shift,s,E,i,n;

 shift := proc(u)
  if type(u,list) then
   return map(shift,u);
  else
   return u - 1;
  fi;
 end:

 n := nops(T["vertices"]);

 s := [
  ["vertices",                    sprintf("%A",[seq(i,i=0..n-1)])],
  ["edges",                       sprintf("%A",shift(T["edges"]))],
  ["faces",                       sprintf("%A",shift(T["faces"]))],
  ["octahedron_edges",            sprintf("%A",shift(T["octahedron_edges"]))],
  ["icosahedron_edges",           sprintf("%A",shift(T["icosahedron_edges"]))],
  ["octahedron_paths",            sprintf("%A",shift(T["octahedron_paths"]))],
  ["icosahedron_paths",           sprintf("%A",shift(T["icosahedron_paths"]))],
  ["octahedron_basic_vertices",   sprintf("%A",shift(T["octahedron_basic_vertices"]))],
  ["octahedron_basic_edges",      sprintf("%A",shift(T["octahedron_basic_edges"]))],
  ["octahedron_basic_faces",      sprintf("%A",shift(T["octahedron_basic_faces"]))],
  ["icosahedron_basic_vertices",  sprintf("%A",shift(T["icosahedron_basic_vertices"]))],
  ["icosahedron_basic_edges",     sprintf("%A",shift(T["icosahedron_basic_edges"]))],
  ["icosahedron_basic_faces",     sprintf("%A",shift(T["icosahedron_basic_faces"]))]
 ]:

 E := copy(T["sphere_embedding"]):
 for i from 1 to n do E[i] := evalf[3](E[i]): od:
 s := [op(s),["sphere_embedding",sprintf("%A",[seq(E[i],i=1..n)])]]:
 E := copy(T["octahedron_embedding"]):
 for i from 1 to n do E[i] := evalf[3](E[i]): od:
 s := [op(s),["octahedron_embedding",sprintf("%A",[seq(E[i],i=1..n)])]]:
 E := copy(T["icosahedron_embedding"]):
 for i from 1 to n do E[i] := evalf[3](E[i]): od:
 s := [op(s),["icosahedron_embedding",sprintf("%A",[seq(E[i],i=1..n)])]]:

 s := map(u -> sprintf("\"%s\" : %s,\n",u[1],u[2]),s);
 s := cat("{\n",op(s),sprintf("\"num_vertices\" : %d\n}\n",n));

 T["javascript"] := s;
 return s;
end:

ico_oct_complex := make_ico_oct_complex():


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

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

$define gga_x_pbe_params
$include "gga_x_pbe.mpl"
$include "gvt4.mpl"

m06_alpha  := 0.00186726:
m06_coeff_d := params_a_d:

(* there is a factor if 2 in the definition of z, as in Theor. Chem. Account 120, 215 (2008) *)
(* A MINUS was missing in Eq. (7) of the paper *)

m06_f := (x, u, t) ->
  + pbe_f(x)*mgga_series_w(params_a_a, 12, t)
  + gtv4(m06_alpha, m06_coeff_d, x, 2*(t - K_FACTOR_C)):

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

ext: mpl


text:
# Teach Maple (through depends and eval) about our new binding forms.
# lam binds from 1st arg to 3rd arg.
# Branch binds from 1st arg (a pattern) to 2nd arg.
# Bind and ary bind from 2nd arg to 3rd arg.

# note that v _can_ in principle occur in t.
`depends/lam` := proc(v::name, t, e, x, $)
  depends(t, x) or depends(e, x minus {v})
end proc:

`depends/Branch` := proc(p, e, x, $)
  depends(e, x minus {Hakaru:-pattern_binds(p)})
end proc:

# note that v _can_ occur in m1.
`depends/Bind` := proc(m1, v::name, m2, x, $)
  depends(m1, x) or depends(m2, x minus {v})
end proc:

# note that i _can_ occur in n.
`depends/ary` := proc(n, i::name, e, x, $)
  depends(n, x) or depends(e, x minus {i})
end proc:
`depends/Plate` := eval(`depends/ary`):

`eval/lam` := proc(e, eqs, $)
  local v, t, ee;
  v, t, ee := op(e);
  v, ee := BindingTools:-generic_evalat(v, ee, eqs);
  t := eval(t, eqs);
  eval(op(0,e), eqs)(v, t, ee)
end proc:

`eval/Branch` := proc(e, eqs, $)
  local p, ee, vBefore, vAfter;
  p, ee := op(e);
  vBefore := [Hakaru:-pattern_binds(p)];
  vAfter, ee := BindingTools:-generic_evalat(vBefore, ee, eqs);
  eval(op(0,e), eqs)(subs(op(zip(`=`, vBefore, vAfter)), p), ee)
end proc:

`eval/Bind` := proc(e, eqs, $)
  local m1, v, m2;
  m1, v, m2 := op(e);
  eval(op(0,e), eqs)(eval(m1, eqs), BindingTools:-generic_evalat(v, m2, eqs))
end proc:

`eval/ary` := proc(e, eqs, $)
  local n, i, ee;
  n, i, ee := op(e);
  eval(op(0,e), eqs)(eval(n, eqs), BindingTools:-generic_evalat(i, ee, eqs))
end proc:
`eval/Plate` := eval(`eval/ary`):

#############################################################################
Hakaru := module ()
  option package;
  local p_true, p_false, make_piece, Mk_Plus, lift1_piecewise,
        ModuleLoad, ModuleUnload;
  export
     # These first few are smart constructors (for themselves):
         case, app, ary, idx, fst, snd, size, Datum,
     # while these are "proper functions"
         verify_measure, verify_hboolean, pattern_equiv,
         piecewise_And, map_piecewiselike, lift_piecewise, foldr_piecewise,
         flatten_piecewise,
         pattern_match, pattern_binds,
         closed_bounds, open_bounds,
         htype_patterns,
         bool_And, bool_Or, bool_Not,
         UpdateArchive;
  # These names are not assigned (and should not be).  But they are
  # used as global names, so document that here.
  global
     # Basic syntax for composing measures
         Bind, Weight, Ret, Msum, Plate, Context, Pair, _Unit, PARTITION,
     # Primitive (known) measures
         Lebesgue, Uniform, Gaussian, Cauchy, StudentT, BetaD,
         GammaD, ChiSquared,
         Counting, Categorical, NegativeBinomial, PoissonD,
     # Functions, annotated with argument type, applied using "app"
         lam,
     # Term constructors for Datum (algebraic data type)
         Inr, Inl, Et, Done, Konst, Ident,
     # The parts of "case" besides the scrutinee
         Branches, Branch,
     # Pattern constructors
         PWild, PVar, PDatum, PInr, PInl, PEt, PDone, PKonst, PIdent,
     # Verification of alpha-equivalence among measures
         measure,
     # Punctuation for matching free variables with measures for disintegration
         `&M`,
     # Structure types for Hakaru types and Hakaru "case" expressions
         known_continuous, known_discrete, t_Hakaru, t_type, t_case,
     # Structure types for piecewise-like expressions:
     # piecewise, case, and idx into literal array
         t_pw, t_piecewiselike,
     #other types
         `&implies`,
     # Type constructors for Hakaru
         AlmostEveryReal, HReal, HInt, HData, HMeasure, HArray, HFunction,
         Bound, DatumStruct;
  p_true  := 'PDatum(true,PInl(PDone))';
  p_false := 'PDatum(false,PInr(PInl(PDone)))';

  case := proc(e, bs :: specfunc(Branch(anything, anything), Branches), $)
    local ret, b, substs, eSubst, pSubst, p, binds, uncertain;
    if e :: 't_piecewiselike' then
      map_piecewiselike(procname, _passed)
    else
      ret := Branches();
      for b in bs do
        substs := pattern_match(e, e, op(1,b));
        if substs <> NULL then
          eSubst, pSubst := substs;
          p := subs(pSubst, op(1,b));
          binds := {pattern_binds(p)};
          uncertain := remove((eq -> lhs(eq) in binds), eSubst);
          if nops(uncertain) = 0 then p := PWild end if;
          ret := Branches(op(ret),
                          Branch(p, eval(eval(op(2,b), pSubst), eSubst)));
          if nops(uncertain) = 0 then break end if;
        end if
      end do;
      if ret :: Branches(Branch(identical(PWild), anything)) then
        op([1,2], ret)
      elif ret :: Branches(Branch(identical(p_true), anything),
                           Branch({identical(p_false),
                                   identical(PWild),
                                   PVar(anything)}, anything)) then
        piecewise(make_piece(e), op([1,2], ret), op([2,2], ret))
      elif ret :: Branches(Branch(identical(p_false), anything),
                           Branch({identical(p_true),
                                   identical(PWild),
                                   PVar(anything)}, anything)) then
        piecewise(make_piece(e), op([2,2], ret), op([1,2], ret))
      else
        'case'(e, ret)
      end if
    end if
  end proc;

  pattern_match := proc(e0, e, p, $)
    local x, substs, eSubst, pSubst;
    if p = PWild then return {}, {}
    elif p :: PVar(anything) then
      x := op(1,p);
      pSubst := {`if`(depends(e0,x), x=gensym(x), NULL)};
      return {subs(pSubst,x)=e}, pSubst;
    elif p = p_true then
      if e = true then return {}, {}
      elif e = false then return NULL
      end if
    elif p = p_false then
      if e = false then return {}, {}
      elif e = true then return NULL
      end if
    elif p :: PDatum(anything, anything) then
      if e :: Datum(anything, anything) then
        if op(1,e) = op(1,p) then return pattern_match(e0, op(2,e), op(2,p))
        else return NULL
        end if
      end if
    elif p :: PInl(anything) then
      if e :: Inl(anything) then return pattern_match(e0, op(1,e), op(1,p))
      elif e :: Inr(anything) then return NULL
      end if
    elif p :: PInr(anything) then
      if e :: Inr(anything) then return pattern_match(e0, op(1,e), op(1,p))
      elif e :: Inl(anything) then return NULL
      end if
    elif p :: PEt(anything, anything) then
      if e :: Et(anything, anything) then
        substs := pattern_match(e0, op(1,e), op(1,p));
        if substs = NULL then return NULL end if;
        eSubst, pSubst := substs;
        substs := pattern_match(e0, eval(op(2,e),eSubst), op(2,p));
        if substs = NULL then return NULL end if;
        return eSubst union substs[1], pSubst union substs[2];
      elif e = Done then return NULL
      end if
    elif p = PDone then
      if e = Done then return {}, {}
      elif e :: Et(anything, anything) then return NULL
      end if
    elif p :: PKonst(anything) then
      if e :: Konst(anything) then return pattern_match(e0, op(1,e), op(1,p))
      end if
    elif p :: PIdent(anything) then
      if e :: Ident(anything) then return pattern_match(e0, op(1,e), op(1,p))
      end if
    else
      error "pattern_match: %1 is not a pattern", p
    end if;
    pSubst := map((x -> `if`(depends(e0,x), x=gensym(x), NULL)),
                  {pattern_binds(p)});
    eSubst := {e=subsindets(
                   subsindets[nocache](
                     subs(pSubst,
                          p_true=true,
                          p_false=false,
                          PDatum=Datum, PInr=Inr, PInl=Inl, PEt=Et, PDone=Done,
                          PKonst=Konst, PIdent=Ident,
                          p),
                     identical(PWild),
                     p -> gensym(_)),
                   PVar(anything),
                   p -> op(1,p))};
    eSubst, pSubst
  end proc;

  make_piece := proc(rel, $)
    # Try to prevent PiecewiseTools:-Is from complaining
    # "Wrong kind of parameters in piecewise"
    if rel :: {specfunc(anything, {And,Or,Not}), `and`, `or`, `not`} then
      map(make_piece, rel)
    elif rel :: {'`::`', 'boolean', '`in`'} then
      rel
    else
      rel = true
    end if
  end proc;

  pattern_binds := proc(p, $)
    if p = PWild or p = PDone then
      NULL
    elif p :: PVar(anything) then
      op(1,p)
    elif p :: PDatum(anything, anything) then
      pattern_binds(op(2,p))
    elif p :: {PInl(anything), PInr(anything),
               PKonst(anything), PIdent(anything)} then
      pattern_binds(op(1,p))
    elif p :: PEt(anything, anything) then
      pattern_binds(op(1,p)), pattern_binds(op(2,p))
    else
      error "pattern_binds: %1 is not a pattern", p
    end if
  end proc:

  verify_hboolean := proc(a, b, $)
    local x,y,conv_tbl,and_fn,or_fn,not_fn;
    and_fn := ((x->AND({op(x)}))@`bool_And`);
    or_fn  := ((x->OR({op(x)}))@`bool_Or`);
    not_fn := bool_Not;
    conv_tbl := table([`and`=and_fn,`And`=and_fn,`Or`=or_fn,`or`=or_fn,`Not`=not_fn,`not`=not_fn]);
    x,y := subsindets([a,b], Or(`and`,`or`,`not`, specfunc({`And`,`Or`,`Not`})), x->conv_tbl[op(0,x)](op(x)))[];
    evalb(x=y);
  end proc;

  verify_measure := proc(m, n, v:='hboolean', $)
    local mv, x, i, j, k;
    mv := measure(v);
    if m :: specfunc({Bind, Plate}) and n :: specfunc({Bind, Plate}) and
       (verify(m, n, 'Bind'(mv, true, true))
      or verify(m, n, 'Plate'(v, true, true))) then
      x := gensym(cat(op(2,m), "_", op(2,n), "_"));
      thisproc(subs(op(2,m)=x, op(3,m)),
               subs(op(2,n)=x, op(3,n)), v)
    elif m :: 'specfunc(Msum)' and n :: 'specfunc(Msum)'
         and nops(m) = nops(n) then
      k := nops(m);
      verify(k, GraphTheory[BipartiteMatching](GraphTheory[Graph]({
                seq(seq(`if`(thisproc(op(i,m), op(j,n), v), {i,-j}, NULL),
                        j=1..k), i=1..k)}))[1]);
    elif andmap(type, [m,n], 'specfunc(piecewise)') and nops(m) = nops(n) then
      k := nops(m);
      verify(m, n, 'piecewise'(seq(`if`(i::even or i=k, mv, hboolean), i=1..k)));

    elif andmap(type, [m,n], 'specfunc(PARTITION)') then
      Partition:-SamePartition(verify_hboolean, verify_measure, m, n);

    elif m :: specfunc('case') and
        verify(m, n, 'case'(v, specfunc(Branch(true, true), Branches))) then
      # This code unfortunately only handles alpha-equivalence for 'case' along
      # the control path -- not if 'case' occurs in the argument to 'Ret', say.
      k := nops(op(2,m));
      for i from 1 to k do
        j := pattern_equiv(op([2,i,1],m), op([2,i,1],n));
        if j = false then return j end if;
        j := map(proc(eq, $)
                   local x;
                   x := gensym(cat(lhs(eq), "_", rhs(eq), "_"));
                   [lhs(eq)=x, rhs(eq)=x]
                 end proc, j);
        j := thisproc(subs(map2(op,1,j), op([2,i,2],m)),
                      subs(map2(op,2,j), op([2,i,2],n)), v);
        if j = false then return j end if;
      end do;
      true
    elif m :: 'LO(name, anything)' and n :: 'LO(name, anything)' then
      x := gensym(cat(op(1,m), "_", op(1,n), "_"));
      verify(subs(op(1,m)=x, op(2,m)),
             subs(op(1,n)=x, op(2,n)), v)
    elif m :: specfunc('lam') and n :: specfunc('lam') and
        verify(m, n, 'lam'(true, v, true)) then
      # m and n are not even measures, but we verify them anyway...
      x := gensym(cat(op(1,m), "_", op(1,n), "_"));
      thisproc(subs(op(1,m)=x, op(2,m)),
               subs(op(1,n)=x, op(2,n)), v)
    else
      verify(m, n, {v, Ret(v), Weight(v, mv), Context(v, mv)})
    end if
  end proc;

  pattern_equiv := proc(p, q, $) :: {identical(false),set(`=`)};
    local r, s;
    if ormap((t->andmap(`=`, [p,q], t)), [PWild, PDone]) then
      {}
    elif andmap(type, [p,q], PVar(anything)) then
      {op(1,p)=op(1,q)}
    elif andmap(type, [p,q], PDatum(anything,anything)) and op(1,p)=op(1,q) then
      pattern_equiv(op(2,p),op(2,q))
    elif ormap((t->andmap(type, [p,q], t(anything))),
               [PInl, PInr, PKonst, PIdent]) then
      pattern_equiv(op(1,p),op(1,q))
    elif andmap(type, [p,q], PEt(anything, anything)) then
      r := pattern_equiv(op(1,p),op(1,q));
      s := pattern_equiv(op(2,p),op(2,q));
      if map(lhs,r) intersect map(lhs,s) = {} and
         map(rhs,r) intersect map(rhs,s) = {} then
        r union s
      else
        false
      end if
    else
      false
    end if
  end proc;

  piecewise_And := proc(cond::{list,set,specfunc(`And`),`and`}, th, el, $)
    if nops(cond) = 0 or th = el then
      th
    else
      piecewise(And(op(cond)), th, el)
    end if
  end proc;

  map_piecewiselike := proc(f, p::t_piecewiselike)
    local i, g, h;
    if p :: 'specfunc(piecewise)' then
      piecewise(seq(`if`(i::even or i=nops(p), f(op(i,p),_rest), op(i,p)),
                    i=1..nops(p)))
    elif p :: 't_case' then
      # Mind the hygiene
      subsindets(eval(subsop(2 = map[3](applyop, g, 2, op(2,p)), p),
                      g=h(f,_rest)),
                 'typefunc(specfunc(h))',
                 (e -> op([0,1],e)(op(1,e), op(2..-1,op(0,e)))))
    elif p :: 'idx(list, anything)' then
      idx(map(f,op(1,p),_rest), op(2,p))
    else
      error "map_piecewiselike: %1 is not t_piecewiselike", p
    end if
  end proc;

  lift_piecewise := proc(e, extra:={}, $)
    local e1, e2;
    e2 := e;
    while e1 <> e2 do
      e1 := e2;
      e2 := subsindets(e1,
              '{extra,
                And(`+`, Not(specop(Not(specfunc(piecewise)), `+`))),
                And(`*`, Not(specop(Not(specfunc(piecewise)), `*`))),
                And(`^`, Not(specop(Not(specfunc(piecewise)), `^`))),
                exp(specfunc(piecewise))}',
              lift1_piecewise)
    end do
  end proc;

  lift1_piecewise := proc(e, $)
    local i, p;
    if membertype(t_piecewiselike, [op(e)], i) then
      p := op(i,e);
      if nops(p) :: even and not (e :: `*`) and op(-1,p) <> 0 then
        p := piecewise(op(p), 0);
      end if;
      map_piecewiselike((arm->lift1_piecewise(subsop(i=arm,e))), p)
    else
      e
    end if
  end proc;

  foldr_piecewise := proc(cons, nil, pw, $) # pw may or may not be piecewise
    # View pw as a piecewise and foldr over its arms
    if pw :: 'specfunc(piecewise)' then
      foldr(proc(i,x) cons(op(i,pw), op(i+1,pw), x) end proc,
            `if`(nops(pw)::odd, cons(true, op(-1,pw), nil), nil),
            seq(1..nops(pw)-1, 2))
    else
      cons(true, pw, nil)
    end if
  end proc;

  #Take a piecewise: If some of the branches also contain piecewises, attempt
  #to re-express the whole as a single piecewise. In short: fewer piecewises, more
  #branches, more complex conditions.--Carl 2016Sep09
  flatten_piecewise:= proc(PW::specfunc(piecewise), $)::specfunc(piecewise);
  local
     pwL:= [op(PW)],
     nP:= nops(pwL),
     conds:= pwL[[seq(1..nP-1, 2)]],
     branches:= pwL[[seq(2..nP-1, 2), nP]],
     innerPW:= indets(branches, specfunc(piecewise)),
     outerPW, C_O, C_I, B_I, r
  ;
     if indets(conds, specfunc(piecewise)) <> {} then
          userinfo(2, procname, "Case of inner pw in a condition not yet handled.");
          return PW
     end if;
     if innerPW = {} then
          userinfo(3, procname, "No inner pw.");
          return PW
     end if;
     if nops~(innerPW) <> {3} then
          userinfo(2, procname, "Case of inner pw of other than 3 ops not yet handled.");
          return PW
     end if;
     if nops(PW) <> 3 then
          userinfo(2, procname, "Case of outer pw of other than 3 ops not yet handled.");
          return PW
     end if;
     outerPW:= applyop(lift_piecewise, {2,3}, PW);
     if membertype(
          And(
               satisfies(e-> indets(e, specfunc(piecewise)) <> {}),
               Not(specfunc(piecewise))
          ),
          {op(2..3, outerPW)}
     )
     then
          userinfo(1, procname, "Unliftable pw.");
          return PW
     end if;
     C_O:= op(1,outerPW);
     C_I:= map(e-> `if`(e::specfunc(piecewise), op(1,e), true), [op(2..3, outerPW)]);
     B_I:= map(e-> `if`(e::specfunc(piecewise), [op(2..3, e)], [e$2])[], [op(2..3, outerPW)]);
     #For debugging purposes, I want to show the proposed output before it's passed to
     #piecewise.
     r:= [
          And(C_O, C_I[1]), B_I[1],
          And(C_O, KB:-negate_rel(C_I[1])), B_I[2],
          And(KB:-negate_rel(C_O), C_I[2]), B_I[3],
          And(KB:-negate_rel(C_O), KB:-negate_rel(C_I[2])), B_I[4]
     ];
     userinfo(3, procname, "Proposed ouput: ", print(%piecewise(r[])));
     piecewise(r[])
  end proc;

  app := proc (func, argu, $)
    if func :: 'lam(name, anything, anything)' then
      eval(op(3,func), op(1,func)=argu)
    elif func :: 't_piecewiselike' then
      map_piecewiselike(procname, _passed)
    else
      'procname(_passed)'
    end if
  end proc;

  ary := proc (n, i, e, $)
    local j;
    if e :: 'idx'('freeof'(i), 'identical'(i)) then
      # Array eta-reduction. Assume the size matches.  (We should keep array
      # size information in the KB and use it here, but we don't currently.)
      op(1,e)
    elif n :: nonnegint then
      [seq(eval(e,i=j), j=0..n-1)] # Unroll array with literal size
    else
      'procname(_passed)'
    end if
  end proc;

  idx := proc (a, i, $)
    if a :: 'ary(anything, name, anything)' then
      eval(op(3,a), op(2,a)=i)
    elif a :: 'list' and i::nonnegint then
      a[i+1]
    elif a :: 'list' and nops(convert(a,'set')) = 1 then
      a[1] # Indexing into a literal array whose elements are all the same
    elif a :: 't_piecewiselike' then
      map_piecewiselike(procname, _passed)
    else
      'procname(_passed)'
    end if
  end proc;

  #Extract the first member of a Pair.
  fst:= proc(p, $)
    if p :: 'Pair'('anything'$2) then
      op(1,p)
    elif p :: t_piecewiselike then
      map_piecewiselike(fst, p)
    else
      'procname'(p)
    end if;
  end proc;

  #Extract the second member of a Pair.
  snd:= proc(p, $)
    if p :: 'Pair'('anything'$2) then
      op(2,p)
    elif p :: t_piecewiselike then
      map_piecewiselike(snd, p)
    else
      'procname'(p)
    end if;
  end proc;

  size := proc(a, $)
    local res;
    if a :: 'ary(anything, name, anything)' then
      op(1,a)
    elif a :: 'list' then
      nops(a)
    elif a :: 't_piecewiselike' then
      map_piecewiselike(procname, _passed)
    else
      'procname(_passed)'
    end if
  end proc;

  Datum := proc(hint, payload, $)
    # Further cheating to equate Maple booleans and Hakaru booleans
    if hint = true and payload = Inl(Done) or
       hint = false and payload = Inr(Inl(Done)) then
      hint
    else
      'procname(_passed)'
    end if
  end proc;

  closed_bounds := proc(r::range, $)
    Bound(`>=`, lhs(r)), Bound(`<=`, rhs(r))
  end proc;

  open_bounds := proc(r::range, $)
    Bound(`>`, lhs(r)), Bound(`<`, rhs(r))
  end proc;

  # Creates an N-ary operator `Plus` such:
  #  Plus(Plus(a,b),c)=Plus(a,Plus(a,b))
  #  Plus(a,Zero)=a
  #  Plus(a)=a
  #  Plus(x,x)=x
  Mk_Plus := proc(plus,zero,$) proc()
    local as := {args};
    as := map(a->if op(0,a)=plus then op(a) else a end if, as);
    if nops(as)=0 then zero
    elif nops(as)=1 then op(1,as)
    else plus(op(as))
    end if;
  end proc; end proc;

  # Replacements for `and` and `or` which do not
  # evaluate "x = y" to "false" when "x,y" are unbound vars.
  bool_And := Mk_Plus('And','true' );
  bool_Or  := Mk_Plus('Or' ,'false');

  bool_Not := proc(a,$)
    if a :: t_kb_atom then
      subsindets(KB:-negate_rel(a), `not`, Not@op);
    else
      Not(a)
    end if;
  end proc;

  # Enumerate patterns for a given Hakaru type
  htype_patterns := proc(t::t_type, $)
    :: specfunc(Branch(anything, list(t_type)), Branches);
    local struct;
    uses StringTools;
    if t :: specfunc(DatumStruct(anything, list(Konst(anything))), HData) then
      foldr(proc(struct,ps,$) Branches(
              op(map((p -> Branch(PDatum(op(1,struct), PInl(op(1,p))),
                                  op(2,p))),
                     foldr(proc(kt,qs,$)
                             local p, q;
                             Branches(seq(seq(Branch(PEt(PKonst(op(1,p)),
                                                         op(1,q)),
                                                     [op(op(2,p)),op(op(2,q))]),
                                              q in qs),
                                          p in htype_patterns(op(1,kt))))
                           end proc,
                           Branches(Branch(PDone, [])), op(op(2,struct))))),
              op(map[3](applyop, PInr, [1,2], ps)))
            end proc,
            Branches(), op(t))
    else
      Branches(Branch(PVar(gensym(convert(LowerCase(op(-1, ["x", op(
                                            `if`(t::function,
                                              select(IsUpper, Explode(op(0,t))),
                                              []))])),
                                          name))),
                      [t]))
    end if
  end proc;

  ModuleLoad := proc($)
    local g; #Iterator over thismodule's globals
    ModuleUnload();
    VerifyTools[AddVerification](measure = eval(verify_measure));
    VerifyTools[AddVerification](hboolean = eval(verify_hboolean));

    TypeTools:-AddType(
         `&implies`,
         proc(e, t1, t2, $) type(e, Or(Not(t1), t2)) end proc
    );
    TypeTools:-AddType(Name, And(name, Not({constant,undefined})));
    TypeTools:-AddType(known_continuous,
         '{
              Lebesgue(anything, anything), Uniform(anything, anything),
              Gaussian(anything, anything), Cauchy(anything, anything),
              StudentT(anything, anything, anything), ChiSquared(anything),
              BetaD(anything, anything), GammaD(anything, anything),
              InverseGammaD(anything, anything)
         }'
    );
    TypeTools:-AddType(known_discrete,
         '{
              Counting(anything, anything),
              Categorical(anything),
              Binomial(anything,anything),
              NegativeBinomial(anything,anything),
              PoissonD(anything)
          }'
    );
    #This type t_Hakaru is the basic syntax checker for that part of the
    #Maple/Hakaru language accessible by the external user.--Carl 2016Sep20
    TypeTools:-AddType(t_Hakaru,
         #Alternation/conjunction of types via {} or And follows the McCarthy
         #short-cut rule: Proceeding left to right, once satisfaction of the type
         #can be determined, the remaining elements aren't evaluated. Therefore,
         #recursive types are possible by placing the base cases at the beginning.
         #But use Or instead of {} because you can't control the order of
         #expressions with {}.
         Or(
           'known_continuous', 'known_discrete',
           t_pw, #Needs to be more specific!
           t_partition,
           t_case,
          'Ret(anything)',
          'Bind(t_Hakaru, name, t_Hakaru)',
          'specfunc(t_Hakaru, Msum)',
          'Weight(algebraic, t_Hakaru)',
          'Plate(algebraic, name, t_Hakaru)',
          t_lam
         )
    );
    TypeTools:-AddType(t_type,
      '{specfunc(Bound(identical(`<`,`<=`,`>`,`>=`), anything),
                 {AlmostEveryReal, HReal, HInt}),
        specfunc(DatumStruct(anything, list({Konst(t_type), Ident(t_type)})),
                 HData),
        HMeasure(t_type),
        HArray(t_type),
        HFunction(t_type, t_type)}');
    TypeTools:-AddType(t_case,
      'case(anything, specfunc(Branch(anything, anything), Branches))');
    TypeTools:-AddType(t_pw, 'specfunc(piecewise)');
    TypeTools:-AddType(t_partition, 'specfunc(PARTITION)'); #Appropriate to put this here?
    TypeTools:-AddType(t_piecewiselike,
      '{specfunc(piecewise), t_case, idx(list, anything)}');
    TypeTools:-AddType(t_lam, 'lam(name, t_type, t_Hakaru)');

    # A temporary type which should be removed when piecewise is gone
    TypeTools:-AddType(t_pw_or_part, Or(t_pw,t_partition));

    #Protect the keywords of the Hakaru language.
    #op([2,6], ...) of a module is its globals.
    for g in op([2,6], thismodule) do
         if g <> eval(g) then
           unassign(g);
           WARNING("Previous value of Hakaru keyword '%1' erased.", g);
         end if;
         protect(g)
    end do
  end proc;

  ModuleUnload := proc($)
    map(proc(x::uneval) try eval(x) catch: NULL; end try end proc,
        ['TypeTools:-RemoveType(`&implies`)'
        ,'TypeTools:-RemoveType(Name)'
        ,'TypeTools:-RemoveType(known_continuous)'
        ,'TypeTools:-RemoveType(known_discrete)'
        ,'TypeTools:-RemoveType(t_Hakaru)'
        ,'TypeTools:-RemoveType(t_pw)'
        ,'TypeTools:-RemoveType(t_partition)'
        ,'TypeTools:-RemoveType(t_pw_or_part)'
        ,'TypeTools:-RemoveType(t_piecewiselike)'
        ,'TypeTools:-RemoveType(t_case)'
        ,'TypeTools:-RemoveType(t_type)'
        ,'TypeTools:-RemoveType(t_lam)'
        ,'VerifyTools[RemoveVerification](measure)'
        ]);
    unprotect(op([2,6], thismodule)); #See comment in ModuleLoad.
    #Skip restoring the globals to any prior value they had.
  end proc;
  ModuleLoad();
end module; # Hakaru

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


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

(* These numbers are taken from the original reference, but divided by
     two to convert from Rydbergs to Hartrees *)

A_vwn  := [ 0.0310907, 0.01554535, -1/(6*Pi^2)]:
b_vwn  := [ 3.72744,   7.06042,    1.13107  ]:
c_vwn  := [12.9352,   18.0578,    13.0045   ]:
x0_vwn := [-0.10498,  -0.32500,   -0.0047584]:

A_rpa  := [ 0.0310907,  0.01554535,  -1/(6*Pi^2)]:
b_rpa  := [13.0720,    20.1231,      1.06835  ]:
c_rpa  := [42.7198,   101.578,      11.4813   ]:
x0_rpa := [-0.409286,  -0.743294,   -0.228344 ]:

Q_vwn  := (b, c)     -> sqrt(4*c - b^2):
f1_vwn := (b, c)     -> 2*b/Q_vwn(b, c):
f2_vwn := (b, c, x0) -> b*x0/(x0^2 + b*x0 + c):
f3_vwn := (b, c, x0) -> 2*(2*x0 + b)/Q_vwn(b, c):

fpp_vwn := 4/(9*(2^(1/3) - 1)):

fx_vwn := (b, c, rs) -> rs + b*sqrt(rs) + c:

f_aux := (A, b, c, x0, rs) -> A*(
  + log(rs/fx_vwn(b, c, rs))
  + (f1_vwn(b, c) - f2_vwn(b, c, x0)*f3_vwn(b, c, x0))*arctan(Q_vwn(b, c)/(2*sqrt(rs) + b))
  - f2_vwn(b, c, x0)*log((sqrt(rs) - x0)^2/fx_vwn(b, c, rs))
):

DMC  := (rs, z) ->
    + f_aux(A_vwn[2], b_vwn[2], c_vwn[2], x0_vwn[2], rs)
    - f_aux(A_vwn[1], b_vwn[1], c_vwn[1], x0_vwn[1], rs):

DRPA := (rs, z) ->
    + f_aux(A_rpa[2], b_rpa[2], c_rpa[2], x0_rpa[2], rs)
    - f_aux(A_rpa[1], b_rpa[1], c_rpa[1], x0_rpa[1], rs):


text:
# https://www.reddit.com/r/adventofcode/comments/ra3f5i/2021_day_6_part_3_day_googol/
# Get the 10^100th generation of Lantern fish

input := "3,1,1,0,3,7,5,5,2,4,2,1,0,2,6,4,3,0,2,1,5,1,4,2,3,0,6,3,0,5,0,5,6,0,0,6,7,0,1,6,3,2,1,1,2,2,0,1,1,1,6,0,2,1,0,5,1,4,7,6,3,0,7,2,0,0,2,0,2,7,3,7,2,4,6,1,6,6,1,1,6,3,3,1,0,4,5,0,5,1,2,0,2,0,7,4,6,1,6,1,5,0,0,2,3":
tally := table(sparse=0,Statistics:-Tally(StringTools:-Split(input,","))):
lanternfish := Vector([seq(tally[cat("",i)],i=0..8)]):

U := Matrix([
 [0, 1, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 1, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 1, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 1, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 1, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 1, 0, 0],
 [1, 0, 0, 0, 0, 0, 0, 1, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 1],
 [1, 0, 0, 0, 0, 0, 0, 0, 0]]):

minpoly := LinearAlgebra:-MinimalPolynomial(U,x);

q:=10^20:
n := 10^100: # googol

# reduce x^n modulo the minpoly to get g so g(U) = U^n
g := Powmod(x, n, minpoly, x) mod q:

add(eval(g, x=U) . lanternfish) mod q;
# 80595178351181834124


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

pkzb_a2 := 146/2025:
pkzb_a3 := -73/405:
pkzb_a4 := 0.131957187845257783631757384393: (* DD + 100/(81*81*kappa) *)

pkzb_kappa := 0.804:

pkzb_xx := (p, qt) -> MU_GE*p + pkzb_a2*qt*qt + pkzb_a3*qt*p + pkzb_a4*p*p:

pkzb_f := (x, u, t) -> 1 + pkzb_kappa - \
  pkzb_kappa^2/(pkzb_kappa + pkzb_xx(X2S^2*x^2, 6*X2S^2*t - 9/20 - X2S^2*x^2/12)):

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


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

mbeef_n := 5:
mbeef_coefs := [
    [ 1.17114923e+00,  1.15594371e-01, -5.32167416e-02, -2.01131648e-02,  1.41417107e-03],
    [-6.76157938e-02,  4.53837246e-02, -2.22650139e-02,  1.92374554e-02,  9.19317034e-07],
    [ 1.48659502e-02,  3.18024096e-02, -5.21818079e-03,  1.33707403e-07, -5.00749348e-07],
    [ 1.40794142e-03, -6.08338264e-03, -6.57949254e-07, -5.49909413e-08,  5.74317889e-08],
    [ 1.41530486e-04, -1.00478906e-07,  2.01895739e-07,  3.97324768e-09, -3.40722258e-09]
]:

$include "mbeef.mpl"

mbeefvdw_f := (x, u, t) -> mbeef_expansion(x, t):

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


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

# Real numbers

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

`is_equal/RR` := (x,y) -> evalb(simplify(y-x) = 0):

`is_leq/RR` := (x,y) -> is(simplify(y-x) >= 0);

`random_element/RR` := proc() rand(-720..720)()/(360); end:

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

`random_element/RR_plus` := proc() rand(1..720)()/(360); end:

######################################################################
# R(N) is the vector space RR^N

`is_element/R` := (N::posint) -> (x) -> type(x,[realcons$N]);

`is_equal/R` := (N::posint) -> (x,y) -> evalb(simplify(x -~ y) = [0$N]);

`is_leq/R` := (N::posint) -> proc(x,y) local i; `and`(seq(`is_leq/RR`(x[i],y[i]),i=1..N)); end:

`is_leq_lex/R` := (N::posint) -> proc(x,y)
 local i,z;
 for i from 1 to n do
  z := simplify(y[i] - x[i]);
  if is(z > 0) then
   return true;
  elif is(z < 0) then
   return false;
  fi;
 od;
 return true;
end:

`is_zero/R` := (N::posint) -> (x) -> evalb(simplify(x) = [0$N]);

`norm_2/R` := (N::posint) -> proc(x) local i; sqrt(add(x[i]^2,i=1..N)); end:

`dot/R` := (N::posint) -> proc(x,y) local i; add(x[i]*y[i],i=1..N); end:

`d_2/R` := (N::posint) -> (x,y) -> `norm_2/R`(N)(x -~ y);

`norm_infinity/R` := (N::posint) -> (x) -> max(map(abs,x));

`d_infinity/R` := (N::posint) -> (x,y) -> `norm_infinity/R`(N)(x -~ y);

`random_element/R` := (N::posint) -> proc() local i; [seq(`random_element/RR`(),i=1..N)] end;

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

`cross_product` := proc(x,y)
 [x[2]*y[3] - x[3]*y[2],
  x[3]*y[1] - x[1]*y[3],
  x[1]*y[2] - x[2]*y[1]];
end:


text:
######################################################################
# Partitions of the set {1,...,n} into intervals

`is_element/interval_partitions` := (n::posint) -> proc(pi)
 local A,U,i;
 global reason;

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

 if not(`is_element/partitions`(A)(pi)) then
  reason := [convert(procname,string),"pi is not a partition",pi,reason];
  return false;
 fi;

 for U in pi do
  if not(`is_element/posint_intervals`(U)) then
   reason := [convert(procname,string),"U is not a posint interval",U];
   return false;
  fi;
 od;

 return true;
end:

`is_equal/interval_partitions` := (n::posint) -> proc(pi,rho)
 local A,i;

 A := {seq(i,i=1..n)};
 return `is_equal/partitions`(A)(pi,rho);
end:

`is_leq/interval_partitions` := (n::posint) -> proc(pi,rho)
 local A,i;

 A := {seq(i,i=1..n)};
 return `is_leq/partitions`(A)(pi,rho);
end:

# Given an integer n > 0 and a subset C of {1/2,3/2,...,(2n-1)/2}, this
# returns the partition of {1,...,n} obtained by cutting at the points
# of C.

`cut_partition` := proc(n::posint,C)
 local m,C0,C1,i,j;

 m := 1 + nops(C);
 C0 := sort([1,op(map(c->c+1/2,C))]);
 C1 := sort([op(map(c->c-1/2,C)),n]);
 return {seq({seq(j,j=C0[i]..C1[i])},i=1..m)};
end;

`random_element/interval_partitions` := (n::posint) -> proc()
 local C,i;

 C := `random_element/subsets`({seq(i+1/2,i=1..n-1)})();
 return cut_partition(n,C);
end;

`list_elements/interval_partitions` := proc(n)
 local G,P,i;

 G := {seq(i+1/2,i=1..n-1)};
 P := `list_elements/subsets`(G);
 return map(C -> cut_partition(n,C),P);
end;

`count_elements/interval_partitions` := (n) -> 2^(n-1);


text:
    # eval_piecewise has the same calling convention as eval_factor.
    # It simplifies piecewise expressions.
    eval_piecewise := proc(e :: specfunc(piecewise),
                           kb :: t_kb, mode :: identical(`*`,`+`),
                           loops :: list([identical(product,Product,sum,Sum),
                                          name=range]),
                           $)
      local default, kbs, pieces, i, cond, inds, res, s, r, x, b, a;
      default := 0; # the catch-all "else" result
      kbs[1] := kb;
      for i from 1 by 2 to nops(e) do
        if i = nops(e) then
          default := op(i,e);
          pieces[i] := default;
        else
          # Simplify piecewise conditions using KB
          cond := op(i,e);
          # cond := eval_factor(cond, kbs[i], `+`, []);
          kbs[i+1] := assert(cond, kbs[i]);
              # This condition is false in context, so delete this piece
              # by not putting anything inside "pieces"
          if kbs[i+1] :: t_kb then
            # find all of the assertions in "kbs[i+1] - kbs[i]"
            cond := map(proc(cond::[identical(assert),anything], $)
                          op(2,cond)
                        end proc,
                        kb_subtract(kbs[i+1], kbs[i]));
            if nops(cond) = 0 then
              default := op(i+1,e);
              pieces[i] := default;
              break;
            else
              cond := `if`(nops(cond)=1, op(1,cond), And(op(cond)));
            end if;
          end if;
          # TODO: Extend KB interface to optimize for
          #       entails(kb,cond) := nops(kb_subtract(assert(cond,kb),kb))=0
          kbs[i+2] := assert(Not(cond), kbs[i]);
          if not(kb_entails(kbs[i], kbs[i+2])) then
            pieces[i] := cond;
            pieces[i+1] := op(i+1,e);
          else
            # This condition is false in context, so delete this piece
            # by not putting anything inside "pieces"
          end if
        end if
      end do;
      # Combine duplicate branches at end
      inds := [indices(pieces, 'nolist', 'indexorder')];
      for i in ListTools:-Reverse(select(type, inds, 'even')) do
        if Testzero(pieces[i]-default) then
          pieces[i  ] := evaln(pieces[i  ]);
          pieces[i-1] := evaln(pieces[i-1]);
        else
          break;
        end if
      end do;
      # Special processing for when the pieces are few
      res := [entries(pieces, 'nolist', 'indexorder')];
      if nops(res) <= 1 then
        return eval_factor(default, kb, mode, loops);
      end if;
      if nops(res) <= 3
         and op(1,res) :: '{`=`,And(specfunc(And),Not(specfunc(Not(`=`),And)))}'
         and Testzero(default - mode()) then
        # Reduce product(piecewise(i=3,f(i),1),i=1..10) to f(3)
        r := op(1,res);
        r := `if`(r::`=`, And(r), select(type,r,`=`));
        for i from 1 to nops(loops) do
          x := op([i,2,1],loops);
          s, r := selectremove(depends, r, x);
          for cond in s do
            if ispoly(lhs(cond)-rhs(cond), 'linear', x, 'b', 'a') then
              b := Normalizer(-b/a);
              if kb_entails(kb,
                            And(b :: integer,
                                op([i,2,2,1],loops) <= b,
                                b <= op([i,2,2,2],loops))
                            ) then
                kb := assert(x=b, kb);# TODO: why not just use kb?
                ASSERT(type(kb,t_kb), "eval_piecewise{product of pw}: not a kb");
                res := `if`(op(1,res) = cond,
                            op(2,res),
                            piecewise(bool_And(op(remove(`=`, op(1,res), cond))),
                                      op(2..-1,res)));
                return eval_factor(eval(res, x=b),
                                   kb,
                                   mode,
                                   eval(subsop(i=NULL, loops), x=b));
              end if;
            end if;
          end do;
          if nops(r) = 0 then break end if;
        end do;
      end if;
      # Recursively process pieces
      inds := [indices(pieces, 'nolist', 'indexorder')];
      for i in inds do
        if i::even or i=op(-1,inds) then
          # only simplify if the piece is not default;
          # note that kbs[i] could be NotAKB(), but this is still valid
          if not Testzero(pieces[i] - default) then
            pieces[i] := eval_factor(pieces[i], kbs[i], mode, []);
          end if;
        end if;
      end do;
      res := piecewise(entries(pieces, 'nolist', 'indexorder'));
      for i in loops do res := op(1,i)(res, op(2,i)) end do;
      return res;
    end proc;


text:
[1025][1115]You always say that.|The same thing every time.
[1118][1177]"l'm through, never again,|too dangerous."
[1180][1270]I know that's what I always say.|I'm always right too.
[1274][1364]- You forget about it in a day or two.|- The days of me forgetting are over.
[1368][1428]The days of me remembering|have just begun.
[1430][1490]You know, when you go on like this,|what you sound like ?
[1493][1583]- I sound like a sensible fuckin' man.|- You sound like a duck.
[1587][1707]- Quack, quack, quack.|- Take heart, 'cause you're never
[1712][1802]Since I'm never gonna|do it again, you're never
[1806][1956]- After tonight ?|- Correct. I got all tonight to quack.
[1962][2022]- Can I get anyone more coffee ?|- Oh, yes !
[2087][2176]- Thank you.|- You're welcome.


text:
with(LinearAlgebra,
     Multiply, Transpose):
with(Groebner,
     Basis, Reduce, InterReduce, LeadingTerm, LeadingMonomial, TestOrder):

writeto("output.txt"):

truncatePolynomial := proc (f, order_1, order_2)
local poly, leading_coeff_1, leading_mon_1, leading_term,
    curr_index, polys:

    poly := expand(f):

    if evalb(type(poly, `+`)) then
        leading_coeff_1, leading_mon_1 := LeadingTerm(poly, order_1):
        leading_term                   := leading_coeff_1 * leading_mon_1:
        curr_index                     := 1:

        polys := sort([op(poly)], (b, a) -> TestOrder(a, b, order_2)):

        while evalb(polys[curr_index] <> leading_term) do
            curr_index := curr_index + 1:
        end do:

        if curr_index = 1 then
            return polys[1]:
        else
            return add(x, x in polys[1..curr_index]):
        end if:
    else
        return poly:
    end if:
end proc:

basisConversion := proc (basis, order_1, order_2)
local F_original, F_sanity_check, F, F_t, num_iter,
    F_t_gb, multipliers,
    G, repeat, curr_index, num_elements:

    F_original := Basis(basis, order_1):

    lprint("Step 1"):
    F   := Basis(basis, order_1):
    F_t := map(v -> truncatePolynomial(v, order_1, order_2), F):
    lprint("Current F", F):
    lprint("Current F_t", F_t):

    num_iter := 1:

    while true do
        lprint("Iteration ", num_iter):
        num_iter := num_iter + 1:

        lprint("Step 2"):
        F_t_gb, multipliers := Basis(F_t, order_2, output=extended):
        lprint("Matrix of multipliers M': ", multipliers):
        lprint("H (= M' F_t): ", F_t_gb):

        lprint("Step 3"):
        G := convert(simplify(
            Multiply(convert(multipliers, Matrix), Transpose(convert(F, Matrix)))), list):
        lprint("G (= M' * F): ", G):

        lprint("Step 4"):
        repeat       := false:
        curr_index   := 1:
        num_elements := numelems(F_t_gb):
        while curr_index < num_elements do
            if evalb(LeadingMonomial(F_t_gb[curr_index], order_2)
                     <> LeadingMonomial(G[curr_index], order_2)) then
                lprint("There are the witness polynomials that prove"\
                       "G_s is not empty (g_j, h_j) respectively ",
                       F_t_gb[curr_index], G[curr_index]):
                repeat := true:
                break:
            end if:
            curr_index := curr_index + 1:
        end do:

        #F   := G:
        F   := InterReduce(G, order_2):

        F_sanity_check := Basis(F, order_1):
        lprint("GB w.r.t. >1 of current F", F_sanity_check):
        lprint("GB w.r.t. >1 of original input", F_original):
        ASSERT(evalb(F_original = F_sanity_check)):

        if repeat = false then
            lprint("Done"):
            return F:
        else
            F_t := map(v -> truncatePolynomial(v, order_1, order_2), F):
            lprint("Current F", F):
            lprint("Current F_t", F_t):
        end if:
    end do:

end proc:

testBasisConversion := proc (basis, order_1, order_2)
local output_1, output_2:
    lprint():
    lprint("Basis conversion of ", basis, " from ", order_1, " to ", order_2):
    output_1 := basisConversion(basis, order_1, order_2):
    output_2 := Basis(basis, order_2):
    lprint("Result obtained using implemented algorithm ", output_1):
    lprint("Result obtained using Maple Basis algorithm ", output_2):
    lprint("Are these results the same? ", evalb(output_1 = output_2)):
end proc:

# ---------------------------------------------------------------------------
# Testing
testBasisConversion([y^2-x,x^2-y*z-1,z^2-x], grlex(x, y, z), plex(x, y, z)):
#testBasisConversion([y^2-x,x^2-y*z-1,z^2-x], grlex(x, y, z), plex(x, z, y)):
#testBasisConversion([y^2-x,x^2-y*z-1,z^2-x], grlex(x, y, z), plex(y, x, z)):
#testBasisConversion([y^2-x,x^2-y*z-1,z^2-x], grlex(x, y, z), plex(y, z, x)):
#testBasisConversion([y^2-x,x^2-y*z-1,z^2-x], grlex(x, y, z), plex(z, x, y)):
#testBasisConversion([y^2-x,x^2-y*z-1,z^2-x], grlex(x, y, z), plex(z, y, x)):
#testBasisConversion([y^2-x,x^2-y*z-1,z^2-x], plex(x, y, z), grlex(x, y, z)):
#testBasisConversion([x*y+z-x*z, x^2-z, 2*x^3-x^2*y*z-1], tdeg(x, y, z), plex(z, y, x)):
#testBasisConversion([x + y + z, x*y + y*z + z*x, x*y*z-1], grlex(x, y, z), plex(x, y, z)):
#testBasisConversion([x + y + z + u, x*y + y*z + z*u + u*x, x*y*z + y*z*u + z*u*x + u*x*y, x*y*z*u-1], grlex(x, y, z, u), plex(x, y, z, u)):
# ---------------------------------------------------------------------------

max_stars_repo_name: typesAreSpaces/Basis-Conversion

max_issues_repo_name: typesAreSpaces/Basis-Conversion

max_forks_repo_name: typesAreSpaces/Basis-Conversion


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

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

$include "mgga_x_br89.mpl"

(* This is only used in B00 *)
params_a_at := 0:

br89_y := Q -> 2*Pi^(2/3)/(3*Q):

(* lower piece *)
pgk_a1 := 1.5255251812009530:
pgk_a2 := 0.4576575543602858:
pgk_a3 := 0.4292036732051034:

pgk_b  := [0.4771976183772063, -1.7799813494556270, 3.8433841862302150,
       -9.5912050880518490, 2.1730180285916720, -30.425133851603660]:

pgk_c  := [0.7566445420735584, -2.6363977871370960, 5.4745159964232880,
       -12.657308127108290, 4.1250584725121360, -30.425133957163840]:

pgk_d  := [0.00004435009886795587, 0.58128653604457910, 66.742764515940610,
       434.26780897229770, 824.7765766052239000, 1657.9652731582120]:

pgk_e  := [0.00003347285060926091, 0.47917931023971350, 62.392268338574240,
       463.14816427938120, 785.2360350104029000, 1657.962968223273000000]:

pgk_UB := 2.085749716493756:

pgk_x_lower := y -> (-arctan(pgk_a1*y + pgk_a2) + pgk_a3) *
            add(pgk_c[i]*y^(i-1), i=1..6)/add(pgk_b[i]*y^(i-1), i=1..6):

pgk_x_upper := y -> (arccsch(pgk_UB*y) + 2) *
            add(pgk_d[i]*y^(i-1), i=1..6)/add(pgk_e[i]*y^(i-1), i=1..6):

br89_mx := Q -> my_piecewise3(
  br89_y(Q) <= 0,
  pgk_x_lower(m_min(br89_y(Q), -br89_min_Q)),
  pgk_x_upper(m_max(br89_y(Q),  br89_min_Q))
):


text:
with(LinearAlgebra):
with(VectorCalculus):
with(plottools);

NALU := proc(w1,w2,ghat1,ghat2,x,t,epsilon)
	local g1, g2, z1, z2, L, i;
	g1 := 1/~(1 +~ exp~(-ghat1));
	z1 := g1 *~ (w1.x) +~ (1 -~ g1) *~ exp~(w1.log~(abs(x) +~ epsilon));

	g2 := 1/~(1 +~ exp~(-ghat2));
	z2 := g2 *~ (w2.z1) +~ (1 -~ g2) *~ exp~(w2.log~(abs(z1) +~ epsilon));

	L := (z2 -~ t)^~2;
	return add(L[i],i=1..numelems(L));
end proc:

NALUsafe := proc(w1,w2,ghat1,ghat2,x,t)
	local g1, g2, z1, z2, L, i;
	g1 := 1/~(1 +~ exp~(-ghat1));
	z1 := g1 *~ (w1.x) +~ (1 -~ g1) *~ exp~(w1.log~(abs(x -~ 1) +~ 1));

	g2 := 1/~(1 +~ exp~(-ghat2));
	z2 := g2 *~ (w2.z1) +~ (1 -~ g2) *~ exp~(w2.log~(abs(z1 -~ 1) +~ 1));

	L := (z2 -~ t)^~2;
	return add(L[i],i=1..numelems(L));
end proc:

NALU(<<w | w>, <w | w>>, <<w | w>>, <g, g>, <g>, <x[1], x[2]>, t, epsilon);

P := plot3d(
	NALU(<<w | w>, <w | w>>, <<w | w>>, <g, g>, <g>, <2, 2>, 8, 10^(-8)),
	w = -1.5..1.5, g = -3..3,
	view = [-1.5..1.5, -3..3, 0..300],
	axes=boxed):
P;

P := plot3d(
	NALU(<<w | w>, <w | w>>, <<w | w>>, <g, g>, <g>, <2, 2>, 16, 10^(-8)),
	w = -1.5..1.5, g = -3..3,
	view = [-1.5..1.5, -3..3, 0..300],
	axes=boxed):
P;

NALUsafe(<<w | w>, <w | w>>, <<w | w>>, <g, g>, <g>, <x[1], x[2]>, t, epsilon);

P := plot3d(
	NALUsafe(<<w | w>, <w | w>>, <<w | w>>, <g, g>, <g>, <2, 2>, 8),
	w = -1.5..1.5, g = -3..3,
	view = [-1.5..1.5, -3..3, 0..300],
	axes=boxed):
P;

P := plot3d(
	NALUsafe(<<w | w>, <w | w>>, <<w | w>>, <g, g>, <g>, <2, 2>, 16),
	w = -1.5..1.5, g = -3..3,
	view = [-1.5..1.5, -3..3, 0..300],
	axes=boxed):
P;

solveNALUsafe := proc(x, t)
	local i, v, w, g, eq, sol1, sol2, sols;

	eq := NALUsafe(<<w | w>, <w | w>>, <<w | w>>, <g, g>, <g>, x, t);
	sols := [];
     for v from -3 to 3 by 0.1 do
       sol1 := fsolve(eval(eq, g=v) = 0, w);
       sol2 := fsolve(eval(eq, g=v) = 0, w, avoid={{w = sol1}});
       sols := [op(sols), [sol1, v], [sol2, v]];
     end do;

	return sols;
end proc:

solveNALU := proc(x, t, epsilon)
	local i, v, w, g, eq, sols, sol1, sol2, sol3, sol4;

	eq := NALU(<<w | w>, <w | w>>, <<w | w>>, <g, g>, <g>, x, t, epsilon);
	sols := [];
     for v from -3 to 3 by 0.1 do
       sol1 := fsolve(eval(eq, g=v) = 0, w);
       sol2 := fsolve(eval(eq, g=v) = 0, w, avoid={{w = sol1}});
       if v < 0.9 and v >= 0 then
         sol3 := fsolve(eval(eq, g=v) = 0, w, avoid={{w = sol1}, {w = sol2}});
         sols := [op(sols), [sol1, v], [sol2, v], [sol3, v]];
       elif v < 0 then
         sol3 := fsolve(eval(eq, g=v) = 0, w, avoid={{w = sol1}, {w = sol2}});
         sol4 := fsolve(eval(eq, g=v) = 0, w, avoid={{w = sol1}, {w = sol2}, {w = sol3}});
         sols := [op(sols), [sol1, v], [sol2, v], [sol3, v], [sol4, v]];
       else
         sols := [op(sols), [sol1, v], [sol2, v]];
       end if;
     end do;

	return sols;
end proc:

NALUsols := solveNALU(<2, 2>, 8, 10^(-8));
P := plot(NALUsols, style = 'point', view = [-1.5..1.5, -3..3]):
P;

NALUsafesols := solveNALUsafe(<2, 2>, 8);
P := plot(NALUsafesols, style = 'point', view = [-1.5..1.5, -3..3]):
P;


text:
PotraShi := module()

  description "linesearch"
              "Efficient Line Search Algorithm for Unconstrained Optimization"
              "F. A. POTRA AND Y.SHI"
              "JOURNAL OF OPTIMIZATION THEORY AND APPLICATIONS"
              "Vol.85,No.3,pp.677-704,1995";

  # Module defined as a package (i.e.) collection of procedures
  option package, load = ModuleLoad, unload = ModuleUnLoad;

  export lines, linesearch_A2, esatta, AB_list, doplot, doplot2;

  local ModuleLoad,
        ModuleUnLoad,
        _debug,
        min_quadratic,
        min_cubic,
        sigma, rho, tau1, tau2, tau3, Jbig,
        max_iter, A_list, B_list;

  uses LinearAlgebra;

  ModuleUnLoad := proc()
    printf("PotraShi unload\n");
    NULL;
  end proc:

  ModuleLoad := proc()
    printf("PotraShi load\n");
    _debug   := true;
    rho      := 0.25;
    sigma    := 0.75;
    tau1     := 0.1;
    tau2     := 0.49;
    tau3     := 2.5;
    Jbig     := 2;
    max_iter := 40;
    NULL;
  end proc;

  # Explicitly call ModuleLoad here so the type is registered when this
  # code is cut&pasted in.  ModuleLoad gets called when the module is
  # read from the repository, but the code used to define a module
  # (like the command below) is not called at library read time.
  ModuleLoad();

  AB_list := proc()
    A_list, B_list;
  end proc;

  min_quadratic := proc ( a, fa, Dfa, b, fb )
    local t, delta; # A+B*(x-a)+C*(x-a)^2
    # p  = fa + Dfa*(x-a)+t*(x-a)^2
    # p' = Dfa+2*t*(x-a)
    ##t     := ((fb-fa)/(b-a)-Dfa)/(b-a);
    ##delta := -Dfa/(2*t);
    t     := (fb-fa)/(b-a)-Dfa;
    delta := max(tau1,min(1-tau1,-Dfa/(2*t)));
    a+(b-a)*delta;
  end proc;

  lines := proc( f0, df0, L )
    [ [0,f0], [L,f0+L*rho*df0] ], [ [0,f0], [L,f0+L*sigma*df0] ];
  end proc;

  linesearch_A2 := proc( f_in, df_in )
    local f, df, a, b, c, alpha,
          fa, dfa, fb, fc,
          r_df0, s_df0,
          f0, df0, df1, t1, t2, delta, k;

    f  := x->evalf(f_in(x));
    df := x->evalf(df_in(x));

    alpha  := 1;
    a      := 0;
    b      := 1;
    fa     := f(0);
    fb     := f(1);
    df0    := df(0);

    if df0 >= 0 then
      error "[linesearch_A2] df0 = %1 must be < 0!\n", df0;
    end;

    r_df0  := rho*df0;
    s_df0  := sigma*df0;
    f0     := fa;

    A_list := [];
    B_list := [];

    # step 1 (check alpha = 1)
    if fb <= f0 + r_df0 then
      if sigma > 0.5 then
        if fb >= f0 + s_df0 then
          if _debug then
            printf("[step1] return 1 (a case)\n");
          end;
          return 1;
        end;
      else
        df1 := df(1);
        if df1 >= s_df0 then
          if _debug then
            printf("[step1] return 1 (b case)\n");
          end;
          return 1;
        end;
      end;
      # step 2
      if _debug then
        if sigma > 0.5 then
          printf(
            "\n[step2] now\n"
            "f(1) [%g] <= f(0)+rho*f'(0) [%g]\n"
            "and f(1) [%g] < f(0)+sigma*f'(0) [%g]\n",
            fb, f0 + r_df0,
            fb, f0 + s_df0
          );
        else
          printf(
            "\n[step2] now\n"
            "f(1) [%g] <= f(0)+rho*f'(0) [%g]\n"
            "and f'(1) [%g] < sigma*f'(0) [%g]\n",
            fb, f0 + r_df0,
            df1, s_df0
          );
        end;
      end;
      a      := 1;
      b      := Jbig;
      A_list := [op(A_list), a];
      B_list := [op(B_list), b];
      fa := fb;
      fb := f(b);
      while fb <= fa + (b-a)*r_df0 do

        if fb >= fa + (b-a)*s_df0 then
          if _debug then
            printf( "[step2] return b = %g\n", b );
          end;
          return b;
        end;

        a  := b;
        b  := Jbig*b;
        fa := fb;
        fb := f(b);

        if _debug then
          printf( "[step2] a=%g b=%g f(a)=%g f(b)=%g\n", a, b, fa, fb );
        end;

        A_list := [op(A_list), a];
        B_list := [op(B_list), b];

        if a > 1000 then
          error "a too big\n";
        end;
      end;
    else
      printf(
        "\n[step2] case\n"
        "f(1) [%g] > f(0)+rho*f'(0) [%g]\n",
        fb, f0 + r_df0
      );
    end;

    # aggiunta ad algoritmo
    if a > 0 then
      dfa := df(a);
      if dfa >= 0 then
        if _debug then
          printf("[step3] return a = %g (added)\n", a);
        end;
        return a;
      end;
    else
      dfa := df0;
    end;

    # step 3
    if _debug then
      printf(
        "\n[step3] a=%g b=%g with f(0)=%g, f(a)=%g, f'(a)=%g\n"
        "f(a) [%g] <= f(0)+a*rho*f'(0) [%g]\n"
        "and f(b) [%g] > f(a)+(b-a)*sigma*f'(0) [%g]\n",
        a, b, f0, fa, dfa,
        fa, f0 + a*r_df0,
        fb, fa + (b-a)*s_df0
      );
    end;

    for k from 1 to max_iter do

      c := min_quadratic( a, fa, dfa, b, fb );
      #c  := (a+b)/2;
      fc := f(c);

      if fc <= fa + (b-a)*r_df0 and fc >= fa + (b-a)*s_df0 then
        if _debug then
          printf("[step3] return c = %g (case a)\n", c);
        end;
        return c;
      end;

      t1    := (fb-fc)/(b-c);
      t2    := (fc-fa)/(c-a);
      delta := abs( (t1-t2)/(b-a) );

      if fc <= fa + (c-a)*r_df0 then
        if (rho-sigma)*df0 >= tau3*(b-a)*delta then
          if _debug then
            printf("[step3] return c = %g (case b)\n", c);
          end;
          return c;
        end;
        a   := c;
        fa  := fc;
        dfa := df(c);
      else
        if (rho-sigma)*df0 >= tau3*(b-a)*delta and a > 0 then
          if _debug then
            printf("[step3] return a = %g\n", a);
          end;
          return a;
        end;
        b  := c;
        fb := fc;
      end;

      A_list := [op(A_list), a];
      B_list := [op(B_list), b];

      if _debug then
        printf(
          "[step3] a=%.4g b=%.4g f(a)=%.4g f'(a)=%.4g f(b)=%.4g f'(b)=%.4g\n",
          a, b, fa, dfa, fb, df(b)
        );
      end;

      if a > 1000 then
        error "a = %1\n", a;
      end;
    end;

    printf("[step3] search failed\n");
    return NULL;

  end proc;

  esatta := proc( fun )
    local res, alpha;
    #res := Optimization[Minimize](
    res := Optimization[NLPSolve](
      fun(alpha),
      {alpha >= 0, alpha <= 3},
      #initialpoint={alpha=1},
      iterationlimit=50
    ):
    subs(res[2],alpha);
  end proc:

  doplot := proc( f_in, df_in, alpha, amax)
    local f0, df0, x, AA, BB, CC, DD;
    f0  := evalf(f_in(0));
    df0 := evalf(df_in(0));

    AA := plot(f_in(x),x=0..amax);
    BB := plot([[0,f0],[amax,f0+amax*rho*df0]],color="LimeGreen");
    CC := plot([[0,f0],[amax,f0+amax*sigma*df0]],color="blue");
    DD := plot([[alpha,f_in(alpha)]],color="red",style = point);
    display(AA,BB,CC,DD);
  end;

  doplot2 := proc( f_in, df_in, alpha, amax, A_list, B_list )
    local i, f0, df0, x, AA, BB, CC, DD, EE;
    f0  := evalf(f_in(0));
    df0 := evalf(df_in(0));

    AA := plot(f_in(x),x=0..amax);
    BB := plot([[0,f0],[amax,f0+amax*rho*df0]],color="LimeGreen");
    CC := plot([[0,f0],[amax,f0+amax*sigma*df0]],color="blue");
    DD := plot([[alpha,f_in(alpha)]],color="red",style = point);
    EE := plot([seq([[A_list[i],-i/5],[B_list[i],-i/5]],i=1..nops(A_list))],color="black");
    display(AA,BB,CC,DD,EE);
  end;

end module:

max_stars_repo_name: ebertolazzi/NLtoolbox

max_issues_repo_name: ebertolazzi/NLtoolbox

max_forks_repo_name: ebertolazzi/NLtoolbox


text:
 func $band32I (
  var %i i32
  ) i32 {
   return (
     band i32(dread i32 %i,
       constval i32 0x12))}
 func $band32I_2 (
  var %i i32
  ) i32 {
   return (
     band i32(dread i32 %i,
       constval i32 0x112))}
 func $band32RR (
  var %i i32, var %j i32
  ) i32 {
   return (
     band i32(dread i32 %i,
       dread i32 %j))}

 func $band64I (
  var %i i64
  ) i64 {
   return (
     band i64(dread i64 %i,
       constval i64 0x1200000015))}
 func $band64I_2 (
  var %i i64
  ) i64 {
   return (
     band i64(dread i64 %i,
       constval i64 0x1200000115))}
 func $band64RR (
  var %i i64,
  var %j i64
  ) i64 {
   return (
     band i64(dread i64 %i,
              dread i64 %j))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


text:
fileinfo{
  @magic 345,
  @filename "adir/thisisafile.dex",
  @dexoffset 560,
  @flags 0x01ef,
  @linksize 100
}
type $SS <class {
        @real f32,
        @imag f32}>
var $sconst1 <class {
        @classidx 22,
        @f11 i32,
        @f22 <$SS>,
        @f33 i32,
        @accessflag 0x45e,
        @classdataoffset 118,
        @f44 i32,
        @f55 f32,
        @f66 f32,
        @superclassid 3333
           }>
func $fact () i32 {
  funcinfo{
    @funcname "afunc",
    @registers 345,
    @ins 560,
    @outs 444
  }
  return (dread i32 $sconst1 1)
  funcinfo{
    @funcidx 345,
    @classidx 560
  }
}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
restart;
with(LinearAlgebra):
with(VectorCalculus):
with(plottools):
BasisFormat(false):

MNAC := proc(w1,w2,x,t,epsilon)
	local z1, z2, L, i;
	z1 := w1.x;
	z2 := exp~(w2.log~(abs(z1) +~ epsilon));

	L := (z2 -~ t)^~2;
	return add(L[i],i=1..numelems(L));
end proc:

MNACalt := proc(w1,w2,x,t)
	local z1, z2, zhat2, w2ones, L, i;

	z1 := w1.x;
	zhat2 := log~(Transpose(w2) *~ z1 +~ 1 -~ w2);
	z2 := exp~(zhat2 . Matrix(ColumnDimension(w2),1,1));
	L := (z2 -~ t)^~2;
	return add(L[i],i=1..numelems(L));
end proc:

P := plot3d(
	MNAC(<<w[1]|w[1]|0|0>, <w[1]|w[1]|w[1]|w[1]>>, <<w[2]|w[2]>>, <1,1.2,1.8,2>, 13.2, 1E-7),
	w[1] = -1 .. 1, w[2] = -1 .. 1, view = [-1..1, -1..1, 0..500],
	axes=boxed, orientation=[-45, 39, 0]):
Export(FileTools:-JoinPath([currentdir(), "..", "paper", "graphics", "nac-mul-eps-1em7.jpeg"]), P);

P := plot3d(
	MNAC(<<w[1]|w[1]|0|0>, <w[1]|w[1]|w[1]|w[1]>>, <<w[2]|w[2]>>, <1,1.2,1.8,2>, 13.2, 1E-1),
	w[1] = -1 .. 1, w[2] = -1 .. 1, view = [-1..1, -1..1, 0..500],
	axes=boxed, orientation=[-45, 39, 0]):
Export(FileTools:-JoinPath([currentdir(), "..", "paper", "graphics", "nac-mul-eps-1em1.jpeg"]), P);

P := plot3d(
	MNAC(<<w[1]|w[1]|0|0>, <w[1]|w[1]|w[1]|w[1]>>, <<w[2]|w[2]>>, <1,1.2,1.8,2>, 13.2, 1),
	w[1] = -1 .. 1, w[2] = -1 .. 1, view = [-1..1, -1..1, 0..500],
	axes=boxed, orientation=[-45, 39, 0]):
Export(FileTools:-JoinPath([currentdir(), "..", "paper", "graphics", "nac-mul-eps-1.jpeg"]), P);

P := plot3d(
	MNACalt(<<w[1]|w[1]|0|0>, <w[1]|w[1]|w[1]|w[1]>>, <w[2]|w[2]>, <1,1.2,1.8,2>, 13.2),
	w[1] = -1 .. 1, w[2] = 0 .. 1, view = [-1..1, 0..1, 0..500],
	axes=boxed, orientation=[-45, 39, 0]):
Export(FileTools:-JoinPath([currentdir(), "..", "paper", "graphics", "nac-mul-nmu.jpeg"]), P);


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

$define lda_c_pw_params
$include "lda_c_pw.mpl"

(* Parameters from Table 6 *)
gap_par0 := [
   0.04953, 1.07924, 0.07928, (* a1, a2, a3 *)
  -2.504e-2, 7.026e-3, -1.268e-3, 1.136e-4, -3.841e-6, (* b3, b4, b5, b6, b7 *)
   0.031091,  (* params_a_a of lda_c_pw *)
   0.23878    (* pre-factor of C = 0.06483*((9*Pi)/4)^(2/3) *)
]:
gap_par1 := [
   0.0471985, 1.49676, 0.00179054,
  -3.24091e-2, 9.99978e-3, -1.93483e-3, 1.79118e-4, -6.15798e-6,
   0.015545,
   0.064535
]:

(* Equation (20): e'(rs) *)
gap_eps_1 := (rs, par) ->
  par[1]*rs^(3/2)/(1 + sqrt(rs)*(par[2] + par[3]*sqrt(rs) + par[1]*rs)):
(* Equation (21): e''(rs) *)
gap_eps_2 := (rs, par) ->
  add(par[i+1]*rs^i, i=3..7):

(* Equation (19) *)
gap_C := (rs, par) -> par[10]/rs^2:

(* Equation (17) *)
gap_c2 := (rs, z, par) ->
  + (2*f_pw(rs, z)*gap_eps_1(rs, par) - gap_C(rs, par)*gap_eps_2(rs, par))
  / (2*(gap_C(rs, par)*gap_eps_1(rs, par) - f_pw(rs, z)^2)):
(* Equation (18) *)
gap_c3 := (rs, z, par) ->
  - (2*gap_eps_1(rs, par)^2 - f_pw(rs, z)*gap_eps_2(rs, par))
  / (2*(gap_C(rs, par)*gap_eps_1(rs, par) - f_pw(rs, z)^2)):
(* Equation (16) *)
gap_c1 := (rs, z, par) ->
  - gap_C(rs, par) * gap_c3(rs, z, par):

(* after Equation (6): a = 30 is a parameter fixed by minimizing the
variance of the correlation energy error for the noble gas atoms He,
Ne, and Ar *)
gap_par_a := 30:

(* Equation (6) *)
gap_H := (rs, t, par) ->
  (gap_par_a  + par[9]*rs*log(rs)*t^2/beta_Hu_Langreth(rs))/(gap_par_a + t^2):

gap_t := (rs, z, xt) ->
  xt*n_total(rs)^(1/6)/(4*mphi(z)*(3/Pi)^(1/6)):

(* Gap function, Equation (5) *)
gap_G := (rs, z, xt, par) ->
  + mphi(z)^3*beta_Hu_Langreth(rs)*gap_t(rs, z, xt)^2
  * gap_H(rs, gap_t(rs, z, xt), par)
  / (gap_c1(rs, z, par) - gap_c2(rs, z, par)*f_pw(rs, z)):

(* Correlation energy per particle, Equation (4) *)
gap_eps := (rs, z, xt, par) ->
  + (f_pw(rs, z) + gap_c1(rs, z, par)*gap_G(rs, z, xt, par))
  / (1 + gap_c2(rs, z, par)*gap_G(rs, z, xt, par) + gap_c3(rs, z, par)*gap_G(rs, z, xt, par)^2):

(* Total energy, Equation (2) *)
f_gap := (rs, z, xt) ->
  + gap_eps(rs, 0, xt, gap_par0)
  + f_zeta(z)*(gap_eps(rs, 1, xt, gap_par1) - gap_eps(rs, 0, xt, gap_par0)):

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


text:
######################################################################
# nonempty_subsets(A) is the set of all nonempty subsets of A

`is_element/nonempty_subsets` := (A::set) -> proc(B)
 type(B,set) and B minus A = {} and nops(B) > 0;
end;

`is_equal/nonempty_subsets` := (A::set) -> (B,C) -> evalb(B = C):

`is_leq/nonempty_subsets` := (A::set) -> (B,C) -> evalb(B minus C = {}):

`random_element/nonempty_subsets` := (A::set) -> proc()
 local r,B;

 if nops(A) = 0 then return FAIL; fi;

 r := rand(2);
 B := {};

 while nops(B) = 0 do
  B := select(a -> (r() = 1),A);
 od;

 return B;
end;

`list_elements/nonempty_subsets` := (A::set) ->
 sort(map(sort,[op(combinat[powerset](A) minus {{}})])):

`count_elements/nonempty_subsets` := (A::set) -> 2^nops(A) - 1;


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

`is_element/height_functions` := (A::set) -> proc(h)
 local P,U,V,a,u;
 global reason;

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

 P := {op(`list_elements/nonempty_subsets`(A))};
 if map(op,{indices(h)}) <> P then
  reason := [convert(procname,string),"h is not indexed by the nonempty subsets of A",h,A];
  return false;
 fi;

 for a in A do if h[{a}] <> 0 then
  reason := [convert(procname,string),"h({a}) <> 0",h,a];
  return false;
 fi; od;

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

 for U in P do
  for V in P do
   if U intersect V <> {} and
      h[U union V] <> max(h[U],h[V]) then
    reason := [convert(procname,string),"h(U u V) <> max(h(U),h(V))",U,V,h[U],h[V],h[U union V]];
    return false;
   fi;
  od;
 od;

 return true;
end;

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

`is_equal/height_functions` := (A::set) -> proc(h1,h2)
 local P,U;
 P := `list_elements/nonempty_subsets`(A);
 for U in P do
  if h1[U] <> h2[U] then return false; fi;
 od;
 return true;
end;

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

`is_leq/height_functions` := NULL;

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

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

`is_critical/height_functions` := (A::set) -> proc(h,U)
 local V,a;

 V := A minus U;

 if V = {} then return evalb(h[A] < 1); fi;

 for a in V do
  if h[U union {a}] = h[U] then return false; fi;
 od;

 return true;
end;

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

`critical_tree/height_functions` := (A::set) -> proc(h)
 local P;

 P := `list_elements/nonempty_subsets`(A);

 return {op(select(U ->`is_critical/height_functions`(A)(h,U),P))};
end;

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

`theta/partitions/height_functions` := (A::set) -> proc(pi)
 local h,P,T,B,Q;

 h := table();
 P := `list_elements/nonempty_subsets`(A);
 for T in P do h[T] := 1; od;
 for B in pi do
  Q := `list_elements/nonempty_subsets`(B);
  for T in Q do h[T] := 0; od;
 od;

 return eval(h);
end;

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

`theta/realisation/partitions/height_functions` := (A::set) -> proc(x)
 local P,T,C,pi,t,k,h;

 h := table();
 P := `list_elements/nonempty_subsets`(A);
 for T in P do h[T] := 0; od;

 C := map(op,[indices(x)]);

 for pi in C do
  t := x[pi];
  k := `theta/partitions`(A)(pi);
  for T in P do
   h[T] := h[T] + t * k[T];
  od;
 od;

 return eval(h);
end:

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

`theta_inv/realisation/partitions/height_functions` := (A::set) -> proc(h)
 local TT,TTT,S,r,s,t,x,U,X,Y,pi,a,m,i;

 TT := `critical_tree/height_functions`(A)(h);
 S := map(T -> h[T],TT) minus {1};
 r := nops(S) - 1;
 s := table();
 s[-1] := 1;
 for i from 0 to r do s[i] := S[r+1-i]; od;
 TTT := table();
 x := table();
 for i from 0 to r do
  TTT[i] := select(U -> h[U] <= s[i],TT);
  X := TTT[i];
  pi[i] := {};
  while X <> {} do
   a := X[1][1];
   Y := select(T -> member(a,T),X);
   m := max(map(nops,Y));
   U := select(T -> nops(T) = m,Y)[1];
   pi[i] := {op(pi[i]),U};
   X := X minus Y;
  od;
  x[pi[i]] := s[i-1] - s[i];
 od;

 return eval(x);
end:


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

lines := subs(""=NULL,(op@StringTools:-Split)~(StringTools:-Split(input, "\n"),"|"));

# part 1, just count the 1, 4, 7, and 8s
total := 0:
for i from 2 to nops(lines) by 2 do
    tal := table(sparse=0,Statistics:-Tally(map(length, StringTools:-Split(lines[i]))));
    total := total + tal[2] + tal[3] + tal[4] + tal[7];
end do:

answer1 := total;

# part 2

Decode := proc( code, output )

    local decode := table();
    # inputs are sets of letters so we can use set operations on everything
    local displays := { ({op@StringTools:-Explode})~(
        StringTools:-Split(StringTools:-Trim(code)," "))[] };

    # easy cases - unique lengths
    decode["1"] := select(d->nops(d)=2, displays)[1];
    decode["4"] := select(d->nops(d)=4, displays)[1];
    decode["7"] := select(d->nops(d)=3, displays)[1];
    decode["8"] := select(d->nops(d)=7, displays)[1];

    local sixes := select(d->nops(d)=6, displays);
    local fives := select(d->nops(d)=5, displays);

    # length six 0, 6, 9
    # 6 is disjoint from 1
    decode["6"] := select(d->not decode["1"] subset d, sixes)[1];
    # 9 contains 4
    decode["9"] := select(d->decode["4"] subset d, sixes)[1];
    # otherwise 0 - 0 contains 7 also works
    decode["0"] := op(1,sixes minus {decode["6"], decode["9"]});

    # length five 2,3,5
    # 3 contains 1
    decode["3"] := select(d->decode["1"] subset d, fives)[1];
    # 5 is contained in 6
    decode["5"] := select(d->d subset decode["6"], fives)[1];
    # otherwise 2
    decode["2"] := op(1,fives minus {decode["5"], decode["3"]});

    local decodsubs := sort((rhs=lhs)~([entries(decode,pairs)]));

    local outs := ({op@StringTools:-Explode})~(
        StringTools:-Split(StringTools:-Trim(output)," "));

    return parse(StringTools:-Join(subs(decodsubs, outs),""));

end proc: # Decode

answer2 := add( Decode(lines[i], lines[i+1]), i=1..nops(lines), 2);