text:
# Step 3 of 3: from Maple LO (linear operator) back to Hakaru
fromLO := module()
  export ModuleApply :=
  proc(lo :: LO(name, anything), {_ctx :: t_kb := empty}, $)
    local h;
    h := gensym(op(1,lo));
    _Env_HakaruSolve := false;
    unintegrate(h, eval(op(2,lo), op(1,lo) = h), _ctx)
  end proc;

  export
  unintegrate := proc(h :: name, e, kb :: t_kb_mb, $)
    local x, c, lo, hi, make, m, mm, w, w0, w1, recognition, subintegral,
          i, kb1, kb2, loops, subst, hh, pp, t, bnds, br;
    if kb :: t_not_a_kb then return Msum(); end if;
    if e :: 'And'('specfunc({Int,int})',
                  'anyfunc'('anything','name'='range'('freeof'(h)))) then
      (lo, hi) := op(op([2,2],e));
      x, kb1 := genLebesgue(op([2,1],e), lo, hi, kb);
      subintegral := eval(op(1,e), op([2,1],e) = x);
      (w, m) := unweight(unintegrate(h, subintegral, kb1));
      recognition := recognize_continuous(w, x, lo, hi, kb1);
      if recognition :: 'Recognized(anything, anything)' then
        (w, w0) := factorize(op(2,recognition), x, kb1);
        weight(w0, bind(op(1,recognition), x, weight(w, m)))
      else error "recognize_continuous is never supposed to fail" end if
    elif e :: 'And'('specfunc({Sum,sum})',
                    'anyfunc'('anything','name'='range'('freeof'(h)))) then
      (lo, hi) := op(op([2,2],e));
      x, kb1 := genType(op([2,1],e), HInt(closed_bounds(lo..hi)), kb);
      subintegral := eval(op(1,e), op([2,1],e) = x);
      (w, m) := unweight(unintegrate(h, subintegral, kb1));
      recognition := recognize_discrete(w, x, lo, hi, kb1);
      if recognition :: 'Recognized(anything, anything)' then
        (w, w0) := factorize(op(2,recognition), x, kb1);
        weight(w0, bind(op(1,recognition), x, weight(w, m)))
      else error "recognize_discrete is never supposed to fail" end if
    elif e :: 'And'('specfunc({Ints,ints,Sums,sums})',
                    'anyfunc'('anything', 'name', 'range'('freeof'(h)),
                              'list(name=range)')) then
      loops := op(4,e);
      bnds  := op(3,e);
      if op(0,e) in {Ints,ints} then
        t := HReal(open_bounds(bnds));
        make := Int;
      else
        t := HInt(closed_bounds(bnds));
        make := Sum;
      end if;
      x, kb1 := genType(op(2,e), mk_HArray(t, loops), kb);
      if nops(op(4,e)) > 0 then
        kb1 := assert(size(x)=op([4,-1,2,2],e)-op([4,-1,2,1],e)+1, kb1);
        ASSERT(type(kb1,t_kb), "unintegrate{Ints,Sums}: integral bounds invalid");
      end if;
      subintegral := eval(op(1,e), op(2,e) = x);
      (w, m) := unweight(unintegrate(h, subintegral, kb1));
      w := simplify_factor_assuming(peel(w), kb1); # for "Don't be confused by extra iterations" tests
      bnds, loops, kb2 := genLoop(bnds, loops, kb, 'Integrand'(x,[w,m]));
      w, pp := unproducts(w, x, loops, kb2);
      hh := gensym('ph');
      subintegral := make(pp * applyintegrand(hh,x), x=bnds);
      (w1, mm) := unweight(unintegrate(hh, subintegral, kb2));
      mm := foldl(((mmm,loop) ->
                   Plate(op([2,2],loop) - op([2,1],loop) + 1,
                         op(1,loop),
                         eval(mmm, op(1,loop) = op(1,loop) - op([2,1],loop)))),
                  mm, op(loops));
      w := w * foldl(product, w1, op(loops));
      w := simplify_factor_assuming(w, kb1);
      (w, w0) := factorize(w, x, kb1);
      weight(simplify_factor_assuming(w0, kb),
             bind(mm, x, weight(simplify_factor_assuming(w, kb1), m)))
    elif e :: 'applyintegrand'('identical'(h), 'freeof'(h)) then
      Ret(op(2,e))
    elif e = 0 then
      Msum()
    elif e :: `+` then
      map2(unintegrate, h, Msum(op(e)), kb)
    elif e :: `*` then
      (subintegral, w) := selectremove(depends, e, h);
      if subintegral :: `*` then error "Nonlinear integral %1", e end if;
      (w0, w) := op(Domain:-Extract:-Shape(w));
      w0 := Domain:-Shape:-toConstraints(w0);
      kb1 := foldr(assert, kb, op(w0));
      if kb1 :: t_kb then
        m := weight(w, unintegrate(h, subintegral, kb1));
        if m :: Weight(anything, anything) then
          m := weight(simplify_factor_assuming(op(1,m), kb1), op(2,m));
        end if;
        piecewise_And(w0, m, Msum())
      else # if the domain is empty
        Msum()
      end if;
    elif e :: t_pw and not Partition:-ConditionsDepend(Partition:-PWToPartition(e), h) then
        m := kb_piecewise(e, kb, ((lhs, kb)-> lhs), ((rhs, kb)-> unintegrate(h, rhs, kb)), 'no_split_disj');
        if m :: t_pw and nops(m) = 2 then
          piecewise(op(m), Msum());
        else
          m;
        end if;
    elif e :: Partition and not Partition:-ConditionsDepend(e, h) then
        kb_Partition(e, kb, ((lhs, kb)-> lhs), ((rhs, kb)-> unintegrate(h, rhs, kb)));
    elif e :: t_case then
      subsop(2=map(proc(b :: Branch(anything, anything))
                     eval(subsop(2='toLO:-unintegrate'(x,op(2,b),c),b),
                          {x=h, c=kb})
                   end proc,
                   op(2,e)),
             e);
    elif e :: 'Context(anything, anything)' then
      kb1 := assert(op(1,e),kb);
      if kb1 :: t_kb then
          subsop(2=unintegrate(h, op(2,e), kb1), e);
      else# A contradictory `Context' implies anything, so produce 'anything'
          # In particular, 42 :: t_Hakaru = false, so a term under a false
          # assumption should never be inspected in any way.
          42
      end if;

    elif e :: 'toLO:-integrate'('freeof'(h), 'anything', identical([])) then
      x := mk_sym('x', op(2,e));
      # If we had HType information for op(1,e),
      # then we could use it to tell kb about x.
      (w, m) := unweight(unintegrate(h, applyintegrand(op(2,e), x), kb));
      (w, w0) := factorize(w, x, kb);
      weight(w0, bind(op(1,e), x, weight(w, m)))
    elif e :: identical('undefined') then
      undefined
    else
      # Failure: return residual LO
      LO(h, e)
    end if
  end proc;

  export
  recognize_continuous := proc(weight0, x, lo, hi, kb, $)
    local Constant, de, Dx, f, w, res, rng;
    res := FAIL;
    # gfun[holexprtodiffeq] contains a test for {radfun,algfun} that seems like
    # it should test for {radfun(anything,x),algfun(anything,x)} instead.
    # Consequently, it issues the error "expression is not holonomic: %1" for
    # actually holonomic expressions such as exp(x*sum(g(i,j),j=1..n)).
    # Moreover, mysolve has trouble solve-ing constraints involving sum, etc.
    # To work around these weaknesses, we wrap sum(...), etc. in Constant[...].
    # Unlike sum(...), Constant[sum(...)] passes the type test {radfun,algfun},
    # which we need to handle exp(x*sum(...)) using gfun[holexprtodiffeq].
    # Like sum(...i...), Constant[sum(...i...)] depends on i, which we need so
    # that product(sum(...i...),i=1..m) doesn't simplify to ...^m.
    w := subsindets[flat](weight0,
           And(function, Not(specfunc({exp, And, Or, Not})),
               'freeof'(x)),
           proc(e) Constant[e] end);
    w := subsindets[flat](w, {`^`, specfunc(exp)},
           proc(e)
             applyop(proc(e)
                       subsindets[flat](e,
                         And({`^`, specfunc(exp)},
                             Not(radfun), Not(algfun), 'freeof'(x)),
                         proc(e) Constant[e] end)
                     end,
                     -1, e)
             end);
    de := get_de(w, x, Dx, f);
    if de :: 'Diffop(anything, anything)' then
      res := recognize_de(op(de), Dx, f, x, lo, hi, kb)
    end if;
    if res = FAIL then
      res := Recognized(Lebesgue(lo, hi), w);
      rng := hi - lo;
      if not (rng :: 'SymbolicInfinity') then
        w := simplify_factor_assuming(w * rng, kb);
        # w could be piecewise and simplify will hide the problem
        if not (w :: {'SymbolicInfinity', 'undefined'}) then
          res := Recognized(Uniform(lo, hi), w)
        end if
      end if
    end if;
    # Undo Constant[...] wrapping
    res := subsindets[flat](res, 'specindex'(anything, Constant), x -> op(1,x));
    res
  end proc;

  export
  recognize_discrete := proc(w, k, lo, hi, kb, $)
    local se, Sk, f, a0, a1, lambda, r, s, res;
    res := FAIL;
    if lo = 0 and hi = infinity then
      se := get_se(w, k, Sk, f);
      if se :: 'Shiftop(anything, anything, identical(ogf))' and
         ispoly(op(1,se), 'linear', Sk, 'a0', 'a1') then
        lambda := normal(-a0/a1*(k+1));
        if not depends(lambda, k) then
          res := Recognized(PoissonD(lambda), eval(w,k=0)/exp(-lambda));
        end if;
        if ispoly(lambda, 'linear', k, 'b0', 'b1') then
          r := b0/b1;
          res := Recognized(NegativeBinomial(r, b1), eval(w,k=0)/(1-b1)^r);
        end if
      end if;
    elif lo = 0 and not(hi :: 'SymbolicInfinity') then
      s, r := factorize(simplify_factor_assuming(w, kb), k, kb);
      if s <> 1 then
        s := simplify_factor_assuming(s, kb);
        res := ary(hi+1, k, s);
        if res :: 'list' and nops(convert(res,'set')) = 1 then
          res := Recognized(Counting(lo, hi+1), res[1]);
        else
          res := Recognized(Categorical(res), r);
        end if;
      end if;
    end if;
    if res = FAIL then
      res := Recognized(Counting(lo, hi+1), w);
    end if;
    applyop(simplify_assuming, 1,
            applyop(simplify_factor_assuming, 2, res, kb), kb)
  end proc;

  local
  get_de := proc(dens, var, Dx, f, $)
    :: Or(Diffop(anything, set(function=anything)), identical(FAIL));
    local de, init;
    try
      de := gfun[holexprtodiffeq](dens, f(var));
      de := gfun[diffeqtohomdiffeq](de, f(var));
      if not (de :: set) then
        de := {de}
      end if;
      init, de := selectremove(type, de, `=`);
      if nops(de) = 1 then
        if nops(init) = 0 then
          # TODO: Replace {0, 1/2, 1} by PyMC's distribution-specific "testval"
          init := map(proc (val)
                        try f(val) = eval(dens, var=val)
                        catch: NULL
                        end try
                      end proc,
                      {0, 1/2, 1})
        end if;
        return Diffop(DEtools[de2diffop](de[1], f(var), [Dx, var]), init)
      end if
    catch: # do nothing
    end try;
    FAIL
  end proc;

  local
  get_se := proc(dens, var, Sk, u, $)
    :: Or(Shiftop(anything, set(function=anything), name), identical(FAIL));
    local x, de, re, gftype, init, f;
    try
      # ser := series(sum(dens * x^var, var=0..infinity), x);
      # re := gfun[seriestorec](ser, f(var));
      # re, gftype := op(re);
      _EnvFormal := true;
      de := gfun[holexprtodiffeq](sum(dens*x^var, var=0..infinity), f(x));
      re := gfun[diffeqtorec](de, f(x), u(var));
      re := gfun[rectohomrec](re, u(var));
      if not (re :: set) then
        re := {re}
      end if;
      init, re := selectremove(type, re, `=`);
      if nops(re) = 1 then
        if nops(init) = 0 then
          init := {u(0) = eval(rens, var=0)};
        end if;
        re := map(proc(t)
                    local s, r;
                    s, r := selectremove(type, convert(t, 'list', `*`),
                                         u(polynom(nonnegint, var)));
                    if nops(s) <> 1 then
                      error "rectohomrec result nonhomogeneous";
                    end if;
                    s := op([1,1],s) - var;
                    if s :: nonnegint and r :: list(polynom(anything, var)) then
                      `*`(op(r), Sk^s);
                    else
                      error "unexpected result from rectohomrec"
                    end if
                  end proc,
                  convert(re[1], 'list', `+`));
        return Shiftop(`+`(op(re)), init, 'ogf')
      end if
    catch: # do nothing
    end try;
    FAIL
  end proc;

  local
  recognize_de := proc(diffop, init, Dx, f, var, lo, hi, kb, $)
    local dist, ii, constraints, w, a0, a1, a, b0, b1, c0, c1, c2, loc, nu;
    dist := FAIL;
    if lo = -infinity and hi = infinity
       and ispoly(diffop, 'linear', Dx, 'a0', 'a1') then
      a := normal(a0/a1);
      if ispoly(a, 'linear', var, 'b0', 'b1') then
        dist := Gaussian(-b0/b1, sqrt(1/b1))
      elif ispoly(numer(a), 'linear', var, 'b0', 'b1') and
           ispoly(denom(a), 'quadratic', var, 'c0', 'c1', 'c2') then
        loc := -c1/c2/2;
        if Testzero(b0 + loc * b1) then
          nu := b1/c2 - 1;
          if Testzero(nu - 1) then
            dist := Cauchy(loc, sqrt(c0/c2-loc^2))
          else
            dist := StudentT(nu, loc, sqrt((c0/c2-loc^2)/nu))
          end if
        end if
      end if;
    elif lo = 0 and hi = 1
         and ispoly(diffop, 'linear', Dx, 'a0', 'a1')
         and ispoly(normal(a0*var*(1-var)/a1), 'linear', var, 'b0', 'b1') then
      dist := BetaD(1-b0, 1+b0+b1)
    # elif not evalb((hi - lo) :: 'SymbolicInfinity')
    #      and ispoly(diffop, 'linear', Dx, 'a0', 'a1')
    #      and ispoly(a0 - 2*var, 'linear', var, 'b0', 'b1') then
    #   c0 := (lo*b1 + hi + lo + b0) / (hi - lo);
    #   c1 := -(hi*b1 + hi + lo + b0) / (hi - lo);
    #   if c0 = 1 and c1 = 1 then
    #       dist := Uniform(lo, hi)
    #   else
    #       dist := bind(BetaD(c0, c1),x,lo+(hi-lo)*x)
    #   end if
    elif lo = 0 and hi = infinity
         and ispoly(diffop, 'linear', Dx, 'a0', 'a1')
         and ispoly(normal(a0*var/a1), 'linear', var, 'b0', 'b1') then
    #  if Testzero(b1-1/2) then
    #    dist := ChiSquared(2*(1-b0))
    #  else
        dist := GammaD(1-b0, 1/b1)
    #  end if;
    end if;
    if dist <> FAIL then
      try
        ii := map(convert, init, 'diff');
        constraints := eval(ii, f = (x -> w*density[op(0,dist)](op(dist))(x)));
        w := eval(w, mysolve(constraints, w));
        if not (has(w, 'w')) then
          return Recognized(simplify_assuming(dist, kb),
                            simplify_factor_assuming(w, kb));
        end if
      catch: # do nothing
      end try;
      WARNING("recognized %1 as %2 but could not solve %3", f, dist, init)
    end if;
    FAIL
  end proc;

  # (s,r):=factorize(e,var,kb) expresses e in the context kb as s*r,
  # where r doesn't depend on var and s is as simple as possible
  # (and non-negative if possible).
  local
  factorize := proc(e, var, kb, $)
    local res, x, y, kb1, s, r;
    if not depends(e, var) then
      return 1, e;
    end if;
    if e :: `*` then
      res := map(`[]`@factorize, list_of_mul(e), var, kb);
      return `*`(op(map2(op,1,res))),
             `*`(op(map2(op,2,res)));
    end if;
    if e :: 'anything^freeof(var)' then
      s, r := factorize(op(1,e), var, kb);
      return s^op(2,e),
             r^op(2,e);
    end if;
    if e :: 'And(specfunc({product,Product}),
                 anyfunc(anything, name=range(freeof(var))))' then
      x, kb1 := genType(op([2,1],e), HInt(closed_bounds(op([2,2],e))), kb, var);
      s, r := factorize(eval(op(1,e), op([2,1],e)=x), var, kb1);
      return op(0,e)(s, x=op([2,2],e)),
             op(0,e)(r, x=op([2,2],e));
    end if;
    if e :: 'And(specfunc({product,Product}),
                 anyfunc(anything, name=range))'
       and not depends(subsop([2,2,2]=undefined,e), var) then
      s, r := termize(op([2,2,2],e), var, kb);
      x := op([2,1],e);
      y := `if`(depends(r,x), gensym(x), x);
      return op(0,e)(eval(op(1,e),x=r+1+y), y=0..s-1),
             op(0,e)(op(1,e), x=op([2,2,1],e)..r);
    end if;
    e, 1;
  end proc;

  # (s,r):=termize(e,var,kb) expresses e in the context kb as s+r,
  # where r doesn't depend on var and s is as simple as possible.
  local
  termize := proc(e, var, kb, $)
    local res, x, y, kb1, s, r, i, conds, pw;
    if not depends(e, var) then
      return 0, e;
    end if;
    if e :: `+` then
      res := map(`[]`@termize, [op(e)], var, kb);
      return `+`(op(map2(op,1,res))),
             `+`(op(map2(op,2,res)));
    end if;
    if e :: `*` then
      s, r := selectremove(depends, e, var);
      if r <> 1 then return op(map(`*`, [termize(s, var, kb)], r)) end if;
    end if;
    if e :: 'And(specfunc({sum,Sum}),
                 anyfunc(anything, name=range(freeof(var))))' then
      x, kb1 := genType(op([2,1],e), HInt(closed_bounds(op([2,2],e))), kb, var);
      s, r := termize(eval(op(1,e), op([2,1],e)=x), var, kb1);
      return op(0,e)(s, x=op([2,2],e)),
             op(0,e)(r, x=op([2,2],e));
    end if;
    if e :: 'And(specfunc({sum,Sum}),
                 anyfunc(anything, name=range))'
       and not depends(subsop([2,2,2]=undefined,e), var) then
      s, r := termize(op([2,2,2],e), var, kb);
      x := op([2,1],e);
      y := `if`(depends(r,x), gensym(x), x);
      return op(0,e)(eval(op(1,e),x=r+1+y), y=0..s-1),
             op(0,e)(op(1,e), x=op([2,2,1],e)..r);
    end if;
    if e :: 'specfunc(piecewise)' then
      conds := [seq(op(i,e), i=1..nops(e)-1, 2)];
      if depends(conds, var) then
        # Too bad the conditions depend on var.
        # But maybe the conditions depend on var only in certain cases
        # (whose conditions in turn do not depend on var)?
        pw := select(proc(pw, $)
                        local i;
                        if not depends(pw, var) then return false end if;
                        for i from 1 by 2 to nops(pw)-1 do
                          if depends(op(i,pw), var) then return false end if;
                        end do;
                        for i from 2 by 2 to nops(pw)-1 do
                          if not depends(op(i,pw), var) then return true end if;
                        end do;
                        return not depends(op(-1,pw), var);
                      end proc,
                      indets(conds, 'specfunc(piecewise)'));
        if nops(pw) > 0 then
          pw := op(1, pw); # Pick any eligible piecewise to lift
          pw := piecewise(seq(`if`(i::odd and i<nops(pw),
                                   op(i,pw),
                                   subs(pw=op(i,pw), e)),
                              i=1..nops(pw)));
          return termize(pw, var, kb);
        end if;
      else
        # Yay, the conditions don't depend on var.
        # So just map into the piecewise.
        res := kb_piecewise(e, kb, ((cond, kb) -> cond),
                                   ((ee, kb) -> [termize(ee, var, kb)]));
        if res :: 'specfunc(piecewise)'
           and [seq(op(i,res), i=2..nops(res)-1, 2), op(-1,res)]
               :: 'list([anything, anything])' then
          return piecewise(seq(op(`if`(i::odd and i<nops(res), i, [i,1]), res),
                               i=1..nops(res))),
                 piecewise(seq(op(`if`(i::odd and i<nops(res), i, [i,2]), res),
                               i=1..nops(res)));
        elif res :: '[anything, anything]' then
          return op(res)
        end if;
      end if;
    end if;
    e, 0;
  end proc;

  export
  unweight := proc(m, $)
    local total, ww, mm;
    if m :: 'Weight(anything, anything)' then
      op(m)
    elif m :: 'specfunc(Msum)' then
      total := `+`(op(map((mi -> unweight(mi)[1]), m)));
      (total, map((mi -> weight(1/total, mi)), m))
    else
      (1, m)
    end if;
  end proc;

