text:
BinomSums := module()

option package;

local
  gfdict, gfdictcomp, gfnames, multinomial, packvars, geomsum, solvecons,
  linearorder, slopes, asylt, inorout_base, ratres2, geomsum_, simpfacts,
  simpfacts0, SimpleImpl, isconvergent, asy, inorout, sumtores0 ;

export
  sumtores, addnewgf, rser, computesum, hermitered, ratres, geomred,
  geomredall, sumtoct, `BinomSums/version`;

  `BinomSums/version` := 0.11;

  # This module provides a simple implementation of the computation of integral
  # representation of binomial sums.
  SimpleImpl := module()

    export sumtores;

    # Input :
    #   - U a binomial sum, as a Maple expression
    #   - v a name
    #   - num an integer (default: 1)
    # Output :
    #   A rational function R(v[num],...,v[num+k]) such that U = res R.
    sumtores := proc(U, v :: name, num :: integer := 1)
      local L, first, rest, rat_first, num_first, i;
      if type(U, `+`) then
        return normal(map(sumtores, U, v, num));
      elif type(U, `*`) or type(U, `^`(anything, posint)) then
        first := op(1, U);   # the first factor
        if type(U, `*`) then
          rest := subsop(1=1, U);  # the rest of the product
        else
          rest := first^(op(2, U) - 1);
        end if;
        rat_first := sumtores(first, v, num);
        num_first := max(num, op(map(op, indets(rat_first, specindex(v)))));
        return normal(rat_first*sumtores(rest, v, num_first + 1));
      elif type(U, specfunc(Delta)) then
        return v[num]^op(U);
      elif type(U, specfunc(Binomial)) then
        return (1+v[num])^op(1, U)/v[num]^(op(2, U)+1);
      elif type(U, specfunc(Sum)) then
        return normal(sum(expand(sumtores(op(1, U), v, num)), op(2, U)));
      else
        return U/v[num];
      end if;
    end proc;

  end module:

#### SUM TO PERIOD

gfdict := table();
gfdictcomp := table();
gfnames := {};

# Input :
#   - name, a name
#   - a function gf such that name(n) stands for res( gf(n) )
#   - a function comp such that comp( n ) = res( gf(n) )
#
# Effect :
#   Add 'name' as a basic block for constructing binomial sums
addnewgf := proc(name :: name, gf, comp)
  #global gfdict, gfdictcomp, gfnames;
  gfdict[name] := gf;
  gfdictcomp[name] := comp;
  gfnames := gfnames union {name};
end proc;

# There are different definitions of the binomial.
# Binomial(n,k) = coeff of x^k in the power series (1+x)^n
# It is non-zero only if k ≥ 0.
addnewgf(Binomial,
  ((v,n,k) -> (1+v[1])^n/v[1]^(k)),
  (n,k) -> `if`(k>=0, binomial(n,k), 0));

# binomial2(n, k) = binomial(n, n-k)
# If Binomial2(n,k) and Binomial(n,k) are both non zero, their value coincide.
# Maple's binomial evaluated at (n,k) is the non-zero element of {Binomial(n,k),Binomial2(n,k)}, if any.
addnewgf(Binomial2,
  ((v,n,k) -> (1+v[1])^n*v[1]^k/v[1]^n),
  (n,k) -> `if`(n-k>=0, binomial(n,n-k), 0));

# Binomial3(n,k) = Binomial2(n,k)
# alternative definition.
addnewgf(Binomial3,
  ((v,n,k) -> 1/(1-v[1])^(k+1)/v[1]^(n-k)),
  (n,k) -> `if`(n-k>=0, binomial(n,n-k), 0));

# natbinomial(n,k) = binomial(n,k) if n ≥ 0 and k ≥ 0 and 0 otherwise.
addnewgf(NatBinomial,
  ((v,n,k) -> 1/(1-v[1]-v[2])/v[1]^k/v[2]^(n-k)),
  (n,k) -> `if`(n >=0 and k >= 0, binomial(n,k), 0));

addnewgf(Multinomial,
  ((v,L) -> 1/(1-add(v[i],i=1..nops(L)))/mul(v[i]^(L[i]),i=1..nops(L))),
  (L ->multinomial(L)));

addnewgf(Catalan,
  ((v,n) -> (1+v[1])^(2*n)*(1-v[1])/v[1]^(n)),
  (n -> 1/(n+1)*binomial(2*n,n) ));

addnewgf(H,
  ((v,n) -> v[1]^(-n)/(1-v[1])),
  n -> `if`(n >= 0, 1, 0));

addnewgf(Delta,
  ((v,n) -> v[1]^n),
  n -> `if`(n = 0, 1, 0));

#addnewgf(CT,
#  ( (R, v) -> CT(R,v) ),
#  (R, v) -> residue(R/v,v=0));

multinomial := proc( x :: list(numeric))
  if nops(x) <= 1 then
    if x[1] < 0 then
      return 0;
    else
      return 1;
    end if;
  elif x[1] < 0 then
    return 0;
  else
    return binomial(`+`(op(x)), x[1])*multinomial(x[2..-1]);
  end if;
end;

# Outputs a list of sets s[1],...,s[nops(L)] such that the substitutions
# subs(s[i], L[i]) guarantees that the variables of the elements of L does not
# overlap and and are numbered consecutively.
#
# Example : packvars([ 1/(x[4]+x[6]), x[7] ], x);
#    => [{x[4] = x[1], x[6] = x[2]}, {x[7] = x[3]}]
packvars := proc(L :: list, v :: name)
  local Lv, Ls, i, j, l;
  Lv := map2(select, has, map(indets, L, name), v);
  Ls := NULL;
  i := 0;
  for l in Lv do
    Ls := Ls, {seq(l[j]=v[i+j], j=1..nops(l))};
    i := i + nops(l);
  end do;
  return [Ls];
end proc;

#
# (Alternative implementation)
geomsum_ := proc(S, bounds)
  local T, svar, infb, supb, ret, prim:
  T := expand(normal(S));
  if type(T, `+`) then
    return normal(map(geomsum_, T, bounds));
  fi;

  svar := op(1, bounds);
  infb := op([2, 1], bounds);
  supb := op([2, 2], bounds);

  ASSERT(infb <> -infinity);
  if supb = infinity then
    T := _W[_||svar]^svar*T;
  fi;

  prim := SumTools[Hypergeometric][Gosper](T, svar);
  ret := -subs(svar=infb, prim);

  if supb <> infinity then
    ret := ret + subs(svar=supb+1, prim);
  fi;

  return normal(ret);
end proc;

# Input :
#   - S, an expression of the form P(k)*A^k, where P is a polynomial
#   - bounds, an expression in the form k=a..b, where b can be infinity
#
# Returns an expression T without k such that
#   T = sum(S(k), k=a..b)
#
# If the upper bound b is `infinity' then it marks the sum
# with the extra variable _W[_k].
geomsum := proc(S, bounds)
  local T, svar, infb, supb, ret, prim, den, opT, denother, dendep:

  svar := op(1, bounds);
  infb := op([2, 1], bounds);
  supb := op([2, 2], bounds);

  T := normal(S);

  if infb = -infinity or supb = -infinity then
    error "Summation bounds cannot be -infinity.";
  fi;

  if supb = infinity then
    T := _W[_||svar]^svar*T;
  fi;

  # dendep (resp. denother) contains the factors of the denominator of T that
  # depends (resp. do not depend) on the summation variable.
  den := factor(denom(T));
  if type(den, `*`) then
    dendep, denother := selectremove(has, den, svar);
  elif has(den, svar) then
    dendep, denother := den, 1;
  else
    dendep, denother := 1, den;
  end if;

  T := expand(numer(T));
  if type(T, `+`) then
    opT := [op(T)];
  else
    opT := [T];
  end if;
  opT := map(`*`, opT, 1/dendep);

  # Factor simplification, to rewrite things like
  # A := -(-1+u[1])^(-n-1+j)*(-1+u[2])^(-n+j-1)*(1-u[1])^(-j)*(1-u[2])^(-j)*u[1]^(-n)*u[2]^(-n)
  # which actually does not depend on j.
  # On such an input, SumTools[IndefiniteSummation](A, j) does not work.
  opT := map(simpfacts, opT);
  opT := map(SumTools[IndefiniteSummation], opT, svar);

  ret := -convert(subs(svar=infb, opT),`+`);

  if supb <> infinity then
    ret := ret + convert(subs(svar=supb+1, opT), `+`);
  fi;

  return (ret/denother);
end proc;

# Input:
#   - P, a product of expressions of the form F^a, where `F' is a polynomial
#   and `a' an expression that may depend on variables.
#
# Returns an equivalent product (when the variables in the exponent are
# integers) with the garantee that for any two factors F^a and G^b, F is not
# proportional to G.
simpfacts := proc(P)
  local t, i, ret;
  t := simpfacts0(P);
  ret := 1;
  for i in indices(t) do
    if op(i) = -1 then
      ret := ret*op(i)^(t[op(i)] mod 2);
    else
      ret := ret*op(i)^t[op(i)];
    end if;
  end do;
end proc;

# Input:
#   - same input as simpfacts
#
# Return an associative array T such that
# P = mul(F^a, (F, a) in T)
simpfacts0 := proc(P)
  local t, t1, ind, f, i, lc;
  if type(P, `^`) then
    t := simpfacts0(op(1, P));
    for f in indices(t) do
      t[op(f)] := t[op(f)]*op(2, P);
    end do;
    return t;
  elif type(P, `*`) then
    t := table();
    ind := {};
    for f in [op(P)] do
      t1 := simpfacts0(f);
      for i in indices(t1) do
        if not (i in ind) then
          ind := ind union {i};
          t[op(i)] := 0;
        end if;
        t[op(i)] := t[op(i)] + t1[op(i)];
      end do;
    end do;
    return t;
  else
    f := normal(expand(P));
    lc := lcoeff(f);
    f := normal(expand(f/lc));
    return table([ lc = 1, f = 1 ]);
  end if;
end proc;

# Given a binomial sum, returns an expression such that the binomial sum equal
# the constant term of the expression.
#
# The output may contains extravariables _W[_k] to track infinite summations
sumtoct := proc(S, v :: name)
  local L;
  #global gfdict, gfnames;
  if type(S, specfunc(Sum)) then
    return geomsum(sumtoct(op(1, S), v), op(2, S));
  elif type(S, `+`) then
    return map(sumtoct, S, v);
  elif type(S, `*`) then
    L := map(sumtoct, convert(S, list), v);
    return convert(zip(subs, packvars(L, v), L), `*`);
  elif type(S, `^`(anything, posint)) then
    L := [sumtoct(op(1, S), v) $ op(2, S)];
    return convert(zip(subs, packvars(L, v), L), `*`);
  elif type(S, specfunc(CT)) then
    return subs(op(2,S)=v[1], op(1, S));
  elif type(S, specfunc(gfnames)) then
    return eval(gfdict[op(0,S)](v, op(S)));
  else
    return S;
  end if;
end proc;

#### INFINITE SUMS

# Input :
#   - cons, a set of Laurent monomials
#   - G, a directed graph whose vertices are variables and edges are domination
#   relation.
#   - params, set of `small' variables
#
# If possible, returns a directed acyclic graph H extending G such that if the
# variables are ordered according to H, then every monomials in cons is greater
# than 1 (lexicographic ordering).
#
# Raises an error if not possible.
solvecons := proc(cons :: set, G := false, params := {})
  local cons1, rcons, u, v, H;
  uses GraphTheory;

  if G = false then
    H := GraphTheory[Digraph]({seq([1, v], v in indets(cons)), seq(seq([v, u], u in params), v in indets(cons) minus params)});
    return solvecons(cons, H);
  elif not IsAcyclic(G) then error "inconsistent";
  elif nops(cons) = 0 then return G;
  end if;

  cons1 := cons[1];
  rcons := cons minus {cons1};
  for v in indets(numer(cons1)) do
    try
      H := CopyGraph(G);
      AddArc(H, {seq([w, v], w in indets(denom(cons1)))});
      return solvecons(rcons, H);
    catch "inconsistent" :
    end;
  end do;

  error "inconsistent";
end proc;

# Returns a linear order compatible with the DAG G.
linearorder := proc(G)
  local racines, H, rest;
  uses GraphTheory;

  if nops(Vertices(G)) = 0 then
    return [];
  end if;

  racines := select(v -> InDegree(G, v) = 0, Vertices(G));
  H := DeleteVertex(G, racines);
  rest := linearorder(H);

  return [op(rest),op(racines)];
end proc;

# Input :
#   - R, a rational function
#   - params, a list of variables
#
# Assume that R is a Laurent formal series w.r.t. variables _W[1],...,_W[r]
# Let T = sum of coefficients of this series
#
# If there exist an order such that T is convergent, then returns T and that order.
# If not, fails
#
# params is an indication on the way the parameters should be ordered.
isconvergent := proc(R :: ratpoly, params :: list(name) := [])
  local svars, L, ct, co, mord, res, cons, G, ord, den, facts, f;

  svars := indets(R, specindex(_W));
  L := [ op(params), op(indets(R) minus svars minus convert(params, set)) ];
  mord :=  plex(op(L));
  cons := {};

  facts := map2(op, 1, factors(denom(normal(R)))[2]);
  for f in facts do
    ct := subs(map(`=`, svars, 0), f);
    co := remove(`=`, [coeffs(collect(normal(f/ct-1), svars, distributed, normal), svars)], 0);
    cons := cons union {seq( Groebner[TrailingTerm](numer(c), mord)[2]/Groebner[TrailingTerm](denom(c), mord)[2], c in co)};
  end do;

  G := solvecons(cons, false, {});
  ord := remove(`=`, linearorder(G), 1);

  mord :=  plex(op(ord));
  res := not `or`(
    seq( Groebner[TestOrder](
      Groebner[TrailingTerm](numer(c), mord)[2],
      Groebner[TrailingTerm](denom(c), mord)[2],
      mord ), c in co ) );

  if res then
    return normal(subs(map(`=`, svars, 1), R)), ord;
  else
    return FAIL;
  end if;
end proc;

# Input :
#   - S, a binomial sum
#   - name, a name
#
# Output :  R, L
#   - R, a rational function
#   - L, a list of variables
#
# The list L gives an order on the variables.
# S is the residue of R with respect to the variables in L that are not
# parameters
sumtores0 := proc(S, name)
  local R, vars, x;
  R := sumtoct(S, name);
  vars := select(has, indets(R), name);
  return isconvergent(normal(R/mul(x, x in vars)));