end module; # fromLO

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
   {  InformerScen.mpl  }

TITLE
   InformerScen;

STOCHASTIC

TREE
   TwoStage;

SCENARIOS
   scen := 1..10;

RANDOM DATA
   Demand[scen] := (210, 290, 335, 370, 415, 450, 460, 495, 540, 585);

PROBABILITIES
   prob[scen] := (0.03,  0.134, 0.086, 0.188, 0.226, 0.08, 0.056, 0.132, 0.06,  0.008);

DATA
   Budget    := 1500;
   ProdCost  := 0.81;
   SellPrice := 0.90;
   Margin    := SellPrice - ProdCost;

DECISION VARIABLES
   Production;

STAGE2 VARIABLES
   Shortfall;
   Excess;

MODEL

   MAX Profit =  Margin * Production
               - Margin * Shortfall
               - SellPrice * Excess;

SUBJECT TO

   BudgetLimit:    Production * ProdCost <= Budget;

   Balance:        Production + Shortfall - Excess = Demand;

END

--------------------------------------------------------------------------------------
VarCount = 3
ConCount = 2
StageCount = 2
VarStages[vars] := (1,2,2)
ConStages[cons] := (1,2)

ScenCount = 10
RandomCount = 10

TreeType = TWOSTAGE
TreeStage[scen] := (1,2,2,2,2,2,2,2,2,2)

ProbData[scen] := (0.03, 0.134, 0.086, 0.188, 0.226, 0.08, 0.056, 0.132, 0.06,  0.008)

RandomScen[random] := (1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
RandomVar[random]  := (RHS, RHS, RHS, RHS, RHS, RHS, RHS, RHS, RHS, RHS)
RandomCon[random]  := (2, 2, 2, 2, 2, 2, 2, 2, 2, 2)
RandomData[random] := (210, 290, 335, 370, 415, 450, 460, 495, 540, 585)
--------------------------------------------------------------------------------------

ext: mpl

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


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

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

(* type: work_gga_x *)

kappa := 0.8040:

den := s -> kappa + 1*MU_GE*s^2 + MU_GE^2*s^4/kappa:
f0  := s -> 1 + kappa * (1 - kappa/den(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:
(*
 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_mu := MU_GE:
params_a_b  := 0.4:
params_a_c  := 2.14951:
params_a_e  := 1.987:

vt84_gamma := 0.000023:
tpss_ff    := z -> 3:
tpss_kappa := (x, t) ->
  1/(vt84_gamma/params_a_mu^2 + vt84_gamma/params_a_mu + 1):

$include "tpss_x.mpl"

(* Equation (8) *)

vt84_f   := (x, u, t) -> 1
    + tpss_fx(x, t)*exp(-vt84_gamma*tpss_fx(x, t)/params_a_mu)/(1 + tpss_fx(x, t))
    + (1 - exp(-vt84_gamma*tpss_fx(x, t)^2/params_a_mu^2))*(params_a_mu/tpss_fx(x, t) - 1):

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


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

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

(* type: work_gga_c *)
(* prefix:
  gga_x_hjs_b88_v2_params *params;

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

params_a_a := [0.0253933, -0.0673075, 0.0891476, -0.0454168, -0.00765813, 0.0142506]:
params_a_b := [-2.6506, 3.91108, -3.31509, 1.54485, -0.198386,
  -0.136112, 0.0647862, 0.0159586, -0.000245066]:

xi := 1/(exp(20) - 1):
fs := s -> -log((exp(-s) + xi)/(1 + xi)):

$include "gga_x_hjs.mpl"

f := (rs, z, xt, xs0, xs1) ->
  f_lda(rs, z)*f1(rs, z, fs(X2S*xs0)) + f_lda(rs, -z)*f1(rs, -z, fs(X2S*xs1)):


text:
#######################################################################
# This file is part of the crlibm library, and is distributed under
# the  LGPL.
# To use:
# restart; read "log2-td.mpl";
Digits := 120:

interface(quiet=true):

read "common-procedures.mpl":
read "triple-double.mpl":
mkdir("TEMPLOG"):

# We want log2h + log2m + log2l + delta = log(2) such that
# log2h and log2m have at least 11 trailing zeros
# in order to have an exact multiplication with E, which is lower than 1024 in
# magnitude
# The resting accuracy is enough for both quick and accurate phases.

log2acc := log(2):
log2h := round(log2acc * 2**(floor(-log[2](abs(log2acc))) + (53 - 11))) / 2**(floor(-log[2](abs(log2acc))) + (53 - 11)):
log2m := round((log2acc - log2h) * 2**(floor(-log[2](abs((log2acc - log2h)))) + (53 - 11))) / 2**(floor(-log[2](abs((log2acc - log2h)))) + (53 - 11)):
log2l := log2acc - (log2h + log2m):

L := 7: # number of bits used to address the table

MAXINDEX    := round(2^L * (sqrt(2)-1)):

for i from 0 to MAXINDEX-1 do
    center[i] := 1 + i*2^(-L): # center[i] in [1, 2[
    t :=  evalf(1/center[i]):
    r[i] := round(t * 2**(floor(-log[2](abs(t))) + 23)) / 2**(floor(-log[2](abs(t))) + 23):
    (logih[i], logim[i], logil[i]) := hi_mi_lo(evalf(-log(r[i]))):
od:
for i from MAXINDEX to 2^L do
    # y has been divided by two, center[i] in [0.5, 1[
    center[i]:=(1 + i*2^(-L)) / 2:
    t :=  evalf(1/center[i]):
    r[i] := round(t * 2**(floor(-log[2](abs(t))) + 23)) / 2**(floor(-log[2](abs(t))) + 23):
    (logih[i], logim[i], logil[i]) := hi_mi_lo(evalf(-log(r[i]))):
od:

#Computation of ZMax.
for i from 0 to MAXINDEX-1 do
    __x := center[i] + 2^(-L-1) :
    zmax[i] := (__x*r[i]-1) :
    __x := center[i] - 2^(-L-1) :
    zmin[i] := (__x*r[i]-1) :
od:
for i from MAXINDEX to 2^L do
    __x := center[i] + 2^(-L-2) :
    zmax[i] := (__x*r[i]-1) :
    __x := center[i] - 2^(-L-2) :
    zmin[i] := (__x*r[i]-1) :
od:

zmaxmax:=0:
zminmin:=0:
for i from 0 to 2^L do
    tabulated_value := logih[i] + logim[i] + logil[i]:
    poly_approx_min := evalf(log(1+zmin[i])):
    poly_approx_max := evalf(log(1+zmax[i])):

    # Test if we have a case where we cancellate a lot
    # i.e. the polynomial approximation value is greater than half the tabulated value
    # the tabulated value is not exactly zero and we are of opposite sign

    if ((abs(poly_approx_min) > 0.75*abs(tabulated_value)) and (tabulated_value <> 0.0) and (poly_approx_min * tabulated_value < 0)) then
	printf("Polynomial approximation is greater in magnitude in zmin[%d] than half the tabluated value\n",i):

	printf("The tabulated value is %1.50e\n",tabulated_value):
	if (tabulated_value <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(tabulated_value)))) fi:
	printf("The value of polynomial in zmin[%d] is %1.50e\n",i,poly_approx_min):
        if (poly_approx_min <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(poly_approx_min)))) fi:

	summe := poly_approx_min + tabulated_value:
	printf("The exponent of the sum of both is %d\n",floor(log2(abs(summe)))):

    fi:

    if ((abs(poly_approx_max) > 0.75*abs(tabulated_value)) and (tabulated_value <> 0.0) and (poly_approx_max * tabulated_value <0)) then
	printf("Polynomial approximation is greater in magnitude in zmax[%d] than half the tabluated value\n",i):

	printf("The tabulated value is %1.50e\n",tabulated_value):
	if (tabulated_value <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(tabulated_value)))) fi:
	printf("The value of polynomial in zmax[%d] is %1.50e\n",i,poly_approx_max):
        if (poly_approx_max <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(poly_approx_max)))) fi:

	summe := poly_approx_max + tabulated_value:
	printf("The exponent of the sum of both is %d\n",floor(log2(abs(summe)))):

    fi:

    if zmax[i] > zmaxmax then zmaxmax := zmax[i]: fi:
    if zmin[i] < zminmin then zminmin := zmin[i]: fi:
od:
printf("zminmin = -2^(%2f)   zmaxmax = 2^(%2f)\n", log2(-zminmin), log2(zmaxmax) ) :

PolyDegreeQuick:=7:

printf("   degree of the polynomial used in the quick phase is %d\n",PolyDegreeQuick);

DDNumberQuick:=3:

printf("   number of double doubles used for the coefficients is %d\n",DDNumberQuick);

#Keep -zmaxmax..zmaxmax to keep c1=1, which is useful in the proof
#and constrain the first two coefficients to 1 and -1/2 in order to save up a full multiplication and a rounding error
polyQuick:= poly_exact2(x*(1+x*(-0.5+x*(numapprox[minimax]((((log(1+x)/x)-1)/x+0.5)/x,
			x=-zmaxmax..zmaxmax,  [PolyDegreeQuick-3,0], 1 ,  'deltaApprox')))), DDNumberQuick):

#Try to verify the bound for using double double arithmetic.
#Idea: compare the maximum absolute value of the polynomial in zmaxmax (the polynomial has its maxima at the limits)
#with the maximum value of the first term which is calculated in double precision only

p := unapply(polyQuick,x):
printf("   using only %d double doubles should be fine since the hi z ulp should affect the result starting from bit %f\n",
	DDNumberQuick,evalf(53 + log2(p(zmaxmax)) - log2(zmaxmax^(DDNumberQuick)),5)):

epsilonApproxQuick := numapprox[infnorm]( 1-polyQuick/log(1+x), x=zminmin..zmaxmax):
printf("   approximation rel error for the quick phase is 2^(%2f)\n", log2(epsilonApproxQuick) ) :
deltaApproxQuick := numapprox[infnorm]( polyQuick-log(1+x), x=zminmin..zmaxmax):
printf("   approximation abs error for the quick phase is 2^(%2f)\n", log2(deltaApproxQuick) ) :

PolyDegreeAccurate:=14:

printf("   degree of the polynomial used in the accurate phase is %d\n",PolyDegreeAccurate);

DDNumberAccu:=7:
TDNumberAccu:=3:

printf("   number of triple doubles used for the coefficients is %d\n",TDNumberAccu);
printf("   number of double doubles used for the coefficients is %d\n",DDNumberAccu);

#Keep -zmaxmax..zmaxmax to keep c1=1, which is useful in the proof
polyAccurate:= poly_exact32(x*(1+x*(-0.5+x*(numapprox[minimax]((((log(1+x)/x)-1)/x+0.5)/x,
				x=-zmaxmax..zmaxmax,  [PolyDegreeAccurate-3,0], 1 ,  'deltaApprox')))),
				TDNumberAccu, DDNumberAccu):

#Try to verify the bound for using double double arithmetic.
#Idea: compare the maximum absolute value of the polynomial in zmaxmax (the polynomial has its maxima at the limits)
#with the maximum value of the first term which is calculated in double precision only

pp := unapply(polyAccurate,x):
printf("   using only %d triple doubles should be fine since the mi z ulp should affect the result starting from bit %f\n",
	TDNumberAccu,evalf(106 + log2(p(zmaxmax)) - log2(zmaxmax^(TDNumberAccu)),5)):
printf("   using only %d double doubles should be fine since the hi z ulp should affect the result starting from bit %f\n",
	DDNumberAccu,evalf(53 + log2(p(zmaxmax)) - log2(zmaxmax^(TDNumberAccu + DDNumberAccu)),5)):

epsilonApproxAccurate := numapprox[infnorm]( 1-polyAccurate/log(1+x), x=zminmin..zmaxmax):
printf("   approximation rel error for the accurate phase is 2^(%2f)\n", log2(epsilonApproxAccurate) ) :
deltaApproxAccurate := numapprox[infnorm]( polyAccurate-log(1+x), x=zminmin..zmaxmax):
printf("   approximation abs error for the quick phase is 2^(%2f)\n", log2(deltaApproxAccurate) ) :

#Compute now the inverse of ln(10) for the final multiplication of ln(x) with this constant for obtaining log10(x)
#Compute also the relative error of the constant stored as a triple double.
Log10inv := evalf(1 / log(10)):

(log10invh, log10invm, log10invl) := hi_mi_lo(Log10inv):

Log10invhml := log10invh + log10invm + log10invl:

epsilonLog10invhml := evalf(abs((Log10invhml - Log10inv) / Log10inv)):

printf("   Log10inv = 1 / ln(10) stored as a triple-double is exact with a relative error of 2^(%2f)\n",
	evalf(log[2](epsilonLog10invhml))):

#Compute now the constant needed for the unfiltered directed final rounding of the triple-double result
#This constant is supposed to be the reciprocal of the critical accuracy of the function.
#We suppose this critical accuracy to be 2^(-120) because the worst case (for RD) we know of is x = 403ce41d 8fa665fa
#We can easily spend some guard bits since we want simply filter out cases with a theoretical worst case accuracy
#of -infty (exact floating point images) and as we know that we are exact to at least 122 bits.

wca := 2^(122):

#-------------------------------------------------------------------
# Output

filename:="TEMPLOG/log10-td.h":
fd:=fopen(filename, WRITE, TEXT):

fprintf(fd, "#include \"crlibm.h\"\n#include \"crlibm_private.h\"\n"):

fprintf(fd, "\n/*File generated by maple/log10-td.mpl*/\n"):

fprintf(fd, "\n\#define L %d\n\n",L):
fprintf(fd, "\#define MAXINDEX %d\n\n", MAXINDEX):
fprintf(fd, "\#define INDEXMASK %d\n", 2^L-1):
fprintf(fd, "\#define two52 %1.50e\n", 2^(52)):
fprintf(fd, "\#define log2h %1.50e\n", log2h):
fprintf(fd, "\#define log2m %1.50e\n", log2m):
fprintf(fd, "\#define log2l %1.50e\n", log2l):
fprintf(fd, "\#define log10invh %1.50e\n",log10invh):
fprintf(fd, "\#define log10invm %1.50e\n",log10invm):
fprintf(fd, "\#define log10invl %1.50e\n",log10invl):
fprintf(fd, "\#define WORSTCASEACCURACY %1.50e\n",wca):

epsilon_quick_1 := 2^(-61): # The Gappa proof will show this bound
epsilon_quick_2 := 2^(-61): # The Gappa proof will show this bound
fprintf(fd, "\#define ROUNDCST1 %1.50e\n", compute_rn_constant(epsilon_quick_1)):
fprintf(fd, "\#define ROUNDCST2 %1.50e\n", compute_rn_constant(epsilon_quick_2)):
fprintf(fd, "\#define RDROUNDCST1 %1.50e\n", epsilon_quick_1):
fprintf(fd, "\#define RDROUNDCST2 %1.50e\n", epsilon_quick_2):

fprintf(fd, "\n\n"):

# Print the defines for the define statements

for i from 3 to PolyDegreeQuick do
	fprintf(fd, "\#define c%d %1.50e\n",i,coeff(polyQuick,x,i)):
od:

fprintf(fd, "\n\n"):

for i from 3 to (DDNumberAccu + TDNumberAccu -1) do
	(hi,lo) := hi_lo(coeff(polyAccurate,x,i)):
	fprintf(fd, "\#define accPolyC%dh %1.50e\n",i,hi):
	fprintf(fd, "\#define accPolyC%dl %1.50e\n",i,lo):
od:

for i from (DDNumberAccu + TDNumberAccu) to PolyDegreeAccurate do
	fprintf(fd, "\#define accPolyC%d %1.50e\n",i,coeff(polyAccurate,x,i)):
od:

fprintf(fd, "\n\n"):

# Print the table
fprintf(fd, "typedef struct rri_tag {float ri; double logih; double logim; double logil;} rri;  \n"):
fprintf(fd, "static const rri argredtable[%d] = {\n", 2^L):
for i from 0 to 2^L-1 do
      fprintf(fd, "  { \n"):
      fprintf(fd, "    %1.50e,   /* r[%d] */ \n", r[i], i):
      fprintf(fd, "    %1.50e, /* logih[%d] */ \n", logih[i], i):
      fprintf(fd, "    %1.50e, /* logim[%d] */ \n", logim[i], i):
      fprintf(fd, "    %1.50e, /* logil[%d] */ \n", logil[i], i):
      fprintf(fd, "  } "):
      if(i<2^L-1) then  fprintf(fd, ", \n"): fi
od:
fprintf(fd, "}; \n \n"):

fclose(fd):

for j from 0 to 2^L-1 do
    filename:=cat("TEMPLOG/log10-td_",j,".sed"):
    fd:=fopen(filename, WRITE, TEXT):
    fprintf(fd, " s/_log2h/%1.50e/g\n", log2h):
    fprintf(fd, " s/_log2m/%1.50e/g\n", log2m):
    fprintf(fd, " s/_log2l/%1.50e/g\n", log2l):
    fprintf(fd, " s/_logih/%1.50e/g\n", logih[j]):
    fprintf(fd, " s/_logim/%1.50e/g\n", logim[j]):
    fprintf(fd, " s/_logil/%1.50e/g\n", logil[j]):
    fprintf(fd, " s/_zmin/%1.50e/g\n", zmin[j]):
    fprintf(fd, " s/_zmax/%1.50e/g\n", zmax[j]):
    for i from 3 to PolyDegreeQuick do
        fprintf(fd, " s/_c%d/%1.50e/g\n", i, coeff(polyQuick,x,i)):
    od:
    fprintf(fd, " s/_epsilonApproxQuick/%1.50e/g\n", epsilonApproxQuick):
    fprintf(fd, " s/_epsilonLog10invhml/%1.50e/g\n", epsilonLog10invhml):
    fclose(fd):
  od:

for j from 0 to 2^L-1 do
    filename:=cat("TEMPLOG/log10-td-accurate_",j,".sed"):
    fd:=fopen(filename, WRITE, TEXT):
    fprintf(fd, " s/_log2h/%1.50e/g\n", log2h):
    fprintf(fd, " s/_log2m/%1.50e/g\n", log2m):
    fprintf(fd, " s/_log2l/%1.50e/g\n", log2l):
    fprintf(fd, " s/_logih/%1.50e/g\n", logih[j]):
    fprintf(fd, " s/_logim/%1.50e/g\n", logim[j]):
    fprintf(fd, " s/_logil/%1.50e/g\n", logil[j]):
    fprintf(fd, " s/_zmin/%1.50e/g\n", zmin[j]):
    fprintf(fd, " s/_zmax/%1.50e/g\n", zmax[j]):
    for i from 3 to (DDNumberAccu + TDNumberAccu -1) do
	(hi,lo) := hi_lo(coeff(polyAccurate,x,i)):
        fprintf(fd, " s/_accPolyC%dh/%1.50e/g\n", i, hi):
        fprintf(fd, " s/_accPolyC%dl/%1.50e/g\n", i, lo):
    od:
    for i from (DDNumberAccu + TDNumberAccu) to PolyDegreeAccurate do
        fprintf(fd, " s/_accPolyC%d/%1.50e/g\n", i, coeff(polyAccurate,x,i)):
    od:
    fprintf(fd, " s/_epsilonApproxAccurate/%1.50e/g\n", epsilonApproxAccurate):
    fprintf(fd, " s/_epsilonLog10invhml/%1.50e/g\n", epsilonLog10invhml):
    fclose(fd):
  od:

# A shell script to use them
filename:="run-log10-td-proof.sh":
fd:=fopen(filename, WRITE, TEXT):
fprintf(fd, "#!/bin/sh\n"):
fprintf(fd, "# You probably need to edit the path to the gappa executable\n"):
fprintf(fd, "GAPPA=~/ble/gappa-0.4.5/src/gappa\n"):
fprintf(fd, "# Test all the possible table value for E=1\n"):
fprintf(fd, "for num in `seq 0 %d`; do\n", 2^L-1):
fprintf(fd, "  echo $num, E=1:\n"):
fprintf(fd, "  sed -f ./TEMPLOG/log10-td_$num.sed log10-td.gappa | $GAPPA > /dev/null\n"):
fprintf(fd, "  echo\n"):
fprintf(fd, "done\n"):
fprintf(fd, "# For the case E=0 we first handle the cases 0 and %d using log10-td-E0-logir0.gappa\n", 2^L):
fprintf(fd, "echo 0 and %d, E=0:\n", 2^L):
fprintf(fd, "sed -f log10-td_0.sed log10-td-E0-logir0.gappa | $GAPPA > /dev/null\n"):
fprintf(fd, "# then the other cases where logirh <> 0\n"):
fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
fprintf(fd, "  echo $num, E=0:\n"):
fprintf(fd, "  sed -f ./TEMPLOG/log10-td_$num.sed log10-td-E0.gappa | $GAPPA > /dev/null\n"):
fprintf(fd, "  echo\n"):
fprintf(fd, "done\n"):
fprintf(fd, "# Accurate phase: Test all the possible table value for E=1\n"):
fprintf(fd, "for num in `seq 0 %d`; do\n", 2^L-1):
fprintf(fd, "  echo Accurate phase: $num, E=1:\n"):
fprintf(fd, "  sed -f ./TEMPLOG/log10-td-accurate_$num.sed log10-td-accurate.gappa | $GAPPA > /dev/null\n"):
fprintf(fd, "  echo\n"):
fprintf(fd, "done\n"):
fprintf(fd, "# Accurate phase: For the case E=0 we first handle the cases 0 and %d using log10-td-accurate-E0-logir0.gappa\n", 2^L):
fprintf(fd, "echo 0 and %d, E=0:\n", 2^L):
fprintf(fd, "sed -f ./TEMPLOG/log10-td-accurate_0.sed log10-td-accurate-E0-logir0.gappa | $GAPPA > /dev/null\n"):
fprintf(fd, "# Accurate phase: then the other cases where logirh <> 0\n"):
fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
fprintf(fd, "  echo $num, E=0:\n"):
fprintf(fd, "  sed -f ./TEMPLOG/log10-td-accurate_$num.sed log10-td-accurate-E0.gappa | $GAPPA > /dev/null\n"):
fprintf(fd, "  echo\n"):
fprintf(fd, "done\n"):
fclose(fd):

printf("----DONE---\n") :

max_stars_repo_name: squarePenguin/parvsl

max_issues_repo_name: squarePenguin/parvsl

max_forks_repo_name: squarePenguin/parvsl


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

gvt4_alpha  := 0.00186726:
gvt4_coeff_d := [-9.800683e-01, -3.556788e-03, 6.250326e-03, -2.354518e-05, -1.282732e-04, 3.574822e-04]:

$include "gvt4.mpl"

gvt4_f := (x, u, t) -> -gtv4(gvt4_alpha, gvt4_coeff_d, x, 2*(t - K_FACTOR_C))/X_FACTOR_C:

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


text:
kernelopts(printbytes=false);
with(RigidMotionsRecoverNMM);
LaunchFindDistinctSamplePoints([a,b,c], "N2", [-2, -1, 0, 1, 2], "__REPLACE__");
done

max_stars_repo_name: copyme/MapleTools

max_issues_repo_name: copyme/MapleTools

max_forks_repo_name: copyme/MapleTools


text:
######################################################
# Performing the compactificatoin and non-uniform
# coordinate transformation on the tensor expressions
######################################################

all_pdes:=[
   [DPH,0,dn] , [dg,2,dn] ,
         [DDL,2,dn],
   [DTK,0,dn],
         [RR,2,dn],  [P,2,dn],
   [DLAM,1,up],
       [C1,0,dn],
   [HS,0,dn],
   [DB,0,dn],
   [HS3,0,dn],
   [MK,1,up],
   [EQRPIb2,0,dn],
   [EQIPIb2,0,dn],
   [DRPSI,0,dn],
   [DIPSI,0,dn]
             ];

# Going from x -> ctfm(x) using DEtools
with(DEtools);
addcoords(compact,[t,x],[t,ctfm(x)]);

# See built_eqns.mpl for how 'tt' table is defined

for ii from 1 to nops(all_pdes) do
   nm := all_pdes[ii][1]:
   rank := all_pdes[ii][2]:
   tp := all_pdes[ii][3]:
   printf("Transforming %a\n",nm);
   if rank <> 0 then

    for jj from 1 to nops(coord) do
     if (tt[nm][all][coord[jj]] <> 0 ) then
      tt[nm][allcomp][coord[jj]] := PDEchangecoords(tt[nm][all][coord[jj]],[t,x],compact,[t,x]);
      tt[nm][regcomp][coord[jj]] := PDEchangecoords(tt[nm][reg][coord[jj]],[t,x],compact,[t,x]);

	   # Symbolically replacing derivatives of compactification function
      # since we do not want to finite difference them)

      tt[nm][allcomp][coord[jj]] := subs({diff(ctfm(x),x)=ctfmp(x),diff(ctfm(x),x,x)=ctfmpp(x)},tt[nm][allcomp][coord[jj]]);
      tt[nm][regcomp][coord[jj]] := subs({diff(ctfm(x),x)=ctfmp(x),diff(ctfm(x),x,x)=ctfmpp(x)},tt[nm][regcomp][coord[jj]]);

      # DEtools will replace x with ctfm(x) symbolically in our expressions, but ctfm(x) is our new x, fixing that
      for kk from 1 to nops(grid_functions) do
        fn := grid_functions[kk];
        tt[nm][allcomp][coord[jj]] := subs(fn(t,ctfm(x))=fn(t,x),tt[nm][allcomp][coord[jj]]);
        tt[nm][regcomp][coord[jj]] := subs(fn(t,ctfm(x))=fn(t,x),tt[nm][regcomp][coord[jj]]);
      end do;

     end if;

    end do;
   else

     # Same as above, symbolically replacing derivatives of compactification function
     if ( tt[nm][all] <> 0 ) then
      tt[nm][allcomp] := PDEchangecoords(tt[nm][all],[t,x],compact,[t,x]);
      tt[nm][regcomp] := PDEchangecoords(tt[nm][reg],[t,x],compact,[t,x]);
      tt[nm][allcomp] := subs({diff(ctfm(x),x)=ctfmp(x),diff(ctfm(x),x,x)=ctfmpp(x)},tt[nm][allcomp]);
      tt[nm][regcomp] := subs({diff(ctfm(x),x)=ctfmp(x),diff(ctfm(x),x,x)=ctfmpp(x)},tt[nm][regcomp]);

      # Same as above, fixing DEtools mess!
      for kk from 1 to nops(grid_functions) do
       fn := grid_functions[kk];
       tt[nm][allcomp] := subs(fn(t,ctfm(x))=fn(t,x),tt[nm][allcomp]);
       tt[nm][regcomp] := subs(fn(t,ctfm(x))=fn(t,x),tt[nm][regcomp]);
      end do;

     end if;

   end if;
end do:

# Checking if everything is transformed properly at the limit ctfm(x) = x

for ii from 1 to nops(all_pdes) do
   nm := all_pdes[ii][1]:
   rank := all_pdes[ii][2]:
   tp := all_pdes[ii][3]:
   if rank <> 0 then

			 for jj from 1 to nops(coord) do
            if ( tt[nm][all][coord[jj]] <> 0 ) then
				res := simplify(tt[nm][all][coord[jj]] - eval(tt[nm][allcomp][coord[jj]],{ctfm(x) = x, ctfmp(x) = 1, ctfmpp(x) = 0}));
				printf("res = %a\n",res);
				res := simplify(tt[nm][reg][coord[jj]] - eval(tt[nm][regcomp][coord[jj]],{ctfm(x) = x, ctfmp(x) = 1, ctfmpp(x) = 0}));
				printf("res = %a\n",res);
            end if;
			 end do;

   else
      res := simplify(tt[nm][all] - eval(tt[nm][allcomp],{ctfm(x) = x, ctfmp(x) = 1, ctfmpp(x) = 0}));
      printf("res = %a\n",res);
      res := simplify(tt[nm][reg] - eval(tt[nm][regcomp],{ctfm(x) = x, ctfmp(x) = 1, ctfmpp(x) = 0}));
      printf("res = %a\n",res);
   end if;
end do:

##################################################
# Performing a collection operation, i.e.
# powers of ctfm(x) function, this will help the
# intel c and fortran compiler a lot!
#################################################

for ii from 1 to nops(all_pdes) do
   nm := all_pdes[ii][1]:
   rank := all_pdes[ii][2]:
   tp := all_pdes[ii][3]:
   if rank <> 0 then

			 for jj from 1 to nops(coord) do
            if (tt[nm][all][coord[jj]] <> 0 ) then
  			   tt[nm][all][coord[jj]] := collect(collect(eval(expand(tt[nm][allcomp][coord[jj]]),{ctfmpp(x)=octfmp(x)*ctfmp(x)^2}),1/ctfmp(x)),1/ctfm(x));
			   tt[nm][reg][coord[jj]] := collect(collect(eval(expand(tt[nm][regcomp][coord[jj]]),{ctfmpp(x)=octfmp(x)*ctfmp(x)^2}),1/ctfmp(x)),1/ctfm(x));
				print(tt[nm][all][coord[jj]]);
            print(tt[nm][reg][coord[jj]]);
            end if;
			 end do;

   else
      tt[nm][all] :=  collect(collect(eval(expand(tt[nm][allcomp]),{ctfmpp(x)=octfmp(x)*ctfmp(x)^2}),1/ctfmp(x)),1/ctfm(x));
      tt[nm][reg] :=  collect(collect(eval(expand(tt[nm][regcomp]),{ctfmpp(x)=octfmp(x)*ctfmp(x)^2}),1/ctfmp(x)),1/ctfm(x));
	   print(tt[nm][all]);
	   print(tt[nm][reg]);
   end if;
end do:

# Adding new members of grid functions to our set:
grid_functions := grid_functions union {ctfm,ctfmp,octfmp};


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

# Complex numbers

`is_element/CC` := (x) -> type(x,constant):

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

`is_leq/CC` := NULL;

`is_leq_lex/CC` := proc(x,y)
 local z;
 z := simplify(y - x);
 if is(Re(z) > 0) then
  return true;
 elif is(Re(z) < 0) then
  return false;
 elif is(Im(z) > 0) then
  return true;
 elif is(Im(z) < 0) then
  return false;
 fi;
 return true;
end:

`random_element/CC` := proc()
 rand(-720..720)()/(360) + rand(-720..720)()/(360) * I;
end:

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

######################################################################
# C(N) is the vector space CC^N

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

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

`is_leq/C` := NULL;

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

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

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

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

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

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

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


text:
# Example from Section 5.2 of the paper

read "../ComputeIdentifiableFunctions.mpl";

model := [
     diff(xa(t), t) = -k1 * xa(t),
     diff(xb(t), t) = k1 * xa(t) - k2 * xb(t),
     diff(xc(t), t) = k2 * xb(t),
     diff(eA(t), t) = 0,
     diff(eC(t), t) = 0,
     y1(t) = eA(t) * xa(t) + eB * xb(t) + eC(t) * xc(t),
     y2(t) = xc(t),
     y3(t) = eA(t),
     y4(t) = eC(t)
]:

# Computing function identifiable from a single experiment
printf("Single-experiment identifiable functions: %a\n", SingleExperimentIdentifiableFunctions(model)):

me := MultiExperimentIdentifiableFunctions(model, simplified_generators=true):
printf("The multi-experiment identifiable functions are %a; This is not the same as for a single experiment!\n", me[3]):
printf("The bound from the theorem is %a showing that %a experiments are sufficient to identify these functions\n", me[1], me[1]):

max_stars_repo_name: iliailmer/AllIdentifiableFunctions

max_issues_repo_name: iliailmer/AllIdentifiableFunctions

max_forks_repo_name: iliailmer/AllIdentifiableFunctions


text:
type $atyp <[10] f32>
var $fff <i32> = 7789
func $foo ( var %i i32 )i32
type $i32ty <struct {
      @FF1 i32,
      @FF2 f64 }>

func $fact (
 var %n i32, var %m f32) i32 {
 var %alocal <$i32ty>
 dassign %alocal (constval i32 999)
 if (ne i32 i32 (dread i32 %n, constval i32 1)) {
    eval (dread i32 %n)
    call &foo(
      add i32 (dread i32 %n, constval i32 1))
    return (regread i32 %%retval)}
   else {
     return (constval i32 1) } }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

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

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

f := (rs, zeta) -> params_a_C1 + params_a_C2*n_total(rs)^(-1/3) + params_a_C3*n_total(rs)^(-2/3):


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

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

(* type: work_gga_x *)

$define gga_x_rpw86_params
$include "gga_x_pw86.mpl"

malpha := 0.02178:
mbeta  := 1.15:
muLV   := 0.8491/9:

f0 := s ->
   + (1 + muLV*s^2)/(1 + malpha*s^6)
   + malpha*s^6*f0_pw86(s)/(mbeta + malpha*s^6):

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


text:
 func $retypef64toi64 (
  var %i f64
  ) i64 {
   return (
     retype i64 <i64> (dread f64 %i))}

 func $retypei64tof64 (
  var %i i64
  ) f64 {
   return (
     retype f64 <f64> (dread i64 %i))}

 func $retypef32toi32 (
  var %i f32
  ) i32 {
   return (
     retype i32 <i32> (dread f32 %i))}

 func $retypei32tof32 (
  var %i i32
  ) f32 {
   return (
     retype f32 <f32> (dread i32 %i))}
 # 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:
# 基于Maple read 命令的代码读取包。
# 优点：
#     1. Maple原生支持，不需要额外的文件，直接被mla带走。
#     2. 不用产生临时文件。
# 缺点：
#     1. 预处理器必须在行首。
#     2. include的相对路径必须以 . 或 .. 开头。
#     3. 文档编码必须是UTF-8。
MapleCodeReader:=module()
    option package;
    export ReadCode;

    ReadCode:=proc(fname)
        local cd,p,f;
        p:=FileTools:-ParentDirectory(fname);
        f:=FileTools:-Filename(fname);
        cd:=currentdir();
        currentdir(p);
        read(f);
        currentdir(cd);
        return;
    end proc:
end module:


text:
flavor 1
srclang 3
id 65535
numfuncs 2
import "huaweitest.mplt"
import "/home/bravewtz/Desktop/openarkcompiler/libjava-core/libjava-core.mplt"
entryfunc &Lhuaweitest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V
fileinfo {
  @INFO_filename "huaweitest.jar"}
srcfileinfo {
  1 "huaweitest.java"}
type $__class_meta__ <struct {
  @shadow <* void> final,
  @monitor i32 final,
  @classloader u16 final,
  @objsize u16 final,
  @itab <* void> final,
  @vtab <* void> final,
  @gctib <* void> final,
  @classinforo <* void> final,
  @clinitbridge <* void> final}>
type $__class_meta_ro__ <struct {
  @classname <* void> final,
  @ifields <* void> final,
  @methods <* void> final,
  @superclass_or_componentclass <* void> final,
  @numoffields u16 final,
  @numofmethods u16 final,
  @flag u16 final,
  @numofsuperclasses u16 final,
  @padding u32 final,
  @mod i32 final,
  @annotation i32 final,
  @clinitAddr i32 final}>
type $__method_info__ <struct {
  @method_in_vtab_index i64 final,
  @declaringclass i64 final,
  @addr i64 final,
  @mod i32 final,
  @methodname i32 final,
  @signaturename i32 final,
  @annotationvalue i32 final,
  @flag u16 final,
  @argsize u16 final,
  @padding u32 final}>
type $__method_info_compact__ <struct {
  @method_in_vtab_index i32 final,
  @addr i32 final,
  @lebPadding0 u8 final}>
type $__field_info__ <struct {
  @offset u64 final,
  @mod u32 final,
  @flag u16 final,
  @index u16 final,
  @typeName i64 final,
  @fieldname u32 final,
  @annotation u32 final,
  @declaringclass <* void> final}>
type $__field_info_compact__ <struct {
  @offset u32 final,
  @lebPadding0 u8 final}>
type $__superclass_meta__ <struct {
  @superclassinfo <* void> final}>
type $MUIDFuncDefTabEntry <struct {
  @funcUnifiedAddr u64}>
type $MUIDFuncInfTabEntry <struct {
  @funcSize u32,
  @funcName u32}>
type $MUIDFuncDefMuidTabEntry <struct {
  @muidLow u64,
  @muidHigh u64}>
type $MUIDDataDefTabEntry <struct {
  @dataUnifiedAddr u64}>
type $MUIDDataDefMuidTabEntry <struct {
  @muidLow u64,
  @muidHigh u64}>
type $MUIDUnifiedUndefTabEntry <struct {
  @globalAddress u64}>
type $MUIDUnifiedUndefMuidTabEntry <struct {
  @muidLow u64,
  @muidHigh u64}>
type $MUIDRangeTabEntry <struct {
  @tabBegin <* void>,
  @tabEnd <* void>}>