end proc;

# Input :
#   - R, a rational function
#   - vars, a list of variables
#   - ord, a positive integer
#
# Returns the first `ord' terms of res_{vars[2..-1]}(R), computed in the Laurent series field
# K((vars[-1]))...((vars[1])).
#
# Very useful to check the consistency of an integral representation.
rser :=
  (R, vars, ord) ->
    map2(foldl, residue, map(normal, series(R, vars[1], ord)), seq(v=0,v in select(member, vars, indets(R))[2..-1])):

# Replaces infinity by maxn in the expression S, replaces Sum by add, Binomial
# by binomial, etc and evaluates.
computesum := (S, maxn) -> eval(subs([Sum=add, infinity=maxn, op(op(op(gfdictcomp)))], S));

###### GEOMETRIC REDUCTION OF PERIODS

inorout_base := proc(S, i)
	local minx, maxx, miny, left, right, dom, middle, ret;
	minx := min(map2(op, 1, S));
	maxx := max(map2(op, 1, S));

	if minx=maxx then return {}; end if;
 	if nops(S[1]) <= 1 or i <= 1 then return {-1} end if;

 	miny := min(map2(op, 2, S));
 	middle := map2(subsop, 2=NULL, select(m -> m[2] = miny, S));

 	left := min(map2(op, 1, middle));
  	right := max(map2(op, 1, middle));
	ret := inorout_base(middle, i-1);

	if left > minx then ret := ret union {1}; end if;
	if right < maxx then ret := ret union {-1}; end if;

	return ret;
end proc;

inorout := proc(P, T, ord)
  local S, vars;
  vars := remove(has, ord, T);
  coeffs(collect(P, [T, op(vars)], distributed), [T, op(vars)], 'mon');
  S := map2(map2, degree, {mon}, [T, op(vars)]);

  return inorout_base(S, ListTools[Search](T, ord));
end proc;

# Hermite reduction
# Input :
#   - R, a rational function
#   - v, a name
#   - cert, a boolean (default: false)
#
# Output :
#   A rational function S such that R - S = T' for some rational function T'
#   and such that S has poles of order at most 1 w.r.t. v (including at infinity)
#
#   If cert=true, then it also returns T.
hermitered := proc(R, v :: name, cert := false)
  local a, d, g, dm, ds, dm2, dms, b, c, k;
  a := numer(R);
  d := denom(R);

  g := 0;
  dm := gcd(d, diff(d, v));
  dm := normal(dm);
  ds := normal(d/dm);

  while degree(dm, v) > 0 do
    dm2 := gcd(dm, diff(dm, v));
    dms := normal(dm/dm2);
    gcdex(-normal(ds*diff(dm,v)/dm), dms, a, v, 'b', 'c');
    a := normal(c - diff(b,v)*ds/dms);
    g := g + b/dm;
    dm := dm2;
  end do;

  if cert then
    return normal(a/ds/dm), g;
  else
    return normal(a/ds/dm);
  end if;
end proc;

# Input:
#   - R, rational function
#   - v, symbol
#   - ord, list of symbols containing v
#
# Returns FAIL or a rational function without v which is the sum of the
# residues of R at `small' poles.
#
# [Implementation using hermitered and Rothstein-Tragger resutant
#   for the residue computation]
ratres := proc(R :: ratpoly, v :: name, ord :: list(name))
  local Rn, F, tot, f, ioo, Q, res, n;

  if type(R, `+`) then
    return map(ratres, R, v, ord);
  end if;

  Rn := hermitered(normal(R), v);
  F := map2(op, 1, select(has, factors(denom(Rn))[2], v));

  tot := 0;
  for f in F do
    if type(f, `+`) then
      ioo := inorout(f, v, ord);

      if ioo = {1} then
        Q := normal(Rn*f/diff(f, v));
        res := collect(resultant(numer(Q)-Z_*denom(Q), f, v), Z_);
        n := degree(res, Z_);
        tot := tot - normal(coeff(res, Z_, n-1)/coeff(res, Z_, n));
      elif nops(ioo) > 1 then
        return FAIL;
      end if;
    else  # f = v
      tot := tot + normal(subs(v=0,normal(v*Rn)));
    end if;
  end do;

  return normal(tot);
end proc;

# Idem
#
# [Implementation using maple's residue(R, v=infinity)]
ratres2 := proc(R, v, ord)
  local terms, F, tot, f, ioo, Q, res, n;

  terms := convert(R, parfrac, v);
  if type(terms, `+`) then
    terms := [op(terms)];
  else
    terms := [terms];
  end if;

  tot := 0;
  for F in terms do
    f := denom(F);
    if eval(f, v=0) = 0 then
      tot := tot + residue(F, v=0);
    elif degree(f, v) > 0 then
      ioo := inorout(f, v, ord);
      if ioo = {1} then
        tot := tot - residue(F, v=infinity);
      elif nops(ioo) > 1 then
        return FAIL;
      end if;
    end if;
  end do;

  return normal(tot);
end proc;

geomred := proc(R :: ratpoly, ord :: list, params :: set)
  local S, red, v;
  S := R;
  for v in indets(R) minus params do
    red := ratres(S, v, ord);
    if red <> FAIL then S := red; end if;
  end do;
  return collect(factor(S), params, factor, distributed);
end proc;

geomredall := proc(R :: ratpoly, ord :: list, params :: set)
  local tbl, cur, ivars, st, red, v, nst, all, deg;

  tbl := table([ {} = normal(R) ]);
  cur := { {} };
  ivars := indets(R) minus params;

  while nops(cur) > 0 do
    st := cur[1];
    red := tbl[st];
    if red <> FAIL then
      for v in ivars minus st do
        nst := st union {v};
        if type(tbl[nst], indexed) then
          tbl[nst] := factor(ratres(red, v, ord));
        end if;
        cur := cur union {nst};
      end do;
    end if;
    cur := cur minus {st};
  end do;

  all := remove(`=`, map(op, [entries(tbl)]), FAIL);
  cur := min(map(nops@indets, all));
  all := select(f -> nops(indets(f))=cur, all);
  deg := f -> degree(denom(f), ivars) + min(0, degree(numer(f),ivars)-degree(denom(f),ivars) +nops(ivars)+1);
  cur := min(map(deg, all));
  all := select(f -> deg(f)=cur, all);

  return map(f -> collect(factor(f), params, factor, distributed), all);
end proc;

# Input :
#   - S, a binomial sum
#   - name, a name
#
# Output :  R, L
#   - R, a rational function
#   - L, a list of variables
#
# The list L gives an order on the variables.
# S is the residue of R with respect to the variables in L that are not
# parameters
sumtores := proc(S, name)
  local R, ord, params, flag;

  R, ord := sumtores0(S, name);
  params := indets(R) intersect indets(S);

  flag := true;
  hasoption([_rest], 'geomred' = boolean, 'flag');
  if flag then
    R := geomred(R, ord, params);
    ord := select(has, ord, indets(R));
  end if;

  return op(subs(packvars([ord], name)[1], [R,ord]));
end proc;

end module:

ext: mpl

max_stars_repo_name: lairez/binomsum

max_issues_repo_name: lairez/binomsum

max_forks_repo_name: lairez/binomsum


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

params_a_kappa := 4.8827323:
params_a_mu    := 0.3511128:
$include "gga_x_pbe.mpl"

dldf_a := [1, -0.1637571, -0.1880028, -0.4490609, -0.0082359]:
dldf_csi_HF := 1 - 0.6144129:

dldf_f := (x, u, t) ->
  + dldf_csi_HF*pbe_f(x)*mgga_series_w(dldf_a, 5, t):

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


text:
`is_element/itloc/I` := (n::nonnegint) -> (i) ->
 type(i,nonnegint) and i < n:

`list_elements/itloc/I` := (n::nonnegint) -> [seq(i,i=0..n-1)]:

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

`is_element/itloc/P` := (n::nonnegint) -> (A) ->
 type(A,set(nonnegint)) and (nops(A) = 0 or max(0,op(A)) < n):

`list_elements/itloc/P` := (n::nonnegint) ->
 [op(combinat[powerset]({seq(i,i=0..n-1)}))];

`is_leq/itloc/P` := (n::nonnegint) -> (A,B) -> evalb (A minus B = {});

`bot/itloc/P` := (n::nonnegint) -> {};

`top/itloc/P` := (n::nonnegint) -> {seq(i,i=0..n-1)};

`sup/itloc/P` := (n::nonnegint) -> (A,B) -> A union B;

`inf/itloc/P` := (n::nonnegint) -> (A,B) -> A intersect B;

`angle/itloc/P` := (n::nonnegint) -> (A,B) ->
 (A = {} or B = {} or max(op(A)) <= min(op(B)));

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

`is_element/itloc/Q` := (n::nonnegint) -> proc(U)
 local N,A,b;

 if not type(U,set(set(nonnegint))) then
  return false;
 fi;

 if n = 0 then
  return evalb(U minus {{}} = {});
 fi;

 if max(0,op(map(op,U))) >= n then
  return false;
 fi;

 N := `list_elements/itloc/I`(n);

 for A in U do
  for b in N minus A do
   if not(member(A union {b},U)) then
    return false;
   fi;
  od;
 od;

 return true;
end:

`is_leq/itloc/Q` := (n::nonnegint) -> (U,V) -> evalb (V minus U = {});

`itloc/u` := (n::nonnegint) -> proc (A::set(nonnegint))
 local BB;

 BB := `list_elements/itloc`(n) minus A;

 return map(B -> A union B,BB);
end:

`bot/itloc/Q` := (n::nonnegint) -> `itloc/u`(n)({});

`top/itloc/Q` := (n::nonnegint) -> {};

`omul/itloc/Q` := (n::nonnegint) -> proc(U,V)
 local UV;

 UV := {seq(seq([A,B],B in V),A in U)};

 UV := select(AB -> `angle/itloc/P`(n)(op(AB)),UV);

 return UV;
end:

`mul/itloc/Q` := (n::nonnegint) -> (U,V) ->
 map(AB -> AB[1] union AB[2],`omul/itloc/Q`(U,V));

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

`is_element/itloc/M` := (n::nonnegint) -> (AB) ->
 type(AB,list) and nops(AB) = 2 and
  `is_element/itloc/P`(n)(AB[1]) and
  `is_element/itloc/P`(n)(AB[2]) and
  `itloc/angle`(n)(op(AB));

`list_elements/itloc/M` := proc(n::nonnegint)
 local P,Q,i;

 P := `list_elements/itloc/P`(n);
 Q := table();
 for i from 0 to n do
  Q[i] := combinat[powerset]({seq(j,j=i..n-1)});
 end:
 return [seq(seq([A,B],B in Q[max(0,op(A))]),A in P)];
end:

`is_leq/itloc/M` := (n::nonnegint) -> (AB0,AB1) ->
 evalb (`is_leq/itloc/P`(n)(AB0[1],AB1[1]) and
        `is_leq/itloc/P`(n)(AB0[2],AB1[2]));

`bot/itloc/M` := (n::nonnegint) -> [{},{}];

`inf/itloc/M` := (n::nonnegint) -> (AB0,AB1) ->
 [AB0[1] intersect AB1[1],AB0[2] intersect AB1[2]];

`zeta/itloc`  := (n::nonnegint) -> B -> [{},B];
`xi/itloc`    := (n::nonnegint) -> A -> [A,{}];

`sigma/itloc` := (n::nonnegint) -> AB -> AB[1] union AB[2];

`alpha/itloc` := (n::nonnegint) -> (i) -> (AB) ->
 [select(a -> a < i,AB[1]), select(a -> a >= i,AB[1]) union AB[2]];

`beta/itloc`  := (n::nonnegint) -> (i) -> (AB) ->
 [select(a -> a <= i,AB[1]), select(a -> a >= i,AB[1]) union AB[2]];

`gamma/itloc` := (n::nonnegint) -> (i) -> (AB) ->
 [AB[1] union select(b -> b < i,AB[2]), select(b -> b >= i,AB[2])];

`delta/itloc` := (n::nonnegint) -> (i) -> (AB) ->
 [AB[1] union select(b -> b <= i,AB[2]), select(b -> b >= i,AB[2])];

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

`is_element/itloc/L` := (n::nonnegint) -> proc(U)
 local N,AB,A,B,a,b,A1,B1,i,j;

 if not type(U,set(list(set(nonnegint)))) then
  return false;
 fi;

 N := `list_elements/itloc/I`(n);

 for AB in U do
  if not(`is_element/itloc/M`(n)(AB)) then
   return false;
  fi;
 od;

 for AB in U do
  A,B = op(AB);
  a := max(0,op(A));
  b := min(n-1,op(B));
  A1 := {seq(i,i=0..b)} minus A;
  B1 := {seq(i,i=a..n-1)} minus B;
  for i in A1 do
   if not(member([A union {i},B]),U) then return false; fi;
  od;
  for j in B1 do
   if not(member([A,B union {j}]),U) then return false; fi;
  od;
 od;

 return true;
end:

`is_leq/itloc/L` := (n::nonnegint) -> (U,V) -> evalb (V minus U = {})

`itloc/v` := (n::nonnegint) -> proc (AB::[set(nonnegint),set(nonnegint)])
 local a,b,A,A1,A2,A3,B,B1,L;

 A,B = op(AB);
 b := min(n-1,op(B));
 A1 := {seq(i,i=0..b)} minus A;

 L := NULL;
 for A2 in combinat[powerset](A1) do
  A3 := A union A2;
  a := max(0,op(A3));
  B1 := {seq(j,j=a..n-1)} minus B;
  L := L,seq([A3,B union B2],B2 in combinat[powerset](B1));
 od:
 return {L};
end:


text:
type $atyp <[10] f32>
var $fff <i32> = 7789
func $foo ( var %i i32 )i32

func $fact (
 var %n i32, var %m f32) void {
  return () }
 # 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:
`is_element/FFbar` := (N::posint) -> (A::set) -> proc(Q)
 local B,C,P,P1;
 global reason;

 P := {op(`list_elements/big_subsets`(A))};

 if not(is_table_on(P)(Q)) then
  reason := [convert(procname,string),"Q is not a table indexed by big subsets of A",eval(Q)];
  return false;
 fi;

 for B in P do
  if not(`is_element/SCP`(N)(B)(Q[B])) then
   reason := [convert(procname,string),"Q[B] is not in SCP(N)(B)",eval(Q[B]),N,B,reason];
   return false;
  fi;

  P1 := {op(`list_elements/big_subsets`(B))} minus {B};

  for C in P1 do
   if not(`is_element/SCP2`(N)(B,C)([Q[B],Q[C]])) then
    reason := [convert(procname,string),"(Q[B],Q[C]) is not in SCP(N)(B,C)",
               eval(Q[B]),eval(Q[C]),N,B,C,reason];
    return false;
   fi;
  od;
 od;

 return true;
end;

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

`is_equal/FFbar` := (N::posint) -> (A::set) -> proc(Q0,Q1)
 local B,P;
 global reason;

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

 for B in P do
  if not(`is_equal/SCP`(N)(B)(Q0[B],Q1[B])) then
   reason := [convert(procname,string),"Q0[B] <> Q1[B]",B,Q0[B],Q1[B],reason];
   return false;
  fi;
 od;

 return true;
end;

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

`is_leq/FFbar` := (N) -> (A) -> proc(Q0,Q1)
 local B,P;
 global reason;

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

 for B in P do
  if not(`is_leq/SCP`(N)(B)(Q0[B],Q1[B])) then
   reason := [convert(procname,string),"Q0[B] is not <= Q1[B]",B,Q0[B],Q1[B],reason];
   return false;
  fi;
 od;

 return true;

end;

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

`build/FFbar` := (N::posint) -> (A::set) -> proc(RTTu)
 local R,TT,i,j,k,n,m,ij,A0,C,TT1,TTS,TTc,T,Q,Q0,D,u,v,M,L,M0,M1,RT,BS,UU,U,a,b;
 R,TT,u := op(RTTu);
 n := nops(A);
 A0 := {seq(i,i=1..n)};
 C := children_map(A0)(TT);
 TT1 := select(T -> nops(T) > 1,TT);
 TTS := sort(map(T -> sort([op(T)]),[op(TT1)]));
 Q := table():
 for T in TTS do
  D := sort([op(C[{op(T)}])],(U,V) -> (U[1] <= V[1]));
  D := map(U -> sort([op(U)]),D);
  m := nops(D);
  v := [seq(u[max(D[i])],i=1..nops(D)-1)];
  M := [seq(i,i=1..m)];
  Q0 := `build/ICP`(N)({op(M)})([M,v]);
  L := [];
  for k from 1 to N do
   M0 := Q0[k];
   M1 := NULL;
   for ij in M0 do
    i,j := op(ij);
    M1 := M1,seq(seq([R[a],R[b]],a in D[i]),b in D[j]);
   od;
   M1 := {M1};
   L := [op(L),M1];
  od;
  RT := {seq(R[a],a in T)};
  Q[RT] := L;
 od:
 TT1 := map(op,{indices(Q)});
 BS := {op(`list_elements/big_subsets`(A))};
 TTc := BS minus TT1;
 for T in TTc do
  UU := select(U -> T minus U = {},TT1);
  m := min(map(nops,UU));
  UU := select(U -> nops(U) = m,UU);
  U := UU[1];
  Q[T] := `res/ICP`(N)(U,T)(Q[U]);
 od:
 return eval(Q);
end:

`unbuild/FFbar` := (N::posint) -> (A::set) -> proc(Q)
 local n,R,r,i,TT,TT0,u,a,b,UU,d,U,P,k;
 n := nops(A);
 R := `totalise/FFbar/ord`(N)(A)(Q);
 r := table():
 for i from 1 to n do r[R[i]] := i; od;
 TT := `critical_tree/FFbar`(N)(A)(Q);
 TT0 := map(U -> map(a -> r[a],U),TT);
 u := NULL;
 for i from 1 to n-1 do
  a := R[i];
  b := R[i+1];
  UU := select(U -> member(a,U) and member(b,U),TT);
  d := min(map(nops,UU));
  UU := select(U -> nops(U) = d,UU);
  U := UU[1];
  P := Q[U];
  k := 1;
  while k < N and member([b,a],P[k]) do k := k+1; od;
  u := u,k;
 od;
 u := [u];
 return [R,TT0,u];
end:

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

`is_interior/FFbar` := (N::posint) -> (A::set) -> proc(Q)
 return `is_separated/SCP`(N)(A)(Q[N]);
end;

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

`inc/ICP/FFbar` := (N::posint) -> (A::set) -> proc(Q0)
 local TT,T,Q;

 Q := table();
 TT := `list_elements/big_subsets`(A);
 for T in TT do
  Q[T] := `res/ACP`(N)(A,T)(Q0);
 od;

 return eval(Q);
end:

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

`res/FFbar/ICP` := (N::posint) -> (A::set) -> proc(Q)
 if not(`is_interior/FFbar`(N)(A)(Q)) then
  return FAIL;
 fi;

 return Q[A];
end;

######################################################################
# The functions below could be made much more efficient, but for
# the moment we will stick with the most obvious approach.

`is_critical/FFbar` := (N::posint) -> (A::set) -> proc(Q,U)
 local UU,a;

 UU := `top/autorel`(U);

 for a in A minus U do
  if UU minus Q[U union {a}][N] <> {} then
   return false;
  fi;
 od;

 return true;
end;

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

`critical_tree/FFbar` := (N::posint) -> (A::set) -> proc(Q)
 select(U -> `is_critical/FFbar`(N)(A)(Q,U),
        {op(`list_elements/nonempty_subsets`(A))});
end;

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

`rank/FFbar` := (N::posint) -> (A::set) -> proc(Q)
 local TT,TT1,T;

 TT := `critical_tree/FFbar`(N)(A)(Q);
 TT1 := select(T -> nops(T) > 1,TT);
 return add(`rank/ACP`(N)(T)(Q[T])-1,T in TT1);
end;

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

`totalise/FFbar/ord` := (N::posint) -> (A::set) -> proc(Q)
 local C,P,k;
 C := proc(a,b)
  option remember;
  if a = b then return true; fi;
  P := Q[{a,b}];
  for k from 1 to N do
   if not(member([b,a],P[k])) then
    return true;
   fi;
   if not(member([a,b],P[k])) then
    return false;
   fi;
  od;
  return true;
 end:
 return sort([op(A)],C);
end:

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

`random_element/FFbar` := (N::posint) -> (A::set) -> proc()
 local n,R,TT,u,i;

 n := nops(A);
 R := `random_element/ord`(A)();
 TT := `random_element/standard_stasheff_trees`(n)();
 u := [seq(rand(1..N)(),i=1..n-1)];
 return `build/FFbar`(N)(A)([R,TT,u]);
end:

`list_elements/FFbar` := (N::posint) -> proc(A::set)
 local n,RR,R,U,u,i,j,TTT,TT;

 n := nops(A);
 RR := `list_elements/ord`(A);
 U := [[]];
 for i from 1 to n-1 do
  U := [seq(seq([op(u),j],j=1..N),u in U)];
 od;
 TTT := `list_elements/standard_stasheff_trees`(n);

 return([
  seq(seq(seq(
   `build/FFbar`(N)(A)([R,TT,u]),
   u in U),TT in TTT),R in RR)
 ]);
end:

`count_elements/FFbar` := (N::posint) -> (A::set) ->
 nops(A)! * N^(nops(A)-1) *
 `count_elements/standard_stasheff_trees`(nops(A));

`list_ordered_elements/FFbar` := (N::posint) -> proc(A::{set,list})
 local n,A0,R,U,u,i,j,TTT,TT;

 n := nops(A);
 A0 := {op(A)};
 R := [op(A)];
 U := [[]];
 for i from 1 to n-1 do
  U := [seq(seq([op(u),j],j=1..N),u in U)];
 od;
 TTT := `list_elements/standard_stasheff_trees`(n);

 return([
  seq(seq(
   `build/FFbar`(N)(A)([R,TT,u]),
   u in U),TT in TTT)
 ]);
end:

`count_ordered_elements/FFbar` := (N::posint) -> (A::set) ->
 N^(nops(A)-1) *
 `count_elements/standard_stasheff_trees`(nops(A));

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

`eta/FFbar` := (N::posint) -> (A::set) -> `if`(nops(A) = 1,table(),FAIL);

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

`gamma/FFbar` := (N::posint) -> (A::set,B::set) -> (p) -> proc(Q,P)
 local R,TT,T,U,M,i;

 R := table();
 TT := `list_elements/big_subsets`(A);

 for T in TT do
  U := map(a -> p[a],T);
  if nops(U) > 1 then
   M := Q[U];
   R[T] := [seq(select(u -> member([p[u[1]],p[u[2]]],M[i]),`top/autorel`(T)),i=1..N)];
  else
   R[T] := eval(P[op(U)][T]);
  fi;
 od;

 return eval(R);
end;

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

`mu/Fbar/FFbar` := (N::posint) -> (A::set) -> proc(x)
 local Q,TT,T;

 Q := table();

 TT := `list_elements/big_subsets`(A);
 for T in TT do
  Q[T] := `mu/W/ACP`(N)(T)(x[T]);
 od;

 return eval(Q);
end;

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

`sigma/Fbar/FFbar` := (N::posint) -> (A::set) -> proc(Q)
 local x,TT,T;

 x := table();

 TT := `list_elements/big_subsets`(A);
 for T in TT do
  x[T] := `sigma/ACP/W`(N)(T)(Q[T]);
 od;

 return eval(x);
end;

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

`describe/FFbar` := (N::posint) -> (A::set) -> proc(Q)
 local TT,TT1,T,s;
 TT := `critical_tree/FFbar`(N)(A)(Q);
 TT1 := select(T -> nops(T) > 1,TT);
 TT1 := sort([op(TT1)],(U,V) -> nops(U) >= nops(V));
 s := "";
 for T in TT1 do
  if s <> "" then s := cat(s,"\n"); fi;
  s := cat(s,`describe/ACP`(N)(T)(Q[T]));
 od:
 return s;
end:


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


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

f0 := s -> 1 + B1*s*log(1 + s) + B2*s*log(1 + log(1 + s)):
f  := x -> f0(X2S*x):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
# A raw simplicial complex is a set of nonempty finite sets, no two
# of which are  contained in each other.  We take this as representing
# a simplicial complex, with the specified sets as the maximal simplices.
#
# A (cooked) simplicial complex is a table with various fields which
# gives a more structured representation of the data.

`is_element/raw_simplicial_complexes` := proc(K)
 global reason;
 local n,i,j;

 if not(type(K,set(set)) or type(K,set(list)) or type(K,list(set)) or type(K,list(list))) then
  reason := ["is_element/raw_simplicial_complexes","K is not a set/list of set/lists",K];
  return false;
 fi;

 if member({},K) or member([],K) then
  reason := ["is_element/raw_simplicial_complexes","K contains the empty set",K];
  return false;
 fi;

 n := nops(K);

 for i from 1 to n-1 do
  for j from i+1 to n do
   if {op(K[i])} minus {op(K[j])} = {} then
    reason := ["is_element/raw_simplicial_complexes","K[i] and K[j] are nested",i,j,K];
    return false;
   fi;
  od;
 od;

 return true;
end:

`vertices/raw_simplicial_complex` := (K) -> map(op,{op(K)});

`is_simplex/raw_simplicial_complex` := proc(K,s)
 local s0,u;

 if not(type(s,{list,set})) then return false; fi;

 s0 := {op(s)};
 if s0 = {} then return false; fi;

 for u in K do
  if s0 minus {op(u)} = {} then
   return true;
  fi;
 od:
 return false;
end:

`is_element/simplicial_complexes` := proc(T)
 if not(type(T,table)) then
  return false;
 fi;

 if not(`is_element/raw_simplicial_complexes`(T["max_simplices"])) then
  return false;
 fi;

 if T["vertices"] <> `vertices/raw_simplicial_complex`(T["max_simplices"]) then
  return false;
 fi;

 return true;
end:

`is_simplex/simplicial_complex` := proc(T,s)
 return `is_simplex/raw_simplicial_complex`(T["max_simplices"],s);
end:

`all_simplices/raw_simplicial_complex` := proc(K)
 local V,K0;
 V := `vertices/raw_simplicial_complex`(K);
 K0 := map(s -> {op(s)},{op(K)});
 K0 := map(op,map(combinat[powerset],K0)) minus {{}};
 K0 := map(s -> sort([op(s)]),K0);
 return K0;
end;

`normalise/raw_simplicial_complex` := proc(K)
 local K0,K1,n,i,j,ok;
 K0 := [op(K)];
 K0 := map(s -> {op(s)},K0);
 K0 := sort(K0,(s,t) -> nops(s) <= nops(t));
 n := nops(K0);
 K1 := NULL;
 for i from 1 to n do
  ok := true;
  for j from i + 1 to n do
   if K0[i] minus K0[j] = {} then
    ok := false;
    break;
   fi;
  od;
  if ok then K1 := K1,K0[i]; fi;
 od;
 K1 := [K1];
 K1 := map(s -> sort([op(s)]),K1);
 return K1;
end:

`cook/raw_simplicial_complex` := proc(K)
 local T;

 T := table();
 T["vertices"] := `vertices/raw_simplicial_complex`(K);
 T["max_simplices"] := K;

 return eval(T);
end;

`edges/raw_simplicial_complex` := proc(K)
 local i,j,s;
 {seq(seq(seq({s[i],s[j]},j=i+1..nops(s)),i=1..nops(s)),s in K)};
end:

`faces/raw_simplicial_complex` := proc(K)
 local i,j,k,s;
 {seq(seq(seq(seq({s[i],s[j],s[k]},k=j+1..nops(s)),j=i+1..nops(s)),i=1..nops(s)),s in K)};
end:

`components/simplicial_complex` := proc(T)
 local P,Q,R,s,v;

 P := {seq({v},v in T["vertices"])};
 for s in T["max_simplices"] do
  Q,R := selectremove(c -> c intersect {op(s)} = {},P);
  P := {map(op,R),op(Q)};
 od:
 T["components"] := P;
 return P;
end:

`reindex/simplicial_complex` := proc(T,f)
 local T0,v;

 T0 := table();
 T0["vertices"] := map(f,T["vertices"]);
 T0["max_simplices"] := map(s -> map(f,s),T["max_simplices"]);

 if member(["embedding"],[indices(T)]) then
  T0["embedding_dim"] = T["embedding_dim"];
  for v in T["vertices"] do
   T0["embedding"][f(v)] := T["embedding"][v];
  od:
 fi;

 return eval(T0);
end:

`set_edges/simplicial_complex` := proc(T,force := false)
 if type(T["edges"],{list,set}) and not force then
  return T["edges"];
 fi;
 T["edges"] := `edges/raw_simplicial_complex`(T["max_simplices"]);
 return T["edges"];
end:

`set_faces/simplicial_complex` := proc(T,force := false)
 if type(T["faces"],{list,set}) and not force then
  return T["faces"];
 fi;
 T["faces"] := `faces/raw_simplicial_complex`(T["max_simplices"]);
 return T["faces"];
end:

`set_all_simplices/simplicial_complex` := proc(T,force := false)
 local K,d,i;

 if (type(T["all_simplices"],{list,set}) and not(force)) then
  return T["all_simplices"];
 fi;

 K := T["max_simplices"];
 T["all_simplices"] := `all_simplices/raw_simplicial_complex`(K);
 d := `dim/raw_simplicial_complex`(K);
 T["dim"] := d;
 T["simplices_by_dim"] := table():
 for i from 0 to d do
  T["simplices_by_dim"][i] :=
   select(s -> nops(s) = i + 1,T["all_simplices"]);
 od:
 return T["all_simplices"];
end:

`dim/raw_simplicial_complex` := (K) -> max(map(nops,K)) - 1;

`dim/simplicial_complex` := proc(T)
 local d;
 d := `dim/raw_simplicial_complex`(T["max_simplices"]);
 T["dim"] := d;
 return d;
end:

`euler_characteristic/raw_simplicial_complex` := proc(K)
 local L;

 L := `all_simplices/raw_simplicial_complex`(K);
 L := [op(L)];
 L := map(s -> (-1)^(nops(s)+1),L);
 return `+`(op(L));
end:

`euler_characteristic/simplicial_complex` := proc(T)
 local L;

 if not(type(T["all_simplices"],{set,list})) then
  T["all_simplices"] := `all_simplices/raw_simplicial_complex`(T["max_simplices"]);
 fi;

 L := T["all_simplices"];
 L := [op(L)];
 L := map(s -> (-1)^(nops(s)+1),L);
 return `+`(op(L));
end:

`has_embedding/simplicial_complex` := proc(T)
 if type(T["embedding"],table) then
  if not(type(T["embedding_dim"],nonnegint))
     and nops(T["vertices"]) > 0 then
   T["embedding_dim"] := nops(T["embedding"][T["vertices"][1]]);
  fi;
  return true;
 else
  return false;
 fi;
end:

`has_normal/simplicial_complex` := proc(T)
 if not(`has_embedding/simplicial_complex`(T)) then
  return false;
 fi;

 return type(T["normal"],table);
end:

# Condensation replaces a raw simplicial complex by an equivalent one in
# which the vertex set has the form {1,...,n} for some n.

`is_condensed/raw_simplicial_complex` := proc(K)
 local V,n,i;
 V := `vertices/raw_simplicial_complex`(K);
 n := nops(V);
 return evalb(V = {seq(i,i=1..n)});
end;

`condense/raw_simplicial_complex` := proc(K)
 local V,ix,n,i,L;
 V := `vertices/raw_simplicial_complex`(K);
 n := nops(V);
 for i from 1 to n do
  ix[V[i]] := i;
 od:
 L := map(s -> map(v -> ix[v],s),K);
 return L;
end:

`is_condensed/simplicial_complex` := proc(T)
 local V,n,i;
 V := {op(T["vertices"])};
 n := nops(V);
 return evalb(V = {seq(i,i=1..n)});
end;

`condense/simplicial_complex` := proc(T)
 local V,T0,n,i;

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

 T0["max_simplices"] := `condense/raw_simplicial_complex`(T["max_simplices"]);
 T0["vertices"] := {seq(i,i=1..n)};

 if `has_embedding/simplicial_complex`(T) then
  T0["embedding_dim"] := T["embedding_dim"];
  T0["embedding"] := table();
  for i from 1 to n do
   T0["embedding"][i] := T["embedding"][V[i]];
  od;

  if `has_normal/simplicial_complex`(T) then
   T0["normal"] := table():
   for i from 1 to n do
    T0["normal"][i] := T["normal"][V[i]];
   od;
  fi;
 fi;

 return eval(T0);
end:

`clone/simplicial_complex` := proc(T)
 local T0,k,keys,i;

 T0 := table():
 keys := [
  "vertices",
  "edges",
  "faces",
  "max_simplices",
  "all_simplices",
  "dim",
  "embedding_dim",
  "embedding",
  "normal"
 ];

 for k in keys do
  if member([k],[indices(T)]) then
   if type(T[k],table) then
    T0[k] := table():
    for i in [indices(T[k])] do
     T0[k][op(i)] := T[k][op(i)];
    od:
   else
    T0[k] := T[k];
   fi;
  fi;
 od:

 return eval(T0);
end:

`restrict/simplicial_complex` := proc(T,V0)
 local V1,T0,K0,P0,P,v;

 V1 := {op(V0)};
 if V1 minus {op(T["vertices"])} <> {} then
  error "V0 is not a subset of the vertex set";
 fi;

 T0 := table();
 T0["vertices"] := V0;
 K0 := map(s -> {op(s)} intersect V1,T["max_simplices"]);
 K0 := `normalise/raw_simplicial_complex`(K0);
 T0["max_simplices"] := K0;

 P := T["embedding"];
 if type(P,table) then
  T0["embedding_dim"] := T["embedding_dim"];
  P0 := table():
  for v in V0 do
   P0[v] := P[v];
  od:
  T0["embedding"] := eval(P0);
 fi;

 return eval(T0);
end:

`disjoint_union/simplicial_complex` := proc(T,U)
 local n,U0,TU,v;

 if `is_condensed/simplicial_complex`(T) and
    `is_condensed/simplicial_complex`(U) then
  n := nops(T["vertices"]);
  U0 := `reindex/simplicial_complex`(U,i -> i + n);
 else
  if {op(T["vertices"])} intersect {op(U["vertex_set"])} <> {} then
   error "Vertex sets are not disjoint";
  fi;
  U0 := eval(U);
 fi;

 TU := table();
 TU["vertices"] := [op(T["vertices"]),op(U0["vertices"])];
 TU["max_simplices"] := [op(T["max_simplices"]),op(U0["max_simplices"])];

 if `has_embedding/simplicial_complex`(T) and
    `has_embedding/simplicial_complex`(U0) and
    T["embedding_dim"] = U0["embedding_dim"] then
  TU["embedding_dim"] := T["embedding_dim"];
  TU["embedding"] := table();
  for v in T["vertices"] do
   TU["embedding"][v] := T["embedding"][v];
  od;
  for v in U0["vertices"] do
   TU["embedding"][v] := U0["embedding"][v];
  od;
 fi;

 return eval(TU);
end:

`join/simplicial_complex` := proc(T,U)
 local n,m,U0,TU,s,t,v;

 if `is_condensed/simplicial_complex`(T) and
    `is_condensed/simplicial_complex`(U) then
  n := nops(T["vertices"]);
  U0 := `reindex/simplicial_complex`(U,i -> i + n);
 else
  if {op(T["vertices"])} intersect {op(U["vertex_set"])} <> {} then
   error "Vertex sets are not disjoint";
  fi;
  U0 := eval(U);
 fi;

 TU := table();
 TU["vertices"] := [op(T["vertices"]),op(U0["vertices"])];
 TU["max_simplices"] := [
  op(T["max_simplices"]),
  op(U0["max_simplices"]),
  seq(seq([op(s),op(t)],t in U0["max_simplices"]),s in T["max_simplices"])
 ];

 if `has_embedding/simplicial_complex`(T) and
    `has_embedding/simplicial_complex`(U0) then
  n := T["embedding_dim"];
  m := U0["embedding_dim"];
  TU["embedding_dim"] := n+m+1;
  TU["embedding"] := table();
  for v in T["vertices"] do
   TU["embedding"][v] := [op(T["embedding"][v]),0$(m+1)];
  od;
  for v in U0["vertices"] do
   TU["embedding"][v] := [0$n,1,op(U0["embedding"][v])];
  od;
 fi;

 return eval(TU);
end:

# The following function checks whether we have an embedding of a raw
# simplicial complex K in R^d.  The argument d is the embedding dimension.
# The argument a should be a list or table, so that a[i] is an element
# of R^d for each index i appearing in K.

`is_embedding/raw_simplicial_complex` := proc(K,d,a)
 global reason;
 local m,Kl,Ka,i,j,k,eqd,bc,rr,AA,BB,u,v,w,p,M;

 m := nops(K);
 Kl := [op(K)];
 Ka := map(s -> [seq(a[i],i in s)],Kl);

 for u in Ka do
  if not `is_element/simplices`(d)(u) then
   reason := ["is_embedding/simplicial_complex","u is not a simplex in R^d",u,d,reason];
   return false;
  fi;
 od;

 eqd := evalb({d+1,op(map(nops,Kl))} = {d+1});

 bc := map(`barycentre/simplex`(d),Ka):
 rr := map(`radius/simplex`(d),Ka):

 if eqd then
  AA := map(`primal_matrix/simplex`(d),Ka):
  BB := map(A -> 1/A,AA):
 fi;

 for i from 1 to m-1 do
  for j from i+1 to m do
   u := Ka[i];
   v := Ka[j];
   p := {op(K[i])} intersect {op(K[j])};
   w := {seq(a[k],k in p)};
   if not(nops(w) = nops(p) and
          w = {op(u)} intersect {op(v)}) then
    return false;
   fi;
   if evalf(`d_2/R`(d)(bc[i],bc[j]) - rr[i] - rr[j]) <= 0 then
    if eqd then
     M := <BB[i].AA[j],BB[j].AA[j]>;
     if min(seq(max(seq(M[i,j],j=1..d+1)),i=1..2*d+2)) >= 0 then
      if not(`intersect_nicely/simplices`(d)(u,v)) then
       reason := ["is_embedding/simplicial_complex","bad intersection",u,v];
       return false;
      fi;
     fi;
    else
     if not(`intersect_nicely/simplices`(d)(u,v)) then
      reason := ["is_embedding/simplicial_complex","bad intersection",u,v];
      return false;
     fi;
    fi;
   fi;
  od;
 od;

 return true;
end:

`star/simplicial_complex` := proc(T,s)
 local K,P,K0,T0,P0,s0,v;

 K := T["max_simplices"];
 K := map(u -> {op(u)},{op(K)});
 if `is_simplex/raw_simplicial_complex`(K,s) then
  s0 := {op(s)};
 elif member(s,T["vertices"]) then
  s0 := {s};
 else
  error "Second argument is not a vertex or simplex";
 fi;

 K0 := select(u -> s0 minus u = {},K);
 T0 := `cook/raw_simplicial_complex`(K0);

 P := T["embedding"];
 if type(P,table) then
  T0["embedding_dim"] := T["embedding_dim"];
  P0 := table():
  for v in T0["vertices"] do
   P0[v] := P[v];
  od:
  T0["embedding"] := eval(P0);
 fi;

 return eval(T0);
end:

`link/simplicial_complex` := proc(T,s)
 local K,s0,T0,T1;

 K := T["max_simplices"];
 K := map(u -> {op(u)},{op(K)});
 if `is_simplex/raw_simplicial_complex`(K,s) then
  s0 := {op(s)};
 elif member(s,T["vertices"]) then
  s0 := {s};
 else
  error "Second argument is not a vertex or simplex";
 fi;

 T0 := `star/simplicial_complex`(T,s);
 T1 := `restrict/simplicial_complex`(T0,{op(T0["vertices"])} minus s0);
 return eval(T1);
end:

# For a connected one-manifold (possibly with boundary), list the
# vertices in natural order.  If the complex is a circle, then
# the first vertex in the list is repeated at the end.

`track/simplicial_complex` := proc(T,v_)
 local V,K,L,E,v,w,t,xx;

 if nops(T["vertices"]) <= 1 then
  return [op(T["vertices"])];
 fi;

 if `dim/simplicial_complex`(T) > 1 then
  error "Dimension is bigger than one";
 fi;

 if nops(`components/simplicial_complex`(T)) > 1 then
  error "Complex is not connected";
 fi;

 V := T["vertices"];
 if type(V,set) then V := sort([op(V)]); fi;

 K := {op(map(s -> {op(s)},T["max_simplices"]))};
 L := table();

 for v in V do
  L[v] := select(u -> u <> v and member({u,v},K),V);
 od;

 if max(seq(nops(L[v]),v in V)) > 2 then
  error "Complex is branched";
 fi;

 E := select(v -> nops(L[v]) = 1, V);
 if E = [] then
  if nargs > 1 then
   if member(v_,V) then
    v := v_;
   else
    error "Invalid starting point";
   fi;
  else
   v := V[1];
  fi;
 else
  if nargs > 1 then
   if member(v_,E) then
    v := v_;
   else
    error "Invalid starting point";
   fi;
  else
   v := E[1];
  fi;
 fi;

 w := L[v][-1];
 t := [v,w];
 xx := {op(L[w])} minus {op(t)};
 while xx <> {} do
  w := xx[1];
  t := [op(t),w];
  xx := {op(L[w])} minus {op(t)};
 od;

 if E = [] then t := [op(t),t[1]]; fi;

 return t;
end:

`barycentric_subdivision/raw_simplicial_complex` := proc(K)
 local V,n,i,j,P,L,d,s,t,p;

 V := `all_simplices/raw_simplicial_complex`(K);

 n := max(map(nops,V));
 for i from 1 to n do
  P[i] := combinat[permute](i);
 od:

 L := NULL;
 for s in K do
  d := nops(s);
  for p in P[d] do
   t := {seq(sort([seq(s[p[i]],i=1..j)]),j=1..d)};
   L := L,t
  od;
 od;

 return({L});
end:

`barycentric_subdivision/simplicial_complex` := proc(T)
 local K,E,K0,V0,T0,E0,d,s,x,v;

 K := T["max_simplices"];
 K0 := `barycentric_subdivision/raw_simplicial_complex`(K);
 V0 := `vertices/raw_simplicial_complex`(K0);
 T0 := table():
 T0["vertices"] := V0;
 T0["max_simplices"] := K0;

 E := T["embedding"];
 if type(E,table) then
  d := T["embedding_dim"];
  T0["embedding_dim"] := d;
  E0 := table();
  for s in V0 do
   x := [0$d];
   for v in s do x := x +~ E[v]; od;
   x := x /~ nops(s);
   E0[s] := x;
  od;
  T0["embedding"] := eval(E0);
 fi;

 return eval(T0);
end:

`triangular_subdivision/simplicial_complex` := proc(T)
 local V,E,F,V0,E0,E1,E2,F0,T0,v,e,f,d,P,P0;

 if `dim/simplicial_complex`(T) > 2 then
  error "Triangular subdivision is only for complexes of dimension <= 2";
  # There is a version for dimension 3, with an extra simplex at the
  # barycentre of each tetrahedron as well as each edge, but we have
  # not coded it yet.  Perhaps there is a version where we adjoin the
  # barycentre of each simplex of odd dimension?  Note that even in
  # dimension two, this triangular subdivision is different from, and
  # more symmetrical than, the edgewise subdivision.  The edgewise
  # subdivision needs some kind of order on the vertices.
 fi;
 V := T["vertices"];
 E := `edges/raw_simplicial_complex`(T["max_simplices"]);
 F := `faces/raw_simplicial_complex`(T["max_simplices"]);
 E := map(e -> sort([op(e)]),E);
 E := sort([op(E)]);
 F := map(f -> sort([op(f)]),F);
 F := sort([op(F)]);

 V0 := [seq([v],v in V),seq(sort([op(e)]),e in E)];
 E0 := [seq(op([[[e[1]],e],[[e[2]],e]]), e in E),
        seq(op([
	 [[f[1],f[2]],[f[1],f[3]]],
	 [[f[1],f[2]],[f[2],f[3]]],
	 [[f[1],f[3]],[f[2],f[3]]]
	]), f in F)
       ];
 F0 := [seq(op([
         [[f[1]],[f[1],f[2]],[f[1],f[3]]],
         [[f[2]],[f[1],f[2]],[f[2],f[3]]],
         [[f[3]],[f[1],f[3]],[f[2],f[3]]],
	 [[f[1],f[2]],[f[1],f[3]],[f[2],f[3]]]
        ]),f in F)];

 T0 := table():
 T0["vertices"] := V0;
 T0["edges"] := E0;
 T0["faces"] := F0;

 E1 := {op(map(f -> op({[f[1],f[2]],[f[1],f[3]],[f[2],f[3]]}),F0))};
 E2 := {op(E0)} minus E1;
 E2 := [op(E2)];
 T0["max_simplices"] := [op(E2),op(F0)];

 P := T["embedding"];
 if type(P,table) then
  d := T["embedding_dim"];
  T0["embedding_dim"] := d;
  P0 := table();
  for v in V do
   P0[[v]] := P[v];
  od:

  for e in E do
   P0[e] := (P[e[1]] +~ P[e[2]]) /~ 2;
  od:
  T0["embedding"] := eval(P0);
 fi;

 return T0;
end:

`partial_triangular_subdivision/simplicial_complex` := proc(T,S)
 local S0,Si,V,E,F,V0,E0,F0,nV,nS,e,e1,e2,e3,f,is_long,fin,
  long_edges,p,p1,p2,p3,q,i,T0;

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

 E := `edges/raw_simplicial_complex`(T["max_simplices"]);
 F := `faces/raw_simplicial_complex`(T["max_simplices"]);
 E := map(e -> sort([op(e)]),E);
 E := sort([op(E)]);
 F := map(f -> sort([op(f)]),F);
 F := sort([op(F)]);

 is_long := table():
 for e in E do
  is_long[e] := false;
 od:
 for e in S0 do
  is_long[e] := true;
 od:

 fin := false;
 long_edges := table():

 while not fin do
  fin := true;

  for f in F do
   e1 := [f[1],f[2]];
   e2 := [f[1],f[3]];
   e3 := [f[2],f[3]];
   long_edges[f] := select(e -> is_long[e],[e1,e2,e3]);
   if nops(long_edges[f]) = 2 then
    fin := false;
    is_long[e1] := true;
    is_long[e2] := true;
    is_long[e3] := true;
   fi;
  od:
 od:

 S0 := select(e -> is_long[e],E);
 V := T["vertices"];
 nV := nops(V):
 nS := nops(S0):

 Si := table():
 for i from 1 to nS do
  Si[S0[i]] := nV + i;
 od:

 F0 := NULL;
 for f in F do
  e1 := [f[1],f[2]];
  e2 := [f[1],f[3]];
  e3 := [f[2],f[3]];
  if nops(long_edges[f]) = 0 then
  if nops(f) < 3 then print(["err0",f]); fi;
   F0 := F0,f;
  elif nops(long_edges[f]) = 1 then
   e := long_edges[f][1];
   p := ({op(f)} minus {op(e)})[1];
   q := Si[long_edges[f][1]];
   if nops({e[1],p,q}) < 3 then print(["err1",f,e[1],p,q]); fi;
   if nops({e[2],p,q}) < 3 then print(["err2",f,e[2],p,q]); fi;
   F0 := F0,[e[1],p,q],[e[2],p,q];
  elif nops(long_edges[f]) = 2 then
   print(["shouldn't happen",f]);
  else
   p1 := Si[e1]; p2 := Si[e2]; p3 := Si[e3];
   if nops({f[1],p1,p2}) < 3 then print(["err3",f,p1,p2,p3]); fi;
   if nops({f[2],p1,p3}) < 3 then print(["err4",f,p1,p2,p3]); fi;
   if nops({f[3],p2,p3}) < 3 then print(["err4",f,p1,p2,p3]); fi;
   if nops({p1,p2,p3})   < 3 then print(["err5",p1,p2,p3]); fi;
   F0 := F0, [f[1],p1,p2],[f[2],p1,p3],[f[3],p2,p3],[p1,p2,p3];
  fi;
 od:
 V0 := [seq(i,i=1..nV+nS)];
 F0 := sort([op(map(sort,{F0}))]);
 T0 := table():
 T0["vertices"] := V0;
 T0["max_simplices"] := F0;
 T0["embedding_dim"] := T["embedding_dim"]:
 T0["embedding"] := table():
 for i from 1 to nV do T0["embedding"][i] := T["embedding"][i]; od;
 for i from 1 to nS do
  e := S0[i];
  T0["embedding"][nV+i] := 0.5 *~ (T["embedding"][e[1]] +~ T["embedding"][e[2]]);
 od;
 `set_faces/simplicial_complex`(T0,true);
 `set_edges/simplicial_complex`(T0,true);

 return eval(T0);
end:

`star_subdivide/simplicial_complex` := proc(T,s)
 local n,V,K,K1,S0,S1,S2,s0,s1,s2,u,P,d,x,v,i;

 if nops(s) = 1 then return; fi;

 n := nops(T["vertices"]);
 V := [seq(i,i=1..n)];
 if {op(T["vertices"])} <> {op(V)} then
  error "Star subdivision is only implemented for condensed complexes";
 fi;

 K := T["max_simplices"];
 K := map(u -> {op(u)},{op(K)});

 s0 := {op(s)};
 S1,S0 := selectremove(u -> s0 minus {op(u)} = {},K);
 S2 := {seq(seq(u minus {v},v in s),u in S1)};

 K1 := {op(S0),seq({op(u),n+1},u in S2)};
 K1 := map(u -> sort([op(u)]),K1);
 K1 := sort([op(K1)]);

 T["vertices"] := [seq(i,i=1..n+1)];
 T["max_simplices"] := K1;

 P := T["embedding"];
 if type(P,table) then
  d := T["embedding_dim"];
  x := [0$d];
  for v in s do x := x +~ P[v]; od;
  x := evalf(x /~ nops(s));
  P[n+1] := x;
 fi;
end:

# This function can be applied to a simplicial surface T with a
# vertex v whose link is an octagon, thought of as a bisected
# square.  The vertex w is the midpoint of one edge of the square.
# This function performs a subdivision that affects the star of v
# by adding four new vertices at the midpoints of the edges
# joining v to the corners of the square.  It returns a new vertex
# w', so that the process can be repeated using v and w'.
`square_refine/simplicial_complex` := proc(T,v,w)
 local n,K,L,W,P,old_faces,new_faces,i;
 if not(`is_condensed/simplicial_complex`(T)) then
  error "Complex is not condensed";
 fi;

 n := nops(T["vertices"]);
 L := `link/simplicial_complex`(T,v);
 W := `track/simplicial_complex`(L,w);

 if nops(W) <> 9 then
  error "Link is not square";
 fi;

 old_faces := {seq({v,W[i],W[i+1]},i=1..8)};
 new_faces := {
  {v,W[1],n+1},{v,n+1,W[3]},{v,W[3],n+2},{v,n+2,W[5]},
  {v,W[5],n+3},{v,n+3,W[7]},{v,W[7],n+4},{v,n+4,W[1]},
  {W[1],W[2],n+1},{W[2],W[3],n+1},{W[3],W[4],n+2},{W[4],W[5],n+2},
  {W[5],W[6],n+3},{W[6],W[7],n+3},{W[7],W[8],n+4},{W[8],W[1],n+4}
 };

 T["vertices"] := [seq(i,i=1..n+4)];
 K := T["max_simplices"];
 K := {op(map(s -> {op(s)},K))};
 K := (K minus old_faces) union new_faces;
 K := map(s -> sort([op(s)]),K);
 K := sort([op(K)]);
 T["max_simplices"] := K;

 if `has_embedding/simplicial_complex`(T) then
  P := T["embedding"];
  P[n+1] := (P[v] +~ P[W[2]]) /~ 2;
  P[n+2] := (P[v] +~ P[W[4]]) /~ 2;
  P[n+3] := (P[v] +~ P[W[6]]) /~ 2;
  P[n+4] := (P[v] +~ P[W[8]]) /~ 2;
 fi;

 return n+1;
end:

# Here P is expected to be a table indexed by points
# i in the 2-simplex for which 2^m * i is integral.
# The values are either NULL or vectors in R^3, and
# the values at the vertices of the 2-simplex must
# be in R^3.  The third parameter r is supposed to
# be a retraction of R^3 onto a subspace Z, and it
# is expected that the non-NULL values of P should
# lie in that subspace.  The procedure fills in the
# values P[i] for all i so they lie in Z and are
# spaced out in a reasonably regular way.

`midpoint_extend/simplicial_complex` := proc(P,m,r)
 local X0,p,q,i,j,k,x,y,u,N;

 X0 := {[1,-1,0],[-1,1,0],[1,0,-1],[-1,0,1],[0,1,-1],[0,-1,1]};

 for p from 0 to m do
  q := 2^p;
  for i from 0 to q do
   for j from 0 to q - i do
    k := q - i - j;
    x := [i,j,k];
    if P[x /~ q] = NULL then
     N := select(y -> min(op(x +~ y)) >= 0 and
                      min(op(x -~ y)) >= 0 and
                      P[(x +~ y) /~ q] <> NULL and
                      P[(x -~ y) /~ q] <> NULL,X0):
     if N = {} then DEBUG(); fi;
     u := [0,0,0]:
     for y in N do u := u +~ P[(x +~ y)/~ q]; od:
     u := evalf(u /~ nops(N));
     u := r(u);
     P[x /~ q] := u;
    fi:
   od:
  od:
 od:
end:

`normalise_embedding/simplicial_complex` := proc(T)
 local d,P,v,x,i;

 d := T["embedding_dim"];
 P := T["embedding"];
 for v in T["vertices"] do
  x := evalf(P[v]);
  if not(type(x,[numeric $d])) then
   print([v,x]);
   error "Bad vertex position";
  fi;
  x := evalf(x /~ sqrt(add(x[i]^2,i=1..d)));
  P[v] := x;
 od:
end:

`map_in_place/simplicial_complex` := proc(f,T)
 local E,v;

 E := eval(T["embedding"]);
 for v in T["vertices"] do E[v] := f(E[v]); od:
end:

`map/simplicial_complex` := proc(f,T)
 local T0;

 T0 := eval(`clone/simplicial_complex`(T));
 `map_in_place/simplicial_complex`(f,T0);
 return eval(T0);
end:

`plot/raw_simplicial_complex` := proc(K,d,a)
 local V,E,F,v,e,f,P;

 V := `vertices/raw_simplicial_complex`(K);
 E := `edges/raw_simplicial_complex`(K);
 F := `faces/raw_simplicial_complex`(K);

 P :=
  seq(point(a[v]),v in V),
  seq(line(a[e[1]],a[e[2]]),e in E),
  seq(polygon([a[f[1]],a[f[2]],a[f[3]]],colour=red),f in F);

 return display(P,scaling=constrained,axes=none):
end:

`plot/simplicial_complex` := proc(T,a_)
 local V,E,F,v,e,f,a;

 `set_edges/simplicial_complex`(T);
 `set_faces/simplicial_complex`(T);
 V := T["vertices"];
 E := T["edges"];
 F := T["faces"];
 a := `if`(nargs > 1, a_, T["embedding"]);
 T["plot"] := display(
  seq(point(a[v]),v in V),
  seq(line(a[e[1]],a[e[2]]),e in E),
  seq(polygon([a[f[1]],a[f[2]],a[f[3]]],colour=red),f in F),
  scaling=constrained,axes=none):
 return T["plot"];
end:

`surface_plot/simplicial_complex` := proc(T,a_)
 local V,E,F,f,a;

 `set_faces/simplicial_complex`(T);
 F := T["faces"];
 a := `if`(nargs > 1, a_, T["embedding"]);
 T["surface_plot"] := display(
  seq(polygon([a[f[1]],a[f[2]],a[f[3]]],style=patchnogrid),f in F),
  scaling=constrained,axes=none):
 return T["surface_plot"];
end:

`vertex_labels/simplicial_complex` := proc(T)
 local v;
 if T["embedding_dim"] = 2 then
  display(seq(textplot([op(T["embedding"][v]),sprintf("%A",v)]),v in T["vertices"]),
          axes=none,scaling=constrained);
 elif T["embedding_dim"] = 3 then
  display(seq(textplot3d([op(T["embedding"][v]),sprintf("%A",v)]),v in T["vertices"]),
          axes=none,scaling=constrained);
 else
  error "No 2 or 3 dimensional embedding";
 fi;
end:

`set_javascript/simplicial_complex` := proc(T,format_)
 local format,format3,V,F,F1,P,N,V_string,N_string,F_string,T_string,f,f1,n1,n2,v,i;

 format := `if`(nargs > 1,format_,"%f");
 format3 := sprintf("[%s,%s,%s]",format,format,format);
 if not(`has_embedding/simplicial_complex`(T)) then
  error "Complex has no embedding data";
 fi;

 if T["embedding_dim"] <> 3 then
  error "Complex is not in R3";
 fi;

 V := T["vertices"];
 P := eval(T["embedding"]);
 `set_faces/simplicial_complex`(T);
 F := T["faces"];
 F := sort([op(map(f -> sort([op(f)]),F))]);

 V_string :=
  cat("[",StringTools[Join]([seq(sprintf(format3,op(evalf(P[v]))),v in V)],","),"]");

 if type(T["normal"],table) then
  N := eval(T["normal"]);
  N_string :=
   cat("[",StringTools[Join]([seq(sprintf(format3,op(evalf(N[v]))),v in V)],","),"]");
  F1 := NULL;
  for f in F do
   n1 := cross_product(P[f[2]] -~ P[f[1]],P[f[3]] -~ P[f[1]]);
   n2 := N[f[1]] +~ N[f[2]] +~ N[f[3]];
   if add(n1[i] * n2[i],i=1..3) > 0 then
    f1 := [f[1]-1,f[2]-1,f[3]-1];
   else
    f1 := [f[1]-1,f[3]-1,f[2]-1];
   fi;
   F1 := F1,f1;
  od:
  F1 := [F1];
  F_string :=
   cat("[",StringTools[Join]([seq(sprintf("%A",f),f in F1)],","),"]");
  T_string :=
   sprintf("{\npositions:\n%s,\nnormals:\n%s,indices:\n%s }",
           V_string,N_string,F_string);
 else
  F_string :=
   cat("[",StringTools[Join]([seq(sprintf("%A",f),f in F)],","),"]");
  T_string :=
   sprintf("{\npositions:\n%s,indices:\n%s }",
           V_string,F_string);
 fi;

 T["javascript"] := T_string;

 return T_string;
end:


text:
#ifndef SERVER_PROTO_H__
#define SERVER_PROTO_H__

namespace server {

%genproto

}  // namespace server

#endif /* server/proto.h */
%# vim: syntax=mplex

max_stars_repo_name: frsyuki/ccf

max_issues_repo_name: frsyuki/ccf

max_forks_repo_name: frsyuki/ccf


text:
 func $foo (
  var %i i32
#var %i1 i32, var %j1 i32, var %k1 i32
  ) i32 {
   return (
     add i32(dread i32 %i,
       constval i32 -998))}

 func $foo1 (
  var %i i32, var %j i32, var %k i32,
  var %i1 i32, var %j1 i32, var %k1 i32
  ) i32 {
   return (
     add i32(dread i32 %i, dread i32 %j))}

 func $foo2 (
  var %i i32, var %j i32, var %k i32
  ) i32 {
   return (
     add i32(dread i32 %i, constval i32 0x111111111))}

 func $foo3 (
  var %i i64, var %j i64, var %k i32
  ) i64 {
   return (
     add i64(dread i64 %i, constval i64 0x111111111))}

 func $foo4 (
  var %i i64, var %j i64, var %k i32
  )i64  {
   return (
     add i64(dread i64 %i, dread i64 %j))}

 func $foo5 (
  var %i i64, var %j i64, var %k i32
  ) i64 {
   return (
     add i64(dread i64 %i, constval i64 0x11111))}

 func $foo6 (
  var %i f64, var %j f64, var %k i32
  ) f64 {
   return (
     add f64(dread f64 %i, constval f64 2.237))}

 func $foo66 (
  var %i f64, var %j f64, var %k i32
  ) f64 {
   return (
     add f64(dread f64 %i, constval f64 2.237))}
 func $foo7 (
  var %i f32, var %j f32, var %k i32
  ) f32 {
   return (
     add f32(dread f32 %i, dread f32 %j))}

 func $foo8 (
  var %i f32, var %j f32, var %k i32
  ) f32 {
   return (
     add f32(dread f32 %i, constval f32 2.237f))}

 func $foo9 (
  var %i i8, var %j u8, var %k i32) i32 {
   return (
     add i32(dread i32 %i, dread u32 %j))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
program ifst;	{sample25}
var ch : char;
begin
    readln(ch);
    if ch = 'a' then writeln('It is ''a'' ');
    else writeln('It is not ''a'' ')
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing


text:
Sudoku := module ()
    description "Sudoku solver";
    global `type/grid`;
    export string2grid, printGrid, solver;
    local newGrid, blockIndices, block, group;
    local insert, isAny, isOpen, isSolved, minimalCandidates;
    local trialAndErrorSolver, nakedSinglesSolver, hiddenSinglesSolver;
    option package;

    #
    #   A Sudoku grid.
    #   field : a field of the grid. Number 0 indicates an empty field.
    #   candiates : set of numbers, candiadates for the field.
    #
    `type/grid` := 'record( field ::Array, candidates ::Array )';

    printGrid := proc ( x ::grid )
        description "prints Sudoku grid to console";
        local rowIndex ::posint, colIndex ::posint;
        for rowIndex from 1 to 9 do
            for colIndex from 1 to 9 do
                if x:-field[rowIndex,colIndex] > 0 then
                    printf( "%1d", x:-field[rowIndex,colIndex] );
                else
                    printf( "." );
                end if;
            end do;
            printf( "\n" );
        end do;
    end proc;

    newGrid := proc () ::grid;
        description "creates an empty Sudoku grid";
        Record( field = Array(1..9,1..9,fill=0,datatype=integer),
                candidates = Array(1..9,1..9,fill={seq(i,i=1..9)},datatype=set(integer) ) );
    end proc;

    #
    #   Indices of the fields in a block.
    #
    blockIndices := Array( 0..2, 0..2, (row,col)->[seq(seq( [i+row*3,j+col*3] ,j=1..3),i=1..3)] );

    block := proc ( row ::posint, col ::posint ) ::list;
        description "gets all indices of a block";
        ASSERT( 0 < row and row < 10 and 0 < col and col < 10 );
        blockIndices[trunc((row-1)/3),trunc((col-1)/3)];
    end proc;

    group := proc ( i ::posint ) ::list;
        option remember;
        description "1..9 the columns, 10..18 the rows, 19..27 the blocks";
        if 1 <= i and i <= 9 then
            [ [j,i] $ j=1..9 ];
        elif 10 <= i and i <= 18 then
            [ [i-9,j] $ j=1..9 ];
        elif 19 <= i and i <= 27 then
            blockIndices[ trunc((i-19)/3), (i-19) mod 3];
        else
            error "index out of range"
        end if;
    end proc;

    insert := proc ( x ::grid, row ::posint, col ::posint, value ::posint )
        description "sets a field in the Sudoku grid, updates the candidates";
        local index ::posint;
        ASSERT( 1 <= value and value <= 9 );
        ASSERT( value in x:-candidates[row,col] and x:-field[row,col] = 0 );
        x:-field[row,col] := value;
        x:-candidates[row,col] := {};
        for index from 1 to 9 do
            x:-candidates[row,index] := x:-candidates[row,index] minus {value};
            x:-candidates[index,col] := x:-candidates[index,col] minus {value};
        end do;
        for index in block(row,col) do
            x:-candidates[op(index)] := x:-candidates[op(index)] minus {value};
        end do;
    end proc;

    string2grid := proc ( str ) ::grid;
        description "converts a single string or a iterable of strings to a Sudoku grid";
        local result ::grid, rowIndex ::posint, colIndex ::posint, values, v;
        if nops(str) = 1 then
            if length(str) <> 9*9 then
                error "invalid number of characters in string"
            end if;
            result := newGrid();
            values := convert( str, bytes ) -~ convert("0",bytes)[1];
            for rowIndex from 1 to 9 do
                for colIndex from 1 to 9 do
                    v := values[rowIndex*9+colIndex-9];
                    if  1 <= v and v <= 9 then
                        insert( result, rowIndex, colIndex, v );
                    end if;
                end do;
            end do;
        elif nops(str) = 9 then
            result := thisproc( cat(str) );
        else
            error "invalid number of strings"
        end if;
        result;
    end proc;

    isAny := proc ( x ::grid, condition )
        description "checks if a element exists which fulfills the condition";
        local rowIndex ::posint, colIndex ::posint;
        for rowIndex from 1 to 9 do
            for colIndex from 1 to 9 do
                if condition( x:-field[rowIndex,colIndex], x:-candidates[rowIndex,colIndex] ) then
                    return true;
                end if;
            end do;
        end do;
        false;
    end proc;

    isSolved := proc ( x ::grid ) ::boolean;
        description "checks if all fields are solved";
        not isAny( x, (field,candidates) -> field < 1 );
    end proc;

    isOpen := proc ( x ::grid ) ::boolean;
        description "checks if open fields are in the grid";
        not isAny( x, (field,candidates) -> field = 0 and candidates = {}  ) and isAny( x, (field,candidates) -> field = 0 );
    end proc;

    minimalCandidates := proc( x ::grid )
        description "found index (row,col) of minimal candidates set";
        local index, count, rowIndex ::posInt, colIndex ::posInt;
        count := 10;
        for rowIndex from 1 to 9 do
            for colIndex from 1 to 9 do
                if x:-field[rowIndex,colIndex] = 0 and nops(x:-candidates[rowIndex,colIndex]) < count then
                    count := nops(x:-candidates[rowIndex,colIndex]);
                    index := rowIndex, colIndex;
                end if;
            end do;
        end do;
        ASSERT( count > 0 );
        ASSERT( index <> 'index' );
        index;
    end proc;

    trialAndErrorSolver := proc ( x ::grid ) ::boolean;
        description "solves by trial and error";
        local index, v ::posint, backup;
        ASSERT( isOpen(x) and not isSolved(x) );
        index := minimalCandidates(x);
        backup := Record( field=copy(x:-field), candidates=copy(x:-candidates) );
        for v in x:-candidates(index) do
            insert( x, index, v );
            if solver( x ) then
                return true;
            else
                x:-field := copy(backup:-field);
                x:-candidates := copy(backup:-candidates);
            end if;
        end do;
        # All candidates failed
        false;
    end proc;

    nakedSinglesSolver := proc ( x ::grid ) ::integer;
        description "inserts the naked singles";
        local count ::integer, rowIndex ::posInt, colIndex ::posInt;
        # Count number of solved fields.
        count := 0;
        for rowIndex from 1 to 9 do
            for colIndex from 1 to 9 do
                if nops(x:-candidates[rowIndex,colIndex]) = 1 then
                    # Found a naked single. Only one candidate.
                    insert( x, rowIndex, colIndex, x:-candidates[rowIndex,colIndex][1] );
                    count := count + 1;
                end if;
            end do;
        end do;
        count;
    end proc;

    hiddenSinglesSolver := proc ( x ::grid ) ::integer;
        description "inserts the hidden singles";
        local count ::integer, groupIndex ::posInt, g ::list, c:: Array, index ::list, candidate ::posInt;
        # Count number of solved fields
        count := 0;
        # Search in all groups for hidden singles
        for groupIndex from 1 to 3*9 do
            g := group(groupIndex);
            c := Array(1..9,fill=0);
            for index in g do
                # Count all candidates.
                for candidate in x:-candidates[op(index)] do
                    c[candidate] := c[candidate] + 1;
                end do;
            end do;
            for candidate from 1 to 9 do
                if c[candidate] = 1 then
                    count := count + 1;
                    # This candidate is a hidden single.
                    for index in g do
                        if candidate in x:-candidates[op(index)] then
                            # insert the hidden single.
                            insert( x, op(index), candidate );
                        end if;
                    end do;
                end if;
            end do;
        end do;
        count;
    end proc;

    solver := proc ( x ::grid ) ::boolean;
        description "solves sudokus, finds one solution";
        local changes ::integer;
        changes := 1;
        # Use the algorithmic solver until there are no open fields or the solver don't work
        while isOpen(x) and changes > 0 do
            changes :=   nakedSinglesSolver(x)
                       + hiddenSinglesSolver(x);
        end do;
        if isOpen(x) then
            # Last chance. Solve by trial and error (backtracking algorithm)
            trialAndErrorSolver(x);
        end if;
        isSolved(x);
    end proc;

end module;

max_stars_repo_name: UlrichBerntien/Uebungen-Maple

max_issues_repo_name: UlrichBerntien/Uebungen-Maple

max_forks_repo_name: UlrichBerntien/Uebungen-Maple


text:
# An RZ_set is a subset of the circle R/Z that can be expressed as
# a finite union of half-open intervals [u,v).  Any such set can
# be expressed uniquely as
# [t[1],t[2]) u [t[3],t[4]) u ... u [t[2*n-1],t[2*n]) with
# 0 <= t[0] < t[1] < ... < t[2*n] <= 1.

`is_element/RZ_set` := proc(T)
 local n,i;

 if not type(T,list(realcons)) then return false; fi;
 n := nops(T)/2;
 if not(type(n,nonnegint)) then return false; fi;
 if n = 0 then return true; fi;
 if T[1] < 0 or T[2*n] > 1 then return false; fi;
 for i from 1 to 2*n-1 do
  if T[i+1] <= T[i] then return false; fi;
 od;
 return true;
end:

`is_equal/RZ_set` := proc(T,U)
 if nops(T) <> nops(U) then return false; fi;
 return evalb(simplify({0,op(T -~ U)}) = {0});
end:

`is_leq/RZ_set` := proc(T,U)
 evalb(`intersect/RZ_set`(T,`complement/RZ_set`(U)) = `empty/RZ_set`);
end:

`list_elements/RZ_set` := NULL:
`count_elements/RZ_set` := NULL:

`random_element/RZ_set` := proc()
 local n,S,T,i;

 n := rand(0..5)();
 S := {};
 while nops(S) < 2*n do
  S := {seq(rand(0..719)(),i=1..2*n)};
 od:
 T := sort([op(S)] /~ 720);
 T := `shift/RZ_set`(`random_element/RZ`())(T);
 return T;
end:

`length/RZ_set` := proc(T) local i; add((-1)^i * T[i],i=1..nops(T)); end:

`bot/RZ_set`   := []:
`empty/RZ_set` := []:
`top/RZ_set`   := [0,1];
`RZ/RZ_set`    := [0,1];

`complement/RZ_set` := proc(T)
 local n,U;
 n := nops(T)/2;
 if n = 0 then return [0,1]; fi;
 U := NULL;
 if T[1] > 0 then U := U,0,T[1]; fi;
 U := U,op(2..2*n-1,T);
 if T[2*n] < 1 then U := U,T[2*n],1; fi;
 U := [U];
 return U;
end:

`add_interval/RZ_set` := proc(T,u,v)
 local n,i,j,k,a,b;
 if nops(T) = 0 then return [u,v]; fi;
 n := nops(T)/2;
 i := 0;
 while i < n and T[2*i+2] < u do i := i+1; od;
 if i = n then return [op(T),u,v]; fi;
 j := i+1;
 while j <= n and T[2*j-1] <= v do j := j+1; od;
 a := min(u,T[2*i+1]);
 b := max(v,`if`(j = 1,0,T[2*j-2]));
 return [seq(T[k],k=1..2*i),a,b,seq(T[k],k=2*j-1..2*n)];
end:

`union/RZ_set` := proc()
 local U,T,n,i,j;
 if nargs = 0 then
  return `bot/RZ_set`;
 else
  U := args[1];
  for i from 2 to nargs do
   T := args[i];
   n := nops(T)/2;
   for j from 1 to n do
    U := `add_interval/RZ_set`(U,T[2*j-1],T[2*j]);
   od;
  od:
  return U;
 fi;
end:

`intersect/RZ_set` := proc()
 local T,U,V,n,m,i,j,a,b;
 if nargs = 0 then
  return `top/RZ_set`;
 elif nargs = 1 then
  return args[1];
 elif nargs > 2 then
  return `intersect/RZ_set`(args[1],`intersect/RZ_set`(args[2..-1]));
 else
  T := args[1];
  U := args[2];
  n := nops(T)/2;
  m := nops(U)/2;
  V := NULL;
  for i from 1 to n do
   for j from 1 to m do
    a := max(T[2*i-1],U[2*j-1]);
    b := min(T[2*i],U[2*j]);
    if a < b then
     V := V,[a,b];
    fi;
   od;
  od;
  V := sort([V],(x,y) -> (x[1] < y[1]));
  V := map(op,V);
  return V;
 fi;
end:

`symdiff/RZ_set` := proc(T,U)
 local A,B;

 A := `intersect/RZ_set`(T,`complement/RZ_set`(U));
 B := `intersect/RZ_set`(U,`complement/RZ_set`(T));
 return(`union/RZ_set`(A,B));
end:

`dist/RZ_set` := proc(T,U) `length/RZ_set`(`symdiff/RZ_set`(T,U)); end:

`shift/RZ_set` := (t) -> proc(T)
 local U,n,i,j;
 if T = [] then return []; fi;

 U := T +~ t;
 U := U -~ floor(U[1]);
 n := nops(U)/2;
 if U[2*n] <= 1 then return U; fi;

 i := 0;
 while i < n and U[2*i+2] <= 1 do i := i+1; od;
 if U[2*i+1] < 1 then
  if U[2*n] - 1 = U[1] then
   return [0,seq(U[j]-1,j=2*i+2..2*n-1),seq(U[j],j=2..2*i+1),1];
  else
   return [0,seq(U[j]-1,j=2*i+2..2*n),seq(U[j],j=1..2*i+1),1];
  fi;
 else
  if U[2*n]-1 = U[1] then
   return [seq(U[j]-1,j=2*i+1..2*n-1),seq(U[j],j=2..2*i)];
  else
   return [seq(U[j]-1,j=2*i+1..2*n),seq(U[j],j=1..2*i)];
  fi;
 fi;
end:

`interval/RZ_set` := proc(a,b)
 local a0,b0;

 a0 := a - floor(a);
 b0 := b - floor(b);
 if a0 < b0 then
  return [a0,b0];
 elif a0 = b0 then
  return [];
 else
  return [0,b0,a0,1];
 fi;
end:

`unroll/RZ_set` := proc(T,k::posint := 2)
 local n,i,j,t;
 n := nops(T)/2;
 if n=0 then return []; fi;
 if k = 1 then return T; fi;

 if T[1] = 0 and T[2*n] = 1 then
  return [seq(T[i],i=1..2*n-1),
          seq(seq(T[i]+j,i=2..2*n-1),j=1..k-2),
          seq(T[i]+k-1,i=2..2*n)];
 else
  return [seq(seq(t+j,t in T),j=0..k-1)];
 fi;
end:

`starts/RZ_set` := proc(T)
 local n,i;
 n := nops(T)/2;

 if n = 0 then return {}; fi;

 if T[1] = 0 and T[2*n] = 1 then
  return {seq(T[2*i-1],i=2..n)};
 else
  return {seq(T[2*i-1],i=1..n)};
 fi;
end:

`boundary/RZ_set` := proc(T)
 local n,i;
 n := nops(T)/2;

 if n = 0 then return {}; fi;

 if T[1] = 0 and T[2*n] = 1 then
  return {seq(T[i],i=2..2*n-1)};
 else
  return {seq(T[i],i=1..2*n)};
 fi;
end:

`are_crossed/RZ_set` := proc(T,U)
 local A,B,n,i;
 if `intersect/RZ_set`(T,U) <> `empty/RZ_set` then
  return true;
 fi;

 if nops(T) = 0 or nops(U) = 0 then
  return false;
 fi;

 A := `shift/RZ_set`(-T[1])(T);
 B := `shift/RZ_set`(-T[1])(U);
 i := 1;
 n := nops(A)/2;
 while i < n and A[2*i+2] <= B[1] do
  i := i+1;
 od;
 if i = n then
  return false;
 else
  return evalb(B[nops(B)] > A[2*i+1]);
 fi;
end:

`centroid_C/RZ_set` := proc(T)
 local n,l,t,i;
 n := nops(T)/2;
 if n = 0 then error("T is empty"); fi;
 l := `length/RZ_set`(T);
 evalf(add(int(exp(2*Pi*I*t),t=T[2*i-1]..T[2*i]),i=1..n)/l);
end:

`centroid_R2/RZ_set` := proc(T)
 local z;
 z := `centroid_C/RZ_set`(T);
 return([Re(z),Im(z)]);
end:

`flat_plot/RZ_set` := proc(T,y)
 local i;

 if T = [] then
  return NULL;
 else
  return
   display(seq(line([T[2*i-1],y],[T[2*i],y],args[3..-1]),i=1..nops(T)/2),
           scaling=constrained,axes=none);
 fi;
end:

`plot/RZ_set` := proc(T,r := 1,c := [0,0])
 local A,i,n,z,t;
 if nargs > 3 then
  A := args[4..-1];
 else
  A := NULL;
 fi;
 if T = [] then
  return NULL;
 else
  n := nops(T)/2;
  z := [seq(c +~ r *~ [cos(2*Pi*T[2*i-1]),sin(2*Pi*T[2*i-1])],i=1..n)];
  if T[1] = 0 and T[2*n] = 1 then
   z := [op(2..-1,z)];
  fi;
  return
   display(
    seq(plot([c[1]+r*cos(2*Pi*t),
	      c[2]+r*sin(2*Pi*t),
	      t = T[2*i-1]..T[2*i]],
	      A),
	i=1..nops(T)/2),
    seq(line(0.95 *~ z[i],z[i],A),i=1..nops(z)),
    scaling=constrained,axes=none);
 fi;
end:


text:
  ReparamDetermined := proc(lo :: LO(name, anything), kb::t_kb)
    local h;
    h := op(1,lo);
    LO(h,
       evalindets(op(2,lo),
                  'And'('specfunc({Int,int})',
                        'anyfunc'(anything, 'name=anything')),
                  g -> `if`(determined(op(1,g),h), reparam(g,h,kb), g)))
  end proc;

  determined := proc(e, h :: name)
    local i;
    for i in indets(e, 'specfunc({Int,int})') do
      if hastype(IntegrationTools:-GetIntegrand(i),
           'applyintegrand'('identical'(h),
             'dependent'(IntegrationTools:-GetVariable(i)))) then
        return false
      end if
    end do;
    return true
  end proc;

  #Beginning of Carl's code devoted to disintegration and the reparametrization (aka change
  #of variables) of integrals and sums.

  #Finds the innermost Ints or Sums in an expression, that is, those which
  #don't contain further Ints or Sums
  innermostIntSum:= proc(e::anything, $)::set(specfunc({Int,Sum}));
    select(IS-> nops(indets(IS, specfunc({Int,Sum}))) = 1, indets(e, specfunc({Int,Sum})))
  end proc;

  #my substitute for IntegrationTools:-Change
  ChangeVarInt:= proc(J::specfunc(Int), target::algebraic, $)
    local
    newJ, #What J will become.
    x::name:= op([2,1], J),
    u::name:= gensym(x),    #new variable of integration
    s,                      #What x will become
    F                       #Boolean: flip integral?
    ;
    s:= {solve({u = target}, {x})};
    if nops(s) = 0 then
      userinfo(1, 'reparam', "Can't solve substitution target.");
      return J
    end if;
    if nops(s) > 1 then
      userinfo(1, 'reparam', "Case of multi-branched substitution not handled.");
      return J
    end if;
    s:= s[];

    newJ:= Int(
      eval(implicitdiff(u = target, x, u)*op(1,J), s),
      u=
      limit(target, x= op([2,2,1], J), right).. #lower limit
      limit(target, x= op([2,2,2], J), left),   #upper limit
      op(3.., J) #optional Int args (unlikely to be used)
      );

    #If lower limit > upper limit, then flip limits of integration.
    F:= is(op([2,2,1], newJ) > op([2,2,2], newJ));
    if F::truefalse then
      if F then
        userinfo(2, reparam, "Switching limits:", op([2,2], newJ));
        newJ:= IntegrationTools:-Flip(newJ)
      end if
    else #If inequality can't be decided, then don't reverse.
      userinfo(1, reparam, "Can't order new limits:", op([2,2], newJ))
    end if;
    newJ
  end proc;

  #main procedure for int/sum reparamterizations
  reparam:= proc(e::algebraic, h::symbol, ctx::{list,t_kb}:= [], $)
    local
    J:= innermostIntSum(e),   #the integral or sum
    newJ, #What J will become possible subs target(s)
    oldarg:=
    map2(op, 2, indets(e, specfunc(applyintegrand))),
    kb := build_kb(ctx,"reparam"),
    newarg, Jbnds, del
    ;
    if not(oldargs::{algebraic, set(algebraic), list(algebraic)}) then
      userinfo(2, 'procname', "Unexpected arguements to applyintegrand");
      return e;
    end if;
    if nops(J) = 0 then
      userinfo(2, 'procname', "No sum or integral found.");
      return e
    end if;
    if nops(J) > 1 then
      userinfo(1, 'procname', "Case of multiple innermost Int or Sums not yet handled.");
      return e
    end if;
    J:= J[];
    if J::specfunc(Sum) then
      userinfo(1, 'procname', "Sums not yet handled.");
      return e
    end if;
    if hastype(op(1,J), t_pw_or_part) then
      userinfo(1, 'procname', "Reparam of piecewise body not yet handled.");
      return e;
    end if;
    if nops(oldarg) <> 1 then
      userinfo(1, 'procname', "More than 1 reparam possible:", oldarg);
      return e
    end if;
    oldarg:= oldarg[];   #Extract the reparam target.

    Jbnds := op([2,2], J);
    del := simplify_assuming(op(2,Jbnds)-op(1,Jbnds), kb);
    if del in {-1,1} then
      userinfo(2, 'procname', "No need for a reparameterization (bound diff is 1).");
      return e;
    end if;
    if has(Jbnds, infinity) then
      userinfo(2, 'procname', "Reparam of infinity bounds probably won't work.");
      return e;
    end if;

    # if the arg to applyintegrand is not a variable or a datum or an index
    # into an array
    if oldarg::And(Not({symbol, function, indexed, list})) then
      newarg := oldarg;
    # if the bounds don't contain infinity, do not have a difference of 1
    # and contain a product of that variable
    elif oldarg::symbol and
         op(1,J)::`*`   and
         select(type, [op(op(1,J))], identical(oldarg))<>[]
    then
      newarg := oldarg/del;
    else
      userinfo(2, 'procname', "No need for a reparameterization.");
      return e
    end if;

    #Make the change of vars.
    newJ:= kb_assuming_mb(ChangeVarInt@op, [J, newarg], kb, _->FAIL);

    if newJ = 0 or not(newJ::specfunc(Int)) or newJ=FAIL then
      userinfo(
        1, 'procname',
        printf(
          "Invalid reparam, likely due to improper handling of an infinity issue.\n"
          "u subs:%a",
          #Reformat the ChangeVarInt command for readability.
          subs(
            x= ':-x',
            'ChangeVarInt'(
              subsindets(
                J,
                   specfunc(applyintegrand),
                f-> ':-h'(op(2,f)))),
            newarg)));
      return e;
    end if;
    subs(J= newJ, e)
  end proc;

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
func $fact (
 var %n i32, var %p1 a64, var %p2 a64) i32 {
  var %r1 i32
  var  %r2 i32
  var  %r3 i32
 if (ne i32 i32 (dread i32 %n, constval i32 1)) {
    # Not correct, but for parse and dump test.
    intrinsiccallassigned JSOP_ADD (
      dread a64 %p1, dread a64 %p2, dread i32 %n) {
      dassign %r1 1
      dassign %r2 2
      dassign %r3 3
    }
    return (constval i32 0)
  }
}

max_stars_repo_name: MapleSystem/OpenArkCompiler

max_issues_repo_name: MapleSystem/OpenArkCompiler

max_forks_repo_name: MapleSystem/OpenArkCompiler


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 *)
(* replace: "BesselK\(0," -> "xc_bessel_K0(" *)
(* replace: "BesselK\(1," -> "xc_bessel_K1(" *)
(* replace: "int1\(" -> "xc_integrate(func1, NULL, 0.0, " *)
(* replace: "int2\(" -> "xc_integrate(func2, NULL, 0.0, " *)
(* prefix:
  lda_x_1d_exponential_params *params;

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

$include "lda_x_1d_exponential.mpl"

x1d_inter := x -> 2.0*BesselK(0, x):


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

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

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

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

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

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

  is_minimal := true;
  for U in TT do
   if U <> V and U minus V = {} then
    is_minimal := false;
    if h[U] > h[V] then return false; fi;
   fi;
  od;
  if is_minimal and h[V] <> 0 then
   reason := [convert(procname,string),"V is minimal in TT but h(V) <> 0",V,TT,h[V]];
   return false;
  fi;
 od;

 return true;
end;

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

`is_equal/tree_height_functions` := (A::set) -> (TT) -> proc(h1,h2)
 local T;
 global reason;

 for T in TT do
  if h1[T] <> h2[T] then
   reason := [convert(procname,string),"h1[T] <>  h2[T]",T,h1[T],h2[T]];
   return false;
  fi;
 od;

 return true;
end;

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

`is_leq/tree_height_functions` := NULL;

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

`random_element/tree_height_functions` := (A::set) -> (TT) -> proc(d::posint := 12)
 local e,TS,T,h,parent,children;

 TS := sort([op(TT)],(U,V) -> nops(U) > nops(V));
 h := table();

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

 for T in TS do
  if children[T] = {} then
   h[T] := 0;
  elif parent[T] = FAIL then
   h[T] := rand(0..d)()/d;
  else
   e := h[parent[T]] * d;
   h[T] := rand(0..e)()/d;
  fi;
 od;

 return eval(h);
end;

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

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

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

`is_singular/tree_height_functions` := (A::set) -> (TT) -> proc(h)
 local TS,n,i,j;

 TS := sort([op(TT)],(U,V) -> nops(U) < nops(V));
 n := nops(TS);
 for i from 1 to n do
  if h[TS[i]] = 1 then return true; fi;
  for j from i+1 to n do
   if TS[i] minus TS[j] = {} and h[TS[i]] = h[TS[j]] then
    return true;
   fi;
  od;
 od;

 return false;
end;

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

`extend/tree_height_functions` := (A::set) -> (TT) -> proc(h)
 local he,TS,P,T,U;

 he := table;
 TS := sort([op(TT)],(U,V) -> nops(U) < nops(V));

 P := `list_elements/nonempty_subsets`(A);
 for T in P do
  he[T] := 1;
  for U in TS do
   if T minus U = {} then
    he[T] := h[U];
    break;
   fi;
  od;
 od;

 return eval(he);
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 *)

mu1     := 0.042:
mu2     := 0.26 - mu1:
nu_MGE4 := -0.195:
k2      := -mu2^2/nu_MGE4:
k1      := 0.804 - k2:

f0 := s -> 1 + k1 + k2
   - k1*(1 - mu1*s^2/k1)/(1 - (mu1*s^2/k1)^5)
   - k2/(1 + mu2*s^2/k2):

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


text:
 func $cvtu64tof32 (
  var %i u64
  ) f32 {
   return (
     cvt f32 u64(dread u64 %i))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
# Berechne kinematische Zwangsbedingungen für aufgestellte Viergelenkkette
# Einleitung

# Die kinematischen Zwangsbedingungen werden als Ersetzungsausdruck für die abhängigen Winkel aufgestellt.
#
# fourbar1turnTE -> Viergelenkkette, Modellierung der Zwangsbedingungen mit Ausdrücken für trigonometrische Elimination. Drehbar aufgestellt mit zusätzlichem Gelenk
# kinematic_constraint -> Kinematische Zwangsbedingungen
#
# TODO: Fertige Ergebnisse von Viergelenkkette laden und nicht hier neu berechnen.
# Quelle
# SA Bejaoui: Bejaoui2018_S749; "Modellierung kinematischer Zwangsbedingungen für hybride serielle Roboter mit planaren Parallelmechanismen"
# SA Shan: Shan2019_S828, "Reduktion der Modellkomplexität von seriell-hybriden Palettierrobotern"
# Autor
# Weipu Shan (Studienarbeit), Abderahman Bejaoui (Vorlage fourbar1), Moritz Schappler (Korrekturen)
# Studienarbeit bei: Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2018-08
# (C) Institut fuer Mechatronische Systeme, Leibniz Universitaet Hannover
#
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
#with(ListTools):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../helper/proc_intersect_circle":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):

# Variable mit Winkeln der Nebenstruktur nur in Abhängigkeit der verallgemeinerten Koordinaten
kintmp_qs := Matrix(RowDimension(kintmp_s),1):
kintmpD_t:=diff~(kintmp_t,t):
kintmpD_s:=Matrix(< etaD_s, xiD_s, rhoD_s>):
# Konstante Winkel bereits hineinschreiben ( Keine für Viergelenkkette)
for i from 1 to RowDimension(kintmp_s) do
  if diff(kintmp_s(i,1), t) = 0 then
    kintmp_qs(i,1) := kintmp_s(i,1):
  end if:
end do:
# Variablen definieren für die Hilfswinkel
#
# Ersetzungsausdrücke definieren.
# Speichere Sinus und Cosinus der Winkel direkt ab, da diese in den Rotationsmatrizen direkt auftreten.
# Spalte 1: Zu suchender Ausdruck (sin oder cos eines Winkels)
# Spalte 2: Einzusetzender Ausdruck.
# Dadurch werden arctan-Ausdrücke in der direkten Kinematik reduziert.
# Ähnliches Vorgehen wie in [1].
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
for i from 1 to RowDimension(kintmp_s) do
  kintmp_subsexp(2*i-1, 1) := sin(kintmp_s(i,1)):
  kintmp_subsexp(2*i,   1) := cos(kintmp_s(i,1)):
  # Initialisierung der rechten Spalte mit gleichen Werten. Später nur Ersetzung, wenn Vorteilhaft.
  kintmp_subsexp(2*i-1, 2) := kintmp_subsexp(2*i-1, 1):
  kintmp_subsexp(2*i,   2) := kintmp_subsexp(2*i,   1):
end do:
printf("Beginn der Berechnungen. %s\n", FormatTime("%Y-%m-%d %H:%M:%S")):
st := time():
# Reihenfolge, in der die Winkel der Parallelstruktur berechnet werden. Hilft beim exportieren von Code, wenn beim Debuggen in der gleichen Reihenfolge ausgegeben wird.
Reihenfolge_kintmp := <1;3;2>:
# Zwangsbedingung mit Kreisschnittpunkt C zur Bestimmung von rho
r_p1:= < l1 , 0 >: # Punkt B
;
r_p2:= < l2*cos(qJ2s) , l2*sin(qJ2s)> : # Punkt D
;
r1:=l4:
r2:=l3:
A:=intersect_circle(r_p1, r_p2, r1, r2):
r_s1x:=simplify(A(1,1),trig):
r_s1y:=simplify(A(2,1),trig):
r_acx:=l2*cos(qJ2s) + l3*(-cos(rho_s)*cos(qJ2s)+sin(rho_s)*sin(qJ2s)): # aus Homogener Trafo
;
r_acy:=l2*sin(qJ2s)- l3*(cos(rho_s)*sin(qJ2s)+sin(rho_s)*cos(qJ2s)):  # aus Homogener Trafo
;
f1:=r_s1x-r_acx=0: # GL1
;
f2:=r_s1y-r_acy=0: # GL2
;
# Zwangsbedingungen für eta , zweiter Weg zum Kreisschnittpunkt C
r_acx2:=l1-l4*cos(eta_s):
r_acy2:=l4*sin(eta_s):
f11:=r_s1x-r_acx2=0: # GL3
;
f12:=r_s1y-r_acy2=0: # GL4
;
# Zwangsbedingungen für xi , zwei Wege zum Punkt D
r_adx1:=l2*cos(qJ2s):
r_ady1:=l2*sin(qJ2s):
r_adx2:=l1-l4*cos(eta_s)+l3*(cos(eta_s)*cos(xi_s)-sin(eta_s)*sin(xi_s)):
r_ady2:=l4*sin(eta_s)-l3*(cos(eta_s)*sin(xi_s)+sin(eta_s)*cos(xi_s)): # hier schon einen Fehler erkannt
;
f21:=r_adx1-r_adx2=0: # GL5
;
f22:=r_ady1-r_ady2=0: # GL6
;
# Bestimme Winkel
sol:=simplify(solve({f1,f2},[cos(rho_s),sin(rho_s)]),trig):
sol:=op(sol):
cos(rho_s):=simplify(rhs(sol[1]),trig):
sin(rho_s):=simplify(rhs(sol[2]),trig):
# eta bestimmen
sol1:=solve({f11,f12},{cos(eta_s),sin(eta_s)}):
assign(sol1):
cos(eta1_s):=simplify(rhs(sol1[1]),trig):
sin(eta1_s):=simplify(rhs(sol1[2]),trig):
# xi bestimmen
sol2:=solve({f21,f22},{cos(xi_s),sin(xi_s)}):
sol2:=op(sol2):
cos(xi_s):=simplify(rhs(sol2[1]),trig):
convert_s_t(cos(xi_s)):
sin(xi_s):=simplify(rhs(sol2[2]),trig):
cos(eta_s):=-cos(eta1_s):
sin(eta_s):=sin(eta1_s):
# kintmp_subsexp
kintmp_subsexp(1,2):=sin(rho_s):
kintmp_subsexp(2,2):=cos(rho_s):
kintmp_subsexp(3,2):=sin(eta_s):
kintmp_subsexp(4,2):=cos(eta_s):
kintmp_subsexp(5,2):=sin(xi_s):
kintmp_subsexp(6,2):=cos(xi_s):
# Kintmp_qs
kintmp_qs(1,1):=arctan(sin(rho_s),cos(rho_s)):
kintmp_qs(2,1):=arctan(sin(eta_s),cos(eta_s)):
kintmp_qs(3,1):=arctan(sin(xi_s),cos(xi_s)):
# Exportiere Code für folgende Skripte
kintmp_qt := convert_s_t(kintmp_qs):
# Speichere Maple-Ausdruck (Eingabe-Format und internes Format)
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name):
printf("Ausdrücke für kintmp_subsexp gespeichert (Maple). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
for i from 1 to RowDimension(kintmp_s) do
  tmp := kintmp_qs(i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d", robot_name, i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d.m", robot_name, i):
end do:
save kin_constraints_exist, kintmp_qs, kintmp_qt,kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert", robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
save kintmp_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_qs_maple_inert", robot_name):
printf("Ausdrücke mit Inert-Arctan exportiert (Matlab). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp ));
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m",robot_name),2);
printf("Fertig. %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


text:
var $arr1 <[10] i32>
var $arr2 <[10] i32>
var $arr3 <[10] i32>

func $ArrayAdd () void {
 var %i i32
 doloop %i (
  constval i32 0,
  lt i32 i32 (dread i32 %i, constval i32 10),
  add i32(dread i32 %i, constval i32 1)){
   iassign <* [10] i32>(
    array 1 a32 <* [10] i32>(addrof a32 $arr1, dread i32 %i),
    add i32(
     iread i32 <* [10] i32>(
      array 1 a32 <* [10] i32>(addrof a32 $arr2, dread i32 %i)),
     iread i32 <* [10] i32>(
      array 1 a32 <* [10] i32>(addrof a32 $arr3, dread i32 %i))
     ))
  }
  return()
}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
CY001954    843    C    T    resistance to the neuraminidase inhibitors
CY001953    793    T    C    resistance to the adamantanes
CY001953    805    G    A    resistance to the adamantanes

max_stars_repo_name: stuckinaboot/iGenomics

max_issues_repo_name: stuckinaboot/iGenomics

max_forks_repo_name: stuckinaboot/iGenomics


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

cases := [
    [
        [PolynomialIdeal([a * x + b * y + c * y^2, d * y], variables={x, y}), PolynomialIdeal([a - a_aux, b - b_aux], variables={a_aux, b_aux})],
        PolynomialIdeal([x, y])
    ],
    [
        [PolynomialIdeal([a * x + b * z + c * y^2, d * y], variables={x, y, z}), PolynomialIdeal([a - a_aux, b - b_aux], variables={a_aux, b_aux})],
        PolynomialIdeal([a * x + b * z, y], variables={x, y, z})
    ],
    [
        [PolynomialIdeal([a - a_aux], variables={a_aux}), PolynomialIdeal([a - a_aux, b - b_aux], variables={a_aux, b_aux})],
        PolynomialIdeal([a - a_aux], variables={a_aux})
    ],
    [
        [PolynomialIdeal([-a - b + a_aux + b_aux, -a * a_aux + a * b + a_aux^2 - a_aux * b], variables={a_aux, b_aux}), PolynomialIdeal([-a + a_aux, -b + b_aux], variables={a_aux, b_aux})],
        PolynomialIdeal([-a - b + a_aux + b_aux, -a * a_aux + a * b + a_aux^2 - a_aux * b], variables={a_aux, b_aux})
    ],
    [
        [
            PolynomialIdeal([-a - b + a_aux + b_aux, -a * a_aux + a * b + a_aux^2 - a_aux * b], variables={a_aux, b_aux}),
            PolynomialIdeal([-a - b + a_aux + b_aux, -a * a_aux + a * b + a_aux^2 - a_aux * b, c - c_aux], variables={a_aux, b_aux, c_aux})
        ],
        PolynomialIdeal([-a - b + a_aux + b_aux, -a * a_aux + a * b + a_aux^2 - a_aux * b], variables={a_aux, b_aux})
    ]
];

num_passed := 0:
num_failed := 0:

for case in cases do
    input := case[1]:
    correct := case[2]:
    if IdealsEq(FieldCoefficientRestriction(op(input)), correct) then
        printf("PASSED\n");
        num_passed := num_passed + 1:
    else
        printf("FAILED\n");
        num_failed := num_failed + 1:
        print("Expected: ", correct);
        print("Got: ", FieldCoefficientRestriction(op(input)));
    end if:
end do:

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

max_stars_repo_name: iliailmer/AllIdentifiableFunctions

max_issues_repo_name: iliailmer/AllIdentifiableFunctions

max_forks_repo_name: iliailmer/AllIdentifiableFunctions


text:
(*
 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 *)

# Note that the files have to be included in this specific order
$define gga_x_pbe_params
$include "gga_x_pbe.mpl"

$include "op.mpl"

op_qab         := 2.3789:
op_enhancement := xs -> pbe_f(xs):