javaclass $Lhuaweitest_3B <$Lhuaweitest_3B> public
var $_C_STR_907ba9a10d33c74f1c37be3b41544334 fstatic <[10] u64> readonly = [0, 0x6b00000000, 0x20490949984357f2, 0x4909490909490920, 0x909494949494949, 0x2020202020202049, 0x2020200909492020, 0x4949494920204920, 0x4949202049494949, 73]
var $_C_STR_194a027f7d649763903a7dcda50cf5a0 fstatic <[9] u64> readonly = [0, 0x5b00000000, 0x94909491faf02dc, 0x2009490949090949, 0x2020492009094920, 0x4920492020202020, 0x4920202020202020, 0x2020090949202020, 73]
var $_C_STR_901191d3899ab5fc7be3d97fe762f978 fstatic <[9] u64> readonly = [0, 0x5f00000000, 0x49494949fbe6d2fb, 0x4909490909490949, 0x2020090949202009, 0x2049202020202049, 0x2020202020492020, 0x909492020200949, 0x492020]
var $_C_STR_4a20a28a5a71422be90661b1072ca29e fstatic <[9] u64> readonly = [0, 0x6500000000, 0x94909494f3904e4, 0x2009490949090949, 0x4920202009094920, 0x2020202049202020, 0x2009492020204920, 0x4949494949492020, 0x492020204949]
var $_C_STR_e4c52869a15445bee0db530b07e9b98a fstatic <[8] u64> readonly = [0, 0x5100000000, 0x949094949ddeb72, 0x2009490949090949, 0x4920490909094920, 0x4920202020202020, 0x949202020094920, 0x49202009]
var $_C_STR_58a8c3b0d70b92b64324906ef9c02c90 fstatic <[10] u64> readonly = [0, 0x6d00000000, 0x94909497041ce40, 0x4949494949494949, 0x4949494949490949, 0x2049200909494949, 0x2020202020202020, 0x4949492020200949, 0x4920094949494949, 0x4949]
var $Ljava_2Flang_2FSystem_3B_7Cout extern <* <$Ljava_2Fio_2FPrintStream_3B>> final public static
var $__cinf_Ljava_2Flang_2FString_3B <$__class_meta__>
func &MCC_GetOrInsertLiteral () <* <$Ljava_2Flang_2FString_3B>>
var $__vtb_Lhuaweitest_3B fstatic <[11] <* void>> = [16, 24, 32, 36, 8, 4, 48, 12, 20, 44, 28]
var $__cinf_Lhuaweitest_3B <$__class_meta__> public
var $__methods_info__Lhuaweitest_3B fstatic <[2] <$__method_info__>> public = [[1= 0xfff6, 2= addrof ptr $__cinf_Lhuaweitest_3B, 3= addroffunc ptr &Lhuaweitest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V, 4= 9, 5= 56, 6= 76, 7= 168, 8= 0xa201, 9= 1, 10= 0], [1= 0xfff6, 2= addrof ptr $__cinf_Lhuaweitest_3B, 3= addroffunc ptr &Lhuaweitest_3B_7C_3Cinit_3E_7C_28_29V, 4= 0x10001, 5= 184, 6= 212, 7= 168, 8= 0xad81, 9= 1, 10= 0]]
var $__cinf_Ljava_2Flang_2FObject_3B extern <$__class_meta__> public
var $__superclasses__Lhuaweitest_3B fstatic <[1] <$__superclass_meta__>> public = [[1= 0x4000000000000000]]
var $__classinforo__Lhuaweitest_3B fstatic <$__class_meta_ro__> public = [1= 4, 2= 0, 3= addrof ptr $__methods_info__Lhuaweitest_3B, 4= addrof ptr $__superclasses__Lhuaweitest_3B, 5= 0, 6= 2, 7= 0, 8= 1, 9= 0, 10= 33, 11= 168, 12= 0]
var $MCC_GCTIB__Lhuaweitest_3B fstatic <* void> public
var $classStateInitialized u64
var $__cinf_Lhuaweitest_3B <$__class_meta__> public = [1= 0x1bf435e944cf11ae, 2= 0, 3= 0xffff, 4= 0, 5= 0, 6= addrof ptr $__vtb_Lhuaweitest_3B, 7= addrof ptr $MCC_GCTIB__Lhuaweitest_3B, 8= addrof ptr $__classinforo__Lhuaweitest_3B, 9= addrof ptr $classStateInitialized]
var $__muid_classmetadata_bucket$$huaweitest_jar <[1] <* void>> public = [addrof ptr $__cinf_Lhuaweitest_3B]
func &MCC_Reflect_ThrowCastException nosideeffect () void
func &MCC_Reflect_Check_Casting_NoArray nosideeffect () void
func &MCC_Reflect_Check_Casting_Array nosideeffect () void
func &MCC_CheckThrowPendingException nosideeffect () void
func &MCC_PreNativeCall (var %caller ref) <* void>
func &MCC_PostNativeCall (var %env <* void>) void
func &MCC_DecodeReference nosideeffect (var %obj ref) ref
func &MCC_CallFastNative (var %func <* void>) <* void>
func &MCC_CallSlowNative0 (var %func <* void>) <* void>
func &MCC_CallSlowNative1 (var %func <* void>) <* void>
func &MCC_CallSlowNative2 (var %func <* void>) <* void>
func &MCC_CallSlowNative3 (var %func <* void>) <* void>
func &MCC_CallSlowNative4 (var %func <* void>) <* void>
func &MCC_CallSlowNative5 (var %func <* void>) <* void>
func &MCC_CallSlowNative6 (var %func <* void>) <* void>
func &MCC_CallSlowNative7 (var %func <* void>) <* void>
func &MCC_CallSlowNative8 (var %func <* void>) <* void>
func &MCC_CallFastNativeExt (var %func <* void>) <* void>
func &MCC_CallSlowNativeExt (var %func <* void>) <* void>
func &MCC_SetReliableUnwindContext nosideeffect () void
var $__reg_jni_func_tab$$huaweitest_jar <[0] <* void>>
var $__cinf_Ljava_2Flang_2FSystem_3B extern <$__class_meta__>
func &MCC_getFuncPtrFromItabSecondHash64 nosideeffect () ptr
var $__muid_func_def_tab$$huaweitest_jar fstatic <[2] <$MUIDFuncDefTabEntry>> = [[1= addroffunc ptr &Lhuaweitest_3B_7C_3Cinit_3E_7C_28_29V], [1= addroffunc ptr &Lhuaweitest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V]]
var $__muid_func_inf_tab$$huaweitest_jar fstatic <[2] <$MUIDFuncInfTabEntry>> = [[1= addroffunc ptr &Lhuaweitest_3B_7C_3Cinit_3E_7C_28_29V, 2= addroffunc ptr &Lhuaweitest_3B_7C_3Cinit_3E_7C_28_29V], [1= addroffunc ptr &Lhuaweitest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V, 2= addroffunc ptr &Lhuaweitest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V]]
var $__muid_func_def_muid_tab$$huaweitest_jar fstatic <[2] <$MUIDFuncDefMuidTabEntry>> = [[1= -1722806870279642204, 2= -946015414860200947], [1= -1343471055290971758, 2= -2119983759179258436]]
var $__muid_func_muid_idx_tab$$huaweitest_jar fstatic <[2] u32> = [1, 0]
var $__muid_data_def_tab$$huaweitest_jar fstatic <[1] <$MUIDDataDefTabEntry>> = [[1= addrof ptr $__cinf_Lhuaweitest_3B]]
var $__muid_data_def_muid_tab$$huaweitest_jar fstatic <[1] <$MUIDDataDefMuidTabEntry>> = [[1= -6287124398786487729, 2= -1011506809175871035]]
var $__muid_func_undef_tab$$huaweitest_jar fstatic <[12] <$MUIDUnifiedUndefTabEntry>> = [[1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0]]
var $__muid_func_undef_muid_tab$$huaweitest_jar fstatic <[12] <$MUIDUnifiedUndefMuidTabEntry>> = [[1= 0x3e32352aee789835, 2= -3887705395317205813], [1= -4187412136968710015, 2= -3800091941095621250], [1= 0x6742c234127e0a27, 2= -3762262047879347071], [1= 0x7230554331c55d92, 2= -3676689525926909155], [1= -715372855679083712, 2= -2647497990906227723], [1= -7464356948810446352, 2= -2259485500590180091], [1= 0x783627f2afd1cbde, 2= -2046851302095768916], [1= -2701934576591406938, 2= -1693831364093527548], [1= 0x126777a7fe39e1fb, 2= -1314856249532362766], [1= 0xbf40578f3343f7a, 2= -1198421541845410999], [1= 0x7ca2bdf69a6c7c94, 2= -801329978528900548], [1= 0x477aafa4d7dd102b, 2= -442561182569419835]]
var $__muid_data_undef_tab$$huaweitest_jar fstatic <[3] <$MUIDUnifiedUndefTabEntry>> = [[1= addrof ptr $__cinf_Ljava_2Flang_2FObject_3B], [1= addrof ptr $Ljava_2Flang_2FSystem_3B_7Cout], [1= addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B]]
var $__muid_data_undef_muid_tab$$huaweitest_jar fstatic <[3] <$MUIDUnifiedUndefMuidTabEntry>> = [[1= -567417612161374449, 2= -3298852447504547670], [1= 0x191283ac418c4bb9, 2= -1676204161023949463], [1= -5921653145571052587, 2= -171150348656858163]]
var $__muid_range_tab$$huaweitest_jar fstatic <[29] <$MUIDRangeTabEntry>> = [[1= 0x5b15548b999b8d93, 2= -222031886891092138], [1= -5449160688458699489, 2= -3307888487811147755], [1= 2, 2= 2], [1= 3, 2= 3], [1= 4, 2= 4], [1= 5, 2= 5], [1= 6, 2= 6], [1= 7, 2= 7], [1= 8, 2= 8], [1= 9, 2= 9], [1= 0, 2= 0], [1= 11, 2= 11], [1= 12, 2= 12], [1= 13, 2= 13], [1= 14, 2= 14], [1= 15, 2= 15], [1= addrof ptr $__muid_func_def_tab$$huaweitest_jar, 2= addrof ptr $__muid_func_def_tab$$huaweitest_jar], [1= 0, 2= 0], [1= addrof ptr $__muid_func_inf_tab$$huaweitest_jar, 2= addrof ptr $__muid_func_inf_tab$$huaweitest_jar], [1= addrof ptr $__muid_func_undef_tab$$huaweitest_jar, 2= addrof ptr $__muid_func_undef_tab$$huaweitest_jar], [1= addrof ptr $__muid_data_def_tab$$huaweitest_jar, 2= addrof ptr $__muid_data_def_tab$$huaweitest_jar], [1= 0, 2= 0], [1= addrof ptr $__muid_data_undef_tab$$huaweitest_jar, 2= addrof ptr $__muid_data_undef_tab$$huaweitest_jar], [1= addrof ptr $__muid_func_def_muid_tab$$huaweitest_jar, 2= addrof ptr $__muid_func_def_muid_tab$$huaweitest_jar], [1= addrof ptr $__muid_func_undef_muid_tab$$huaweitest_jar, 2= addrof ptr $__muid_func_undef_muid_tab$$huaweitest_jar], [1= addrof ptr $__muid_data_def_muid_tab$$huaweitest_jar, 2= addrof ptr $__muid_data_def_muid_tab$$huaweitest_jar], [1= addrof ptr $__muid_data_undef_muid_tab$$huaweitest_jar, 2= addrof ptr $__muid_data_undef_muid_tab$$huaweitest_jar], [1= addrof ptr $__muid_func_muid_idx_tab$$huaweitest_jar, 2= addrof ptr $__muid_func_muid_idx_tab$$huaweitest_jar], [1= 0, 2= 0]]
var $__reflection_strtab$$huaweitest_jar fstatic <[57] u8> = [0, 76, 104, 117, 97, 119, 101, 105, 116, 101, 115, 116, 59, 0, 109, 97, 105, 110, 0, 40, 91, 76, 106, 97, 118, 97, 47, 108, 97, 110, 103, 47, 83, 116, 114, 105, 110, 103, 59, 41, 86, 0, 48, 33, 48, 0, 60, 105, 110, 105, 116, 62, 0, 40, 41, 86, 0]
var $__compilerVersionNum$$huaweitest_jar <[0] <* void>> = [1, 0]
func &Lhuaweitest_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$Lhuaweitest_3B>>) void {
  funcid 48153
  var %Reg1_R43694 <* <$Lhuaweitest_3B>> localrefvar
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>> localrefvar
  var %__muid_symptr <* void>

  intrinsiccall MCCIncRef (dread ref %_this)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43694)
  dassign %Reg1_R43694 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  regassign ptr %1 (dread ref %Reg1_R57)
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43694))
  intrinsiccall MCCIncRef (dread ref %Reg1_R57)
  intrinsiccall MCCDecRef (regread ptr %1)
  #Call function:Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V
  dassign %__muid_symptr 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[12] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_func_undef_tab$$huaweitest_jar, constval i64 9)))
  icallassigned (dread ptr %__muid_symptr, dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  intrinsiccall MPL_CLEANUP_LOCALREFVARS (dread ref %Reg1_R43694, dread ref %Reg1_R57)
  return ()
}
func &Lhuaweitest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg2_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void {
  funcid 48154
  var %Reg0_R562 <* <$Ljava_2Fio_2FPrintStream_3B>> localrefvar
  var %Reg1_R43 <* <$Ljava_2Flang_2FString_3B>> localrefvar
  var %L_STR_161334 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161335 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161336 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161337 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161338 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161339 <* <$Ljava_2Flang_2FString_3B>>

  #INSTIDX : 0||0000:  getstatic
  intrinsiccall MPL_CLINIT_CHECK (addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B)
  regassign ptr %1 (dread ref %Reg0_R562)
  #Read from: Ljava_2Flang_2FSystem_3B_7Cout
  dassign %Reg0_R562 0 (iread ref <* <* <$Ljava_2Fio_2FPrintStream_3B>>> 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[3] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_data_undef_tab$$huaweitest_jar, constval i64 1))))
  intrinsiccall MCCIncRef (dread ref %Reg0_R562)
  intrinsiccall MCCDecRef (regread ptr %1)
  #INSTIDX : 3||0003:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_907ba9a10d33c74f1c37be3b41544334) { dassign %L_STR_161334 0 }
  intrinsiccall MCCIncRef (dread ptr %L_STR_161334)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43)
  dassign %Reg1_R43 0 (dread ptr %L_STR_161334)
  #INSTIDX : 5||0005:  invokevirtual
  icallassigned (
    iread u64 <* u64> 0 (add ptr (
      iread ptr <* <$__class_meta__>> 6 (iread ref <* <$Ljava_2Flang_2FObject_3B>> 1 (dread ref %Reg0_R562)),
      constval u32 312)),
    dread ref %Reg0_R562,
    dread ref %Reg1_R43) {}
  #INSTIDX : 8||0008:  getstatic
  intrinsiccall MPL_CLINIT_CHECK (addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B)
  regassign ptr %2 (dread ref %Reg0_R562)
  #Read from: Ljava_2Flang_2FSystem_3B_7Cout
  dassign %Reg0_R562 0 (iread ref <* <* <$Ljava_2Fio_2FPrintStream_3B>>> 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[3] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_data_undef_tab$$huaweitest_jar, constval i64 1))))
  intrinsiccall MCCIncRef (dread ref %Reg0_R562)
  intrinsiccall MCCDecRef (regread ptr %2)
  #INSTIDX : 11||000b:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_194a027f7d649763903a7dcda50cf5a0) { dassign %L_STR_161335 0 }
  intrinsiccall MCCIncRef (dread ptr %L_STR_161335)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43)
  dassign %Reg1_R43 0 (dread ptr %L_STR_161335)
  #INSTIDX : 13||000d:  invokevirtual
  icallassigned (
    iread u64 <* u64> 0 (add ptr (
      iread ptr <* <$__class_meta__>> 6 (iread ref <* <$Ljava_2Flang_2FObject_3B>> 1 (dread ref %Reg0_R562)),
      constval u32 312)),
    dread ref %Reg0_R562,
    dread ref %Reg1_R43) {}
  #INSTIDX : 16||0010:  getstatic
  intrinsiccall MPL_CLINIT_CHECK (addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B)
  regassign ptr %3 (dread ref %Reg0_R562)
  #Read from: Ljava_2Flang_2FSystem_3B_7Cout
  dassign %Reg0_R562 0 (iread ref <* <* <$Ljava_2Fio_2FPrintStream_3B>>> 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[3] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_data_undef_tab$$huaweitest_jar, constval i64 1))))
  intrinsiccall MCCIncRef (dread ref %Reg0_R562)
  intrinsiccall MCCDecRef (regread ptr %3)
  #INSTIDX : 19||0013:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_901191d3899ab5fc7be3d97fe762f978) { dassign %L_STR_161336 0 }
  intrinsiccall MCCIncRef (dread ptr %L_STR_161336)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43)
  dassign %Reg1_R43 0 (dread ptr %L_STR_161336)
  #INSTIDX : 21||0015:  invokevirtual
  icallassigned (
    iread u64 <* u64> 0 (add ptr (
      iread ptr <* <$__class_meta__>> 6 (iread ref <* <$Ljava_2Flang_2FObject_3B>> 1 (dread ref %Reg0_R562)),
      constval u32 312)),
    dread ref %Reg0_R562,
    dread ref %Reg1_R43) {}
  #INSTIDX : 24||0018:  getstatic
  intrinsiccall MPL_CLINIT_CHECK (addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B)
  regassign ptr %4 (dread ref %Reg0_R562)
  #Read from: Ljava_2Flang_2FSystem_3B_7Cout
  dassign %Reg0_R562 0 (iread ref <* <* <$Ljava_2Fio_2FPrintStream_3B>>> 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[3] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_data_undef_tab$$huaweitest_jar, constval i64 1))))
  intrinsiccall MCCIncRef (dread ref %Reg0_R562)
  intrinsiccall MCCDecRef (regread ptr %4)
  #INSTIDX : 27||001b:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_4a20a28a5a71422be90661b1072ca29e) { dassign %L_STR_161337 0 }
  intrinsiccall MCCIncRef (dread ptr %L_STR_161337)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43)
  dassign %Reg1_R43 0 (dread ptr %L_STR_161337)
  #INSTIDX : 29||001d:  invokevirtual
  icallassigned (
    iread u64 <* u64> 0 (add ptr (
      iread ptr <* <$__class_meta__>> 6 (iread ref <* <$Ljava_2Flang_2FObject_3B>> 1 (dread ref %Reg0_R562)),
      constval u32 312)),
    dread ref %Reg0_R562,
    dread ref %Reg1_R43) {}
  #INSTIDX : 32||0020:  getstatic
  intrinsiccall MPL_CLINIT_CHECK (addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B)
  regassign ptr %5 (dread ref %Reg0_R562)
  #Read from: Ljava_2Flang_2FSystem_3B_7Cout
  dassign %Reg0_R562 0 (iread ref <* <* <$Ljava_2Fio_2FPrintStream_3B>>> 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[3] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_data_undef_tab$$huaweitest_jar, constval i64 1))))
  intrinsiccall MCCIncRef (dread ref %Reg0_R562)
  intrinsiccall MCCDecRef (regread ptr %5)
  #INSTIDX : 35||0023:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_e4c52869a15445bee0db530b07e9b98a) { dassign %L_STR_161338 0 }
  intrinsiccall MCCIncRef (dread ptr %L_STR_161338)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43)
  dassign %Reg1_R43 0 (dread ptr %L_STR_161338)
  #INSTIDX : 37||0025:  invokevirtual
  icallassigned (
    iread u64 <* u64> 0 (add ptr (
      iread ptr <* <$__class_meta__>> 6 (iread ref <* <$Ljava_2Flang_2FObject_3B>> 1 (dread ref %Reg0_R562)),
      constval u32 312)),
    dread ref %Reg0_R562,
    dread ref %Reg1_R43) {}
  #INSTIDX : 40||0028:  getstatic
  intrinsiccall MPL_CLINIT_CHECK (addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B)
  regassign ptr %6 (dread ref %Reg0_R562)
  #Read from: Ljava_2Flang_2FSystem_3B_7Cout
  dassign %Reg0_R562 0 (iread ref <* <* <$Ljava_2Fio_2FPrintStream_3B>>> 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[3] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_data_undef_tab$$huaweitest_jar, constval i64 1))))
  intrinsiccall MCCIncRef (dread ref %Reg0_R562)
  intrinsiccall MCCDecRef (regread ptr %6)
  #INSTIDX : 43||002b:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_58a8c3b0d70b92b64324906ef9c02c90) { dassign %L_STR_161339 0 }
  intrinsiccall MCCIncRef (dread ptr %L_STR_161339)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43)
  dassign %Reg1_R43 0 (dread ptr %L_STR_161339)
  #INSTIDX : 45||002d:  invokevirtual
  icallassigned (
    iread u64 <* u64> 0 (add ptr (
      iread ptr <* <$__class_meta__>> 6 (iread ref <* <$Ljava_2Flang_2FObject_3B>> 1 (dread ref %Reg0_R562)),
      constval u32 312)),
    dread ref %Reg0_R562,
    dread ref %Reg1_R43) {}
  #INSTIDX : 48||0030:  return
  intrinsiccall MPL_CLEANUP_LOCALREFVARS (dread ref %Reg0_R562, dread ref %Reg1_R43)
  return ()
}

max_stars_repo_name: Bravewtz/OpenArkCompiler

max_issues_repo_name: Bravewtz/OpenArkCompiler

max_forks_repo_name: Bravewtz/OpenArkCompiler


text:
# Kinematik-Berechnung 3. Arm KAS5
# Beschreibung
#
# Modell ohne Zwangsbedingungen. Diese Definitionsdatei erzeugt Dummy-Variablen, damit die Code-Generierung funktioniert.
#
# Die Gelenkvariablen der Schnittgelenke sollen nicht berechnet werden und sind als konstante Ausdrücke in kintmp_s vorgegeben und werden später mit Null ersetzt.
#

# Autor
# Moritz Schappler, schappler@imes.uni-hannover.de, 2018-02
# Institut fuer mechatronische Systeme, Leibniz Universitaet Hannover
# Initialisierung
restart:
kin_constraints_exist := true: # Für Speicherung
;
with(LinearAlgebra):
with(StringTools): # Für Zeitausgabe
;
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
codegen_act := true:
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):

kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
kin_constraints_exist := true:
# Dummy-Einträge für Schnitt-Gelenke
# Die Winkel der Schnittgelenke sind egal, daher keine Berechnung.
# Ersetze alle mit Null
kintmp_qs := kintmp_s*0:
kintmp_qt := kintmp_t*0:
# Speichern der Ergebnisse
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):
# Liste der Parameter speichern
# 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(..,2) )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
kintmp_s:


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

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

$define xc_dimensions_1d

f := (rs, zeta) -> \
 (params_a_a1 + params_a_a2*n_total(rs) + params_a_a3*n_total(rs)^2) * n_total(rs)^params_a_alpha:


text:
# Bluetooth Audio - Media player configuration options

#
# Copyright (c) 2022 Nordic Semiconductor
#
# SPDX-License-Identifier: Apache-2.0
#

config BT_MPL
	bool "Support for media player"
	select BT_CCID
	select EXPERIMENTAL
	help
	  Enables support for media player
	  Note that the provided media player is a sample that only provides a
	  mock-up with no actual media being played.
	  For real media playback, the sample must be extended, hooked up to a
	  real media player or replaced with a real media player.

if BT_MPL

config BT_MPL_MEDIA_PLAYER_NAME
	string "Media Player Name"
	default "Player0"
	help
	  Use this option to set the name of the media player.

config BT_MPL_MEDIA_PLAYER_NAME_MAX
	int "Max length of media player name"
	default 20
	range 1 255
	help
	  Sets the maximum number of bytes (including the null termination) of
	  the name of the media player.

config BT_MPL_ICON_URL
	string "Media player Icon URL"
	default "http://server.some.where/path/icon.png"
	help
	  Use this option to set the URL of the Media Player Icon.

config BT_MPL_ICON_URL_MAX
	int "Max length of media player icon URL"
	default 40
	range 1 255
	help
	  Sets the maximum number of bytes (including the null termination) of
	  the media player icon URL.

config BT_MPL_TRACK_TITLE_MAX
	int "Max length of the title of a track"
	default 40
	range 1 255
	help
	  Sets the maximum number of bytes (including the null termination) of
	  the title of any track in the media player.

config BT_MPL_SEGMENT_NAME_MAX
	int "Max length of the name of a track segment"
	default 25
	range 1 255
	help
	  Sets the the maximum number of bytes (including the null termination)
	  of the name of any track segment in the media player.

config BT_MPL_GROUP_TITLE_MAX
	int "Max length of the title of a group of tracks"
	default BT_MPL_TRACK_TITLE_MAX
	range 1 255
	help
	  Sets the maximum number of bytes (including the null termination) of
	  the title of any group in the media player.

config BT_MPL_OBJECTS
	bool "Support for media player objects"
	depends on BT_OTS
	# TODO: Temporarily depends also on BT_MCS, to avoid issues with the
	# bt_mcs_get_ots() call
	depends on BT_MCS
	help
	  Enables support for objects in the media player
	  Objects are used to give/get more information about e.g. media tracks.
	  Requires the Object Transfer Service

if BT_MPL_OBJECTS

config BT_MPL_MAX_OBJ_SIZE
	int "Total memory size to use for storing the content of objects"
	default 127
	range 0 65536
	help
	  Sets the total memory size (in octets) to use for storing the content of objects.
	  This is used for the total memory pool buffer size from which memory
	  is allocated when sending object content.

config BT_MPL_ICON_BITMAP_SIZE
	int "Media player Icon bitmap object size"
	default 127
	help
	  This option sets the maximum size (in octets) of the icon object.

config BT_MPL_TRACK_MAX_SIZE
	int "Maximum size for a track object"
	default 127
	help
	  This option sets the maximum size (in octets) of a track object.

endif # BT_MPL_OBJECTS

config BT_DEBUG_MPL
	bool "Media player debug"
	help
	  Enables debug logs for the media player

endif # BT_MPL

max_stars_repo_name: ldalek/zephyr

max_issues_repo_name: ldalek/zephyr

max_forks_repo_name: ldalek/zephyr


text:
program sample16;	/* prime numbers */
var furui : array[20000] of boolean;
    i, j : integer;
begin
  i := 2;
  while i < 20000 do begin
    furui[i] := true;
    i := i + 1
  end;
  furui[0] := false;
  furui[1] := false;
  i := 2;
  while i < 20000 do begin
    if furui[i] then begin
      writeln(i, ' is a prime number');
      j := i;
      if i < 16384 then
        while j < 20000 do begin
          furui[j] := false;
          j := j + i
        end
    end;
    i := i + 1
  end
end.

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing


text:
3 + 5 ~

5 - 4 ~

5 * 5 ~
2 / 3 ~

jedi ( 5 < 6 )
{ 2 + 3 }
sith
{ 1 + 6 } ~

x = 5

yell ( x < 6 )
{ 6 + 6 } ~
jedi ( 5 > 6 )
{ 2 + 3 }
sith
{ 1 + 6 } ~
x = 3 ~
yell ( x > 6 )
{ 6 + 6 } ~
quit

max_stars_repo_path: MPL/Example Inputs/errorTest.mpl

max_stars_repo_name: ytugba/My-Otiose-Projects

max_issues_repo_path: MPL/Example Inputs/errorTest.mpl

max_issues_repo_name: ytugba/My-Otiose-Projects

max_forks_repo_path: MPL/Example Inputs/errorTest.mpl

max_forks_repo_name: ytugba/My-Otiose-Projects


text:
 func $floorf64toi32 (
  var %i f64
  ) i32 {
   return (
     floor i32 f64(dread f64 %i))}

 func $floorf64toi64 (
  var %i f64
  ) i64 {
   return (
     floor i64 f64(dread f64 %i))}

 func $floorf32toi32 (
  var %i f32
  ) i32 {
   return (
     floor i32 f32(dread f32 %i))}

 func $floorf32toi64 (
  var %i f32
  ) i64 {
   return (
     floor i64 f32(dread f32 %i))}

# todo float floor
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
// Spawn particles at every ACV_Main that did not have a successful validation.
// Params: area_effect_cloud named ACV_Main tagged with ACV_handleFailedPortals_PARAM
process ACV_handleFailedPortals {
/execute @e[type=area_effect_cloud,name=ACV_Main,tag=ACV_handleFailedPortals_PARAM,score_ACV_COLOR_min=0,score_ACV_COLOR=0] ~ ~ ~ /particle iconcrack ~ ~ ~ 0.1 0.1 0.1 0 20 force @a 22
/execute @e[type=area_effect_cloud,name=ACV_Main,tag=ACV_handleFailedPortals_PARAM,score_ACV_COLOR_min=1,score_ACV_COLOR=1] ~ ~ ~ /particle iconcrack ~ ~ ~ 0.1 0.1 0.1 0 20 force @a 159 1
/execute @e[type=area_effect_cloud,name=ACV_Main,tag=ACV_handleFailedPortals_PARAM,score_ACV_COLOR_min=2,score_ACV_COLOR=2] ~ ~ ~ /particle iconcrack ~ ~ ~ 0.1 0.1 0.1 0 20 force @a 159 11
/execute @e[type=area_effect_cloud,name=ACV_Main,tag=ACV_handleFailedPortals_PARAM,score_ACV_COLOR_min=3,score_ACV_COLOR=3] ~ ~ ~ /particle iconcrack ~ ~ ~ 0.1 0.1 0.1 0 20 force @a 159 14

// TODO: playsound
}


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

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

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

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

$ifdef gga_x_b88_params
params_a_beta  := 0.0042:
params_a_gamma := 6.0:
$endif

b88_f := x -> 1 + params_a_beta/X_FACTOR_C*x^2/(1 + params_a_gamma*params_a_beta*x*arcsinh(x)):

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


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

sys := [
diff(x2(t), t) = alpha*x1(t) - beta*x2(t),
diff(x4(t), t) = (gama*sigma*x2(t)*x4(t) - delta*sigma*x3(t)*x4(t)) / (x3(t)),
diff(x1(t), t) = (-b*c*x1(t) - b*x1(t)*x4(t) + 1) / (c + x4(t)),
diff(x3(t), t) = gama*x2(t) - delta*x3(t),
y(t) = x1(t)
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


text:
`bar/EA` := proc ()
 local i;
 for i from 2 to nargs do
  if args[i] = 1 then return 0 fi;
 od;
 return '`bar/EA`'(args);
end:

# The standard differential in the bar complex is `d/EA`

`d/EA` := (u) -> sort(expand(eval(subs(`bar/EA`=`d_bar/EA`,u))));

`d_bar/EA` := proc()
 local n,i,q,a;
 a := args;
 n := nargs;
 if n = 0 then
  return FAIL;
 elif n = 1 then
  return 0;
 else
  return
   add((-1)^i*`bar/EA`(a[1..i-1],a[i]*a[i+1],a[i+2..n]),i=1..n-1) +
    (-1)^n*`bar/EA`(a[1..n-1]);
 fi;
end:

# The shuffle product in the bar complex is `mu/EA`

`mu/EA` := proc()
 apply_linear_assoc(`mu0/EA`,`bar/EA`())(args);
end:

`mu0/EA` := proc(u,v)
 local n,m,uv,c,L,s,i;
 n := nops(u)-1;
 m := nops(v)-1;
 uv := [op(2..n+1,u),op(2..m+1,v)];
 c := op(1,u) * op(1,v);
 L := `list_elements/inverse_shuffles`(n,m);

 return add(`sgn/shuffles`(n,m)(s) *
             `bar/EA`(c,seq(uv[s[i]],i=1..n+m)),s in L):
end:

`mu_bar/EA` := proc()
 `mu/EA`(op(map(`bar/EA`,[args])));
end:

# We now give a chain homotopy of EA for the case where A is cyclic
# of infinite order generated by a.  This gives a retraction onto
# a subcomplex concentrated in degrees 0 and 1.

`s0/EA` := (a) -> proc(u)
 local n,k,s,i;
 if type(u,specfunc(anything,`bar/EA`)) then
  n := nops(u);
  if n <= 1 then return 0; fi;
  k := degree(op(n,u),a);
  if op(n,u) <> a^k then
   return FAIL;
  fi;
  if k >= 0 then
   return add((-1)^n*`bar/EA`(op(1..n-1,u),a^i,a),i=0..k-1);
  else
   return add((-1)^(n+1)*`bar/EA`(op(1..n-1,u),a^(-i),a),i=1..abs(k));
  fi:
 else
  return FAIL;
 fi;
end:

`s/EA` := (a) -> apply_linear(`s0/EA`(a));

`p0/EA` := (a) -> proc(u)
 local n,i,k;
 if type(u,specfunc(anything,`bar/EA`)) then
  n := nops(u);
  if n <= 1 then
   return u;
  elif n = 2 then
   k := degree(op(2,u),a);
   if op(n,u) <> a^k then
    return FAIL;
   fi;
   if k >= 0 then
    return add(`bar/EA`(a^i*op(1,u),a),i=0..k-1);
   else
    return add(-`bar/EA`(a^(-i)*op(1,u),a),i=1..abs(k));
   fi;
  else
   return 0;
  fi;
 else
  return FAIL;
 fi;

end:

`p/EA` := (a) -> apply_linear(`p0/EA`(a));

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

# Given elements a[i] in A we have elements 1 - a[i] in the group ring.
# The tensor product of these is represented by `bar/BA`(a[1],...,a[n])

`bar/BA` := proc ()
 local i;
 for i to nargs do
  if args[i] = 1 then return 0 fi;
 od;
 return '`bar/BA`'(args);
end:

# The standard differential in the bar complex is `d/BA`

`d/BA` := (u) -> sort(expand(eval(subs(`bar/BA`=`d_bar/BA`,u))));

`d_bar/BA` := proc()
 local n,i,q;
 n := nargs;
 if n = 0 then
  return 0;
 else
  return -`bar/BA`(args[2..n]) +
         (-1)^(n+1)*`bar/BA`(args[1..n-1]) +
         add((-1)^(i+1)*`bar/BA`(args[1..i-1],args[i]*args[i+1],args[i+2..n]),
             i=1..n-1);
 fi;
end:

# The shuffle product in the bar complex is `mu/BA`

`mu/BA` := proc()
 apply_linear_assoc(`mu0/BA`,`bar/BA`())(args);
end:

`mu0/BA` := proc(u,v)
 local n,m,uv,L,s,i;
 n := nops(u);
 m := nops(v);
 uv := [op(u),op(v)];
 L := `list_elements/inverse_shuffles`(n,m);

 return add(`sgn/shuffles`(n,m)(s) * `bar/BA`(seq(uv[s[i]],i=1..n+m)),s in L):
end:

`mu_bar/BA` := proc()
 `mu/BA`(op(map(`bar/BA`,[args])));
end:

`proj_bar/EA` := proc() `bar/BA`(args[2..-1]); end:
`lift_bar/BA` := proc() `bar/EA`(1,args); end:
`proj/EA/BA` := (u) -> eval(subs(`bar/EA` = `proj_bar/EA`,u));
`lift/BA/EA` := (u) -> eval(subs(`bar/BA` = `lift_bar/BA`,u));

######################################################################
# The element u = `tau/BA`(n,a) has du = n[a] - [a^n].
# Thus, if a^n = 1 then u is a 2-cycle modulo n.

`tau/BA` := proc(n,a)
 local k;
 add(-`bar/BA`(a^k,a),k=0..n-1);
end:

`zeta/BA` := proc(n,m,a)
 local i,j;
 - add(add(`bar/BA`(a^(m*i),a^j,a),j=0..m-1),i=0..n-1);
end:

# The element `sigma/BA`(n,a,b) is a chain of dimension 3.
# It is symmetric in a and b, and is a cycle if a^n = b^n = 1.

`sigma/BA` := proc(n,a,b)
 `mu/BA`(`tau/BA`(n, a), `bar/BA`(b)) +
 `mu/BA`(`tau/BA`(n, b), `bar/BA`(a))
end:

# The element u = `gamma/BA`(n,a) is a chain of dimension 2n+1.
# It is a cycle if a^n = 1.

`gamma/BA` := proc(n,k,a)
 local u,v,i;
 u := `bar/BA`(a);
 v := `tau/BA`(n,a);
 for i from 1 to k do
  u := `mu/BA`(u,v);
 od;
 u := expand(u/k!);
 return u;
end:

`xi/BA` := proc(n,m,a)
 local i,j;

 add(add(-`bar/BA`(a^i,a,a^j,a),i=1..m-1),j=1..n*m-1) +
 add(add(-`bar/BA`(a^m,a^(i*m),a^j,a),j=1..m-1),i=1..n-1);
end:

`rho/BA` := proc(n,a,b)
 `tau/BA`(n,a) +
 `tau/BA`(n,b) -
 `tau/BA`(n,a*b) +
 n*`bar/BA`(a, b);
end:

`omega/BA` := proc (n,a,b)
 local i,j;

 -add(add(`bar/BA`(a^i*b^j,b,a), j = i .. n-1), i = 0 .. n-1) +
  add(add(`bar/BA`(a^i*b^j,a,b), j = i+1 .. n-1), i = 0 .. n-2);
end:

######################################################################
# We now introduce the standard DGA for H_*(BZ/n).

# The k'th chain group is freely generated over Z[Z/n] by
# `e/EC`(n)(k,0).  The orbit of this consists of elements `e/EC`(n)(k,i).

`de/EC` := (n::posint) -> proc(k::nonnegint,i::nonnegint)
 local j;

 if k = 0 then
  return 0;
 elif modp(k,2) = 0 then
  return(add(`e/EC`(n)(k-1,j),j=0..n-1));
 else
  return `e/EC`(n)(k-1,modp(i+1,n)) - `e/EC`(n)(k-1,i);
 fi;
end:

`d/EC` := (n) -> (u) -> eval(subs(`e/EC` = `de/EC`,u));

`mu/EC` := (n) -> proc()
 apply_linear_assoc(`mu0/EC`(n),`e/EC`(n)(0,0))(args);
end:

`mu0/EC` := (n) -> proc(a,b)
 local ka,kb,ia,ib,la,lb;
 if type(a,specfunc(anything,`e/EC`(n))) and
    type(b,specfunc(anything,`e/EC`(n))) then
  ka,ia := op(a);
  kb,ib := op(b);
  if modp(ka,2) = 1 and modp(kb,2) = 1 then
   return 0;
  fi;
  la := floor(ka/2);
  lb := floor(kb/2);
  return binomial(la+lb,la) * `e/EC`(n)(ka+kb,modp(ia+ib,n));
 else
  return('`mu/EC`'(n)(args))
 fi;
end:

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

`de/BC` := (n::posint) -> proc(k::nonnegint)
 if k = 0 or modp(k,2) = 1 then
  return 0;
 else
  return n * `e/BC`(n)(k-1);
 fi;
end:

`d/BC` := (n) -> (u) -> eval(subs(`e/BC` = `de/BC`,u));

`mu/BC` := (n) -> proc()
 apply_linear_assoc(`mu0/BC`(n),`e/BC`(n)(0))(args);
end:

`mu0/BC` := (n) -> proc(a,b)
 local ka,kb,la,lb;
 if type(a,specfunc(anything,`e/BC`(n))) and
    type(b,specfunc(anything,`e/BC`(n))) then
  ka := op(a);
  kb := op(b);
  if modp(ka,2) = 1 and modp(kb,2) = 1 then
   return 0;
  fi;
  la := floor(ka/2);
  lb := floor(kb/2);
  return binomial(la+lb,la) * `e/BC`(n)(ka+kb);
 else
  return('`mu/BC`'(n)(args))
 fi;
end:

`e/BC/BA` := (n::posint) -> proc(k::nonnegint)
 option remember;

 if k = 0 then
  return `bar/BA`();
 elif k = 1 then
  return `bar/BA`(a);
 else
  return expand(`mu/BA`(`tau/BA`(n,a),`e/BC/BA`(n)(k-2))/floor(k/2));
 fi;
end:

`phi/BC/BA` := proc(u)
 eval(subs(`e/BC` = `e/BC/BA`,u));
end:

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

check_homology_BA := proc()
 local i,j,k,n,m,a,b,c,t,u,v,w,T,ok,err,exp_rel,diff_rel;

 ok := true;
 for i from 0 to 4 do
  for j from 0 to 4 do
   u := `bar/EA`(seq(a[m],m=1..i+1));
   v := `bar/EA`(seq(a[m],m=i+2..i+j+2));
   err := `mu/EA`(u,v) - (-1)^(i*j) * `mu/EA`(v,u);
   if err <> 0 then
    ok := false;
    break;
   fi;
  od:
  if not(ok) then break; fi;
 od:

 _ASSERT(ok,"`mu/EA` is graded-commutative");

 ok := true;
 for i from 0 to 4 do
  for j from 0 to 4 do
   for k from 0 to 4 do
    u := `bar/EA`(seq(a[m],m=1..i+1));
    v := `bar/EA`(seq(a[m],m=i+2..i+j+2));
    w := `bar/EA`(seq(a[m],m=i+j+3..i+j+k+3));
    err := `mu/EA`(`mu/EA`(u,v),w) - `mu/EA`(u,`mu/EA`(v,w));
    if err <> 0 then
     ok := false;
     break;
    fi;
   od:
   if not(ok) then break; fi;
  od:
  if not(ok) then break; fi;
 od:

 _ASSERT(ok,"`mu/EA` is associative");

 ok := true;
 for i from 0 to 3 do
  for j from 0 to 3 do
   u := `bar/EA`(seq(a[m],m=1..i+1));
   v := `bar/EA`(seq(a[m],m=i+2..i+j+2));
   err := `d/EA`(`mu/EA`(u,v)) - `mu/EA`(`d/EA`(u),v) - (-1)^i*`mu/EA`(u,`d/EA`(v));
   if err <> 0 then
    ok := false;
    break;
   fi;
  od:
  if not(ok) then break; fi;
 od:

 _ASSERT(ok,"`d/EA` and `mu/EA` satisfy the Leibniz rule");

 ok := true;
 for i from 0 to 4 do
  for j from 0 to 4 do
   u := `bar/BA`(seq(a[m],m=1..i));
   v := `bar/BA`(seq(a[m],m=i+1..i+j));
   err := `mu/BA`(u,v) - (-1)^(i*j) * `mu/BA`(v,u);
   if err <> 0 then
    ok := false;
    break;
   fi;
  od:
  if not(ok) then break; fi;
 od:

 _ASSERT(ok,"`mu/BA` is graded-commutative");

 ok := true;
 for i from 0 to 4 do
  for j from 0 to 4 do
   for k from 0 to 4 do
    u := `bar/BA`(seq(a[m],m=1..i));
    v := `bar/BA`(seq(a[m],m=i+1..i+j));
    w := `bar/BA`(seq(a[m],m=i+j+1..i+j+k));
    err := `mu/BA`(`mu/BA`(u,v),w) - `mu/BA`(u,`mu/BA`(v,w));
    if err <> 0 then
     ok := false;
     break;
    fi;
   od:
   if not(ok) then break; fi;
  od:
  if not(ok) then break; fi;
 od:

 _ASSERT(ok,"`mu/BA` is associative");

 ok := true;
 for i from 0 to 4 do
  for j from 0 to 4 do
   u := `bar/BA`(seq(a[m],m=1..i));
   v := `bar/BA`(seq(a[m],m=i+1..i+j));
   err := `d/BA`(`mu/BA`(u,v)) - `mu/BA`(`d/BA`(u),v) - (-1)^i*`mu/BA`(u,`d/BA`(v));
   if err <> 0 then
    ok := false;
    break;
   fi;
  od:
  if not(ok) then break; fi;
 od:

 _ASSERT(ok,"`d/BA` and `mu/BA` satisfy the Leibniz rule");

 n := 5;
 err := `d/BA`(`tau/BA`(n,a)) - n * `bar/BA`(a);
 exp_rel := `bar/BA`(a^n);
 err := err + exp_rel;

 _ASSERT(err = 0,"d(tau_n(a))");

 err :=
  `d/BA`(`mu/BA`(`tau/BA`(n,a),`bar/BA`(b))) - n * `mu_bar/BA`(a,b);
 exp_rel := `mu_bar/BA`(a^n,b);
 err := err + exp_rel;

 _ASSERT(err = 0,"d(tau_n(a) b)");

 _ASSERT(
  {seq(seq(
   `tau/BA`(i*j,a) - `tau/BA`(i,a^j) - i * `tau/BA`(j,a) - `d/BA`(`zeta/BA`(i,j,a)),
  j=1..4),i=1..4)} = {0},
  "tau_{ij}(a)"
 );

 err := `d/BA`(`sigma/BA`(n,a,b));
 exp_rel := `mu_bar/BA`(a^n,b) - `mu_bar/BA`(a,b^n);
 err := err + exp_rel;

 _ASSERT(err = 0,"d(sigma_n(a,b))");

 _ASSERT(`sigma/BA`(n,a,b) = `sigma/BA`(n,b,a),"sigma_n is symmetric");

 err :=
  `sigma/BA`(n,a,b*c) - `sigma/BA`(n,a,b) - `sigma/BA`(n,a,c) +
    `mu/BA`(`bar/BA`(a),`rho/BA`(n,b,c));
 diff_rel := - `mu/BA`(`tau/BA`(n,a),`bar/BA`(b,c));
 exp_rel  := - `mu/BA`(`bar/BA`(a^n),`bar/BA`(b,c));
 err := err + `d/BA`(diff_rel) + exp_rel;

 _ASSERT(err = 0,"First formula for sigma_n(a,bc)");

 err :=
  `sigma/BA`(n,a,b*c)-`sigma/BA`(n,a,b)-`sigma/BA`(n,a,c)+
   binomial(n+1,2)*`mu_bar/BA`(a,b,c);

 diff_rel := -(`mu/BA`(`bar/BA`(a),`omega/BA`(n,b,c))+
  	       `mu/BA`(`tau/BA`(n,a),`bar/BA`(b,c)));
 exp_rel :=
  - `mu/BA`(`bar/BA`(a^n),`bar/BA`(b,c)) -
  add(`bar/BA`(a,b^i,b)-`bar/BA`(a,b^i*c^n,b),i=0..n-1)-
  add(`bar/BA`(b^i,b,a)-`bar/BA`(b^i*c^n,b,a),i=0..n-1)+
  add(`bar/BA`(b^i,a,b)-`bar/BA`(b^i*c^n,a,b),i=0..n-1);

 err := err + `d/BA`(diff_rel) + exp_rel;

 _ASSERT(err = 0,"Second formula for sigma_n(a,bc)");

 for k from 1 to 5 do
  err := `d/BA`(`gamma/BA`(n,k,a));
  T := [[]]:
  for i from 1 to k-1 do
   T := [seq(seq([op(t),a,a^j],j=1..n-1),t in T)]:
  od:
  exp_rel := (-1)^k*add(`mu/BA`(`bar/BA`(op(t),a),`bar/BA`(a^n)),t in T);
  err := err + exp_rel;
  _ASSERT(err = 0,sprintf("d(gamma_{n,%d}(u))",k));
 od:

 _ASSERT(simplify(`sigma/BA`(n,a,a) - 2 * `gamma/BA`(n,1,a)) = 0,
         "sigma_n(a,a) = 2 gamma_{n,1}(a)");

 err :=
  `gamma/BA`(n,1,a*b) - `gamma/BA`(n,1,a) - `gamma/BA`(n,1,b) -
    `sigma/BA`(n,a,b);

 diff_rel :=
  n*`bar/BA`(a,b,a,b) -
  `mu/BA`(`bar/BA`(a,b),`tau/BA`(n,a*b)) +
  `mu/BA`(`omega/BA`(n,a,b),`bar/BA`(a)+`bar/BA`(b));

 exp_rel :=
  add(`bar/BA`(a^i,b,a) - `bar/BA`(a^i*b^n,b,a),i=0..n-1) -
  add(`bar/BA`(a^i,a,b) - `bar/BA`(a^i*b^n,a,b),i=0..n-1) -
  add(`bar/BA`(b,a^i,a) - `bar/BA`(b,a^i*b^n,a),i=0..n-1) -
  add(`bar/BA`(a,a^i,a) - `bar/BA`(a,a^i*b^n,a),i=0..n-1) -
 `mu/BA`(`bar/BA`(a^n*b^n),`bar/BA`(a,b));

 err := err + `d/BA`(diff_rel) + exp_rel;

 _ASSERT(err = 0,"gamma_{n,1}(ab)");

  _ASSERT(
  {seq(seq(
   `gamma/BA`(i,1,a^j) - j * `gamma/BA`(i*j,1,a) +
   `mu/BA`(`bar/BA`(a^(i*j)),`tau/BA`(j,a)) -
   `d/BA`(`xi/BA`(i,j,a)),
  j=1..4),i=1..4)} = {0},
  "gamma_{ij,1}(a)"
 );
end:

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

check_homology_BZ := proc()
 local U,i,j,u,du,sdu,su,dsu,pu,ok,err;

 U := [seq(`bar/EA`(a^j),j=-2..2)]:
 for i from 1 to 4 do
  U := [op(U),
	seq(seq(`bar/EA`(op(u),a^j),j=1..2),u in U),
	seq(seq(`bar/EA`(op(u),a^(-j)),j=1..2),u in U)]:
 od:

 ok := true;

 for u in U do
  du := `d/EA`(u);
  sdu := `s/EA`(a)(du);
  su := `s/EA`(a)(u);
  dsu := `d/EA`(su);
  pu := `p/EA`(a)(u);
  err := sdu + dsu - u + pu;
  if err <> 0 then
   ok := false;
   break;
  fi;
 od:

 _ASSERT(ok,"ds + ds = 1 - p",[u,du,sdu,su,dsu,pu]);
end:

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

check_homology_BC := proc()
 local n,ok,u,v,w,ku,kv,kw,err;

 n := 5;

 ok := true;
 for ku from 0 to 4 do
  for kv from 0 to 4 do
   u := `e/BC`(n)(ku);
   v := `e/BC`(n)(kv);
   err := `mu/BC`(n)(u,v) - (-1)^(ku*kv) * `mu/BC`(n)(v,u);
   if err <> 0 then
    ok := false;
    break;
   fi;
  od:
  if not(ok) then break; fi;
 od:

 _ASSERT(ok,"`mu_BC` is graded-commutative");

 ok := true;
 for ku from 0 to 4 do
  for kv from 0 to 4 do
   for kw from 0 to 4 do
    u := `e/BC`(n)(ku);
    v := `e/BC`(n)(kv);
    w := `e/BC`(n)(kw);
    err := `mu/BC`(n)(`mu/BC`(n)(u,v),w) - `mu/BC`(n)(u,`mu/BC`(n)(v,w));
    if err <> 0 then
     ok := false;
     break;
    fi;
   od:
   if not(ok) then break; fi;
  od:
 od:

 _ASSERT(ok,"`mu_BC` is associative",[u,v,w,err]);

 ok := true;
 for ku from 0 to 4 do
  for kv from 0 to 4 do
   u := `e/BC`(n)(ku);
   v := `e/BC`(n)(kv);
   err := `d/BC`(n)(`mu/BC`(n)(u,v)) - `mu/BC`(n)(`d/BC`(n)(u),v) - (-1)^ku * `mu/BC`(n)(u,`d/BC`(n)(v));
   if err <> 0 then
    ok := false;
    break;
   fi;
  od:
  if not(ok) then break; fi;
 od:

 _ASSERT(ok,"Leibniz rule for BC");
end:


text:
make_tetrahedron_complex := proc()
 global tetrahedron_complex;

 local T,V,E,F,v,f,P2,P3,g,s,x,eqs,sol;

 V := [seq(i,i=1..4)];
 E := [seq(seq([i,j],j=i+1..4),i=1..4)];
 F := [seq(seq(seq([i,j,k],k=j+1..4),j=i+1..4),i=1..4)];

 T := table([]);

 T["vertices"] := V;
 T["edges"] := E;
 T["coedges"] := E;
 T["faces"] := F;
 T["max_simplices"] := F;

 T["edge_index"]   := make_index(T["edges"]):
 T["face_index"]   := make_index(T["faces"]):

 T["simplices"] := [map(v -> [v],V),op(E),op(F)];

 T["embedding_dim"] := 3;
 T["embedding"] := table([
  1 = [ 1, 1, 1], 2 = [ 1,-1,-1], 3 = [-1, 1,-1], 4 = [-1,-1, 1]
 ]);

 T["axial_embedding"] :=
  table([1 = simplex_embedding([1,0,0,0]),
         2 = simplex_embedding([0,1,0,0]),
	 3 = simplex_embedding([0,0,1,0]),
	 4 = simplex_embedding([0,0,0,1])]);

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

 T["normalised_embedding"] :=
  table([seq(i = v[i] *~ (sqrt(3)/3), i = 1..4)]);

 T["dual_embedding"] := table([seq(i = -~ v[i],i = 1..4)]);

 `plot/simplicial_complex`(T);

 T["edge_centres"] :=
  map(e -> (v[e[1]] +~ v[e[2]]) /~ 2,T["edges"]);

 T["face_centres"] :=
  map(f -> (v[f[1]] +~ v[f[2]] +~ v[f[3]]) /~ 3, T["faces"]);

 f := u -> combine(simplify(rationalize(u /~ sqrt(add(u[i]^2,i=1..3)))));
 P2 := [seq(seq(f(v[i] +~ v[j]),j=i+1..4),i=1..3)];
 P3 := [seq(f(v[i]),i=1..4),seq(f(-~ v[i]),i=1..4)];

 T["poles"] := table([2 = P2, 3 = P3]);

 T["planes"] := [
  seq(seq(f(v[i] -~ v[j]),j=i+1..4),i=1..3),
  seq(seq(f(v[j] -~ v[i]),j=i+1..4),i=1..3)
 ];

 T["pole_plots"] := table([
  2 = map(u -> line(1.2 *~ u, -1.2 *~ u,color=green,thickness=4),T["poles"][2]),
  3 = map(u -> line(1.2 *~ u, -1.2 *~ u,color=red  ,thickness=4),T["poles"][3])
 ]);

 T["all_poles_plot"] := display(
  op(T["pole_plots"][2]),op(T["pole_plots"][3]),
  scaling=constrained,axes=none
 );

 T["symmetry_group"] := combinat[permute](4);
 T["rotation_group"] :=
  select(s -> mul(mul(s[j] - s[i],j=i+1..4),i=1..3) > 0, T["symmetry_group"]);

 T["symmetry_matrix"] := table():

 for s in T["symmetry_group"] do
  g := Matrix(3,3,[seq(x[i],i=1..9)]);
  eqs := map(op,[seq(convert(g . Vector(v[i]) - Vector(v[s[i]]),list),i = 1..4)]);
  sol := solve(eqs);
  T["symmetry_matrix"][s] := simplify(rationalize(subs(sol,convert(g,listlist))));
 od:

 tetrahedron_complex := eval(T);
 return eval(T);
end:

make_tetrahedron_complex():


text:
`is_element/total_preord` := (A::set) -> proc(R)
 local AA,a,b,ab,U;

 global reason;

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

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

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

 return true;
end;

`is_equal/total_preord` := (A::set) -> (R1,R2) -> `is_equal/preord`(A)(R1,R2);
`is_leq/total_preord`   := (A::set) -> (R1,R2) -> `is_leq/preord`(A)(R1,R2);

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

`build/total_preord` := (A::set) -> proc(pi)
 local m,r,i,a,b,R;

 m := nops(pi);
 r := table();
 for i from 1 to m do
  for a in pi[i] do
   r[a] := i;
  od;
 od;

 R := NULL;
 for a in A do
  for b in A do
   if r[a] <= r[b] then R := R,[a,b]; fi;
  od;
 od;
 R := {R};

 return R;
end;

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

`random_element/total_preord` := (A::set) -> proc()
 local pi,pi0,pi1,m,r,i,a,b,R;

 pi := `random_element/partitions`(A)();
 pi0 := map(u -> sort([op(u)]),pi);
 pi0 := sort([op(pi0)]);
 pi1 := combinat[randperm](pi0);
 return `build/total_preord`(A)(pi1);
end;

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

`list_elements/total_preord` := proc(A::set)
 local PI,pi,pi0,pi1,L;

 PI := `list_elements/partitions`(A);
 L := NULL;
 for pi in PI do
  pi0 := map(u -> sort([op(u)]),pi);
  pi0 := sort([op(pi0)]);
  for pi1 in combinat[permute](pi0) do
   L := L,`build/total_preord`(A)(pi1);
  od;
 od;
 L := [L];
 return L;
end;

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

`count_elements/total_preord` := proc(A::set)
 local n,k;
 n := nops(A);
 add(k! * Stirling2(n,k),k=1..n);
end;


text:
program assigninteger;
var n : integer;
begin
    n := 1
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing


text:
program sample31p;
var a : integer;
procedure p(a : char);
begin
	writeln('proc of p');
	a := 'a'
end;
var b : char;
procedure q(b:integer);
  var a : boolean;
      q : integer;
begin
	writeln('proc of q');
	readln(q);
	a := b = q;
	if a then writeln('true') else writeln('false')
end;
var c : integer;
begin
	a := 1;		b := 'b';
	call p(b);
	call q(a);call q(2*a+1)
end.


text:
program sample11p;
var n : integer;
procedure kazuyomikomi;
begin
	writeln('input the number of data');
	readln(n)
end;
var sum : integer;
procedure wakakidasi;
begin
	writeln('Sum of data = ', sum)
end;
var data : integer;
procedure goukei;
begin
	sum := 0;
	while n > 0 do begin
		readln(data);
		sum := sum + data;
		n := n - 1
	end
end;
begin
	call kazuyomikomi;
	call goukei;
	call wakakidasi
end.


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

func $fact (
 var %n i32, var %m f32) i32 {
 regassign i32 %987 (constval i32 '0')
 if (ne i32 i32 (regread i32 %987, constval i32 10)) {
    eval (regread i32 %987)
    call &foo(
      sub i32 (regread i32 %1, constval i32 1))
    return (regread i32 %%retval)}
   else {
     regassign i32 %987 (constval i32 111)
     return (constval i32 '1') }
 }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl