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

with(StringTools):
scanners := map(l->[parse]~(l[2..-1]),
     map(Split,Split(SubstituteAll(input, "\n\n", "X"),"X"),"\n")):

truelocs := table();
truelocs[1] := scanners[1]; # base 0

sunknown := [seq(2..nops(scanners))];
knownscans := [1=[0,0,0]];

dist := (v,u) -> add((v[i]-u[i])^2,i=1..3):

# used to generate equations
M := Matrix(3,3,(i,j)->a[i,j]):
A := <a[0,1], a[0,2], a[0,3]>:

for k1 to nops(scanners)-1 do

    # find a scanner that overlaps with a known scanner
    for j1 from 1 to nops(knownscans) do
        sn1 := lhs(knownscans[j1]);
        scan0 := truelocs[sn1];

        for i1 from 1 to nops(sunknown) do
            # compute all pairwise distances
                sn2 := sunknown[i1];
                scan1 := scanners[sn2];
                n0 := nops(scan0); n1 := nops(scan1);
                # building pairs"
                M0 := Matrix([seq([seq(dist(scan0[i],scan0[j]),j=1..n0)],i=1..n0)]):
                M1 := Matrix([seq([seq(dist(scan1[i],scan1[j]),j=1..n1)],i=1..n1)]):

            ints01 := {entries(M0,'nolist')} intersect {entries(M1,'nolist')} minus {0};

            if nops(ints01) < 50 then
#                printf("scanner %d can not overlap  scanner %d -- only %d common pairs\n", sn1, sn2, nops(ints01) );
                next;
            end if; # at least 66 needed for 12 in common

                # looking up pair from their common distances
            pairs0 := NULL: pairs1 := NULL:
            for i from 1 to nops(ints01) do
                member(ints01[i], M0, 'loc0'); pairs0 := pairs0, [loc0];
                member(ints01[i], M1, 'loc1'); pairs1 := pairs1, [loc1];nu
            end do:
            pairs0 := [pairs0];
            pairs1 := [pairs1];

            # these are the overlapping sets of beacons
            o0 := [seq](ifelse(numboccur(pairs0,i)=11,i,NULL),i=1..n0);
            o1 := [seq](ifelse(numboccur(pairs1,i)=11,i,NULL),i=1..n1);
            if nops(o0) <> 12 then
                #printf("found only %d=%d overlaps between %d and %d skipping\n", nops(o0), nops(o1), sn1, sn2);
                next;
            end if;

            # determine the pairings
            pot01matches := table(sparse={o1[]});

            for i from 1 to nops(pairs0) do
                if pairs0[i][1] in o0 and pairs0[i][2] in o0 then # check o2 too?
                pot01matches[pairs0[i][1]] := pot01matches[pairs0[i][1]] intersect {pairs1[i][]};
                pot01matches[pairs0[i][2]] := pot01matches[pairs0[i][2]] intersect {pairs1[i][]};
                end if;
            end do;

            if {} in {entries(pot01matches,'nolist')} then
                printf("found a contradiction in the potential pairings\n");
                next;
            end if;

            matches01 := map(p->lhs(p)=rhs(p)[], [entries(pot01matches,'pairs')]);

            # use symbolic computation to build the equations to solve for the rotation
            # and translation
            eqs := {seq(seq((M . Vector(scan1[rhs(matches01[k])]) + A)[j]
                = Vector(scan0[lhs(matches01[k])])[j], j=1..3),k=1..4)};
            sol := solve( eqs );

            Mo := eval(M, sol);
            a0 := eval(<a[0,1], a[0,2], a[0,3]>, sol);

            scan1 := [seq(convert(Mo . Vector(scan1[i]) + a0,list), i=1..nops(scan1))];
            knownscans := [knownscans[], sn2=convert(a0,list)];
            sunknown := subsop(i1=NULL, sunknown); # remove this scan from the unknowns
            truelocs[sn2] := scan1;
            break;

        end do;
    end do;
end do;

ASSERT(nops(knownscans) = nops(scanners));
answer1 := nops( map(op, {entries}(truelocs,'nolist') ) );

# compute the Manhatten distance between scanners
mdist := (u,v) -> add(abs(u[i]-v[i]),i=1..3):
scannerlocs := map(rhs, knownscans);
answer2 := max( seq(seq(mdist(scannerlocs[i],scannerlocs[j]),j=i+1..nops(scannerlocs)),i=1..nops(scannerlocs)));


text:
NewSLO:-Profile(
  proc()
    writeto(`if`(StringTools[Has](kernelopts(system), "Windows"), "NUL", "/dev/null"));
    read("NewSLOTests.mpl");
    writeto(terminal);
  end proc,
  {_prof_name="NewSLO"});

ext: mpl

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


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

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

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

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

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

$define gga_c_pbe_params
$include "gga_c_pbe.mpl"

(* the prefactor of t was chosen to get the right K_FACTOR_C in mgga_series_w *)
m08_f := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  + mgga_series_w(params_a_m08_a, 12, 2^(2/3)*t_total(z, ts0, ts1))
    * f_pw(rs, z)
  + mgga_series_w(params_a_m08_b, 12, 2^(2/3)*t_total(z, ts0, ts1))
    * (f_pbe(rs, z, xt, xs0, xs1) - f_pw(rs, z)):

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


text:
# Elliptic curves in Edwards form
# (This is rudimentary; there is a LaTeX file with much more material.)

`rels/edwards` := (a) -> proc(u)
 [-u[1]^2 + u[2]^2 - u[3]^2 + u[4]^2,a * u[2]*u[4] - u[1]*u[3]];
end:

`is_element/edwards_curve` := (a) -> proc(u)
 if not(`is_element/CP`(3)(u)) then return false; fi;
 return evalb(simplify(`rels/edwards`(a)(u)) = [0,0]);
end:


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

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

(* type: lda_exc *)

(* eq 26 *)
a1 := -0.0357:
a2 :=  0.0562:
b1 := -0.0311:
b2 :=  2.39:

(* eq 25 *)
f := (rs, zeta) -> a1/(1 + a2*rs/RS_FACTOR)
  + b1*log((rs/RS_FACTOR + b2)/(rs/RS_FACTOR)):


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

$include "attenuation.mpl"
$define gga_x_b88_params
$include "gga_x_b88.mpl"

(* Eq. 5 - Note that their K_s = 2*X_FACTOR_C*b88_f
   Note that there is a misprint sqrt(mu) -> sqrt(pi)
*)
lcgau_arg := (rs, z, xs) ->
  sqrt(2*X_FACTOR_C*b88_f(xs)/Pi)/(6*n_spin(rs, z)^(1/3)):

(* Eq. 4 + Eq. 8 *)
lcgau_f := (rs, z, xs) ->  b88_f(xs) * (
  + attenuation_erf(p_a_hyb_omega_0_*lcgau_arg(rs, z, xs))
  +  p_a_hyb_coeff_2_*attenuation_gau(p_a_hyb_omega_2_*lcgau_arg(rs, z, xs))
  +  p_a_hyb_coeff_3_*attenuation_gau(p_a_hyb_omega_3_*lcgau_arg(rs, z, xs))
  ):

f := (rs, zeta, xt, xs0, xs1) -> gga_exchange_nsp(lcgau_f, rs, zeta, xs0, xs1):


text:
# https://www.reddit.com/r/adventofcode/comments/ra88up/2021_day_6_part_4_day_googolplex/
# Calculate the population at the googolplex generation mod nextprime(10^8)
input := "4,1,7,7,4,7,6,2,5,4,3,1,4,7,2,4,5,2,2,1,3,7,4,5,1,3,3,5,5,7,6,3,3,3,7,7,5,4,6,3,1,7,6,1,3,5,1,2,6,6,5,5,4,3,2,6,5,3,7,5,4,2,1,3,6,2,7,2,2,6,5,6,7,6,3,3,1,1,1,3,7,3,3,5,4,7,2,1,4,4,1,2,5,5,4,3,4,4,7,4,2,1,2,2,4":
tally := table(sparse=0,Statistics:-Tally(StringTools:-Split(input,","))):
lanternfish := Vector([seq(tally[cat("",i)],i=0..8)]):

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

minpoly := LinearAlgebra:-MinimalPolynomial(U,x);
# x^9 - x^2 - 1

q := 100000007:

# look for a power of x that is 1 modulo the minpoly in F_q[x]
for i from 1 to 100 do
   qq := q^i-1;
   # x^qq ( mod minpoly ) in F_q[x]
   if Powmod(x, qq, minpoly, x) mod q = 1 then
        break;
   end if;
end do:
i, qq;
# i=8, qq=10000005600001372000192080016807000941192032941720658834405764800

# qq is the period of powers of U mod q

# calculate googolplex mod qq
n := Power(10, 10^100) mod qq;
# 7308711993169798294967234678458568672383628022002873653735230400
# so U^googolplex = U^n mod q

# I like computing powers of U using the minpoly
g := Powmod(x, n, minpoly, x) mod q;
# 46714893*x^8+75968618*x^7+86556894*x^6+35145482*x^5+7013179*x^4
#  +30172454*x^3+81311138*x^2+30235098*x+46689635

# g(U) = U^googolplex mod q

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


text:
program sample34;
var num : integer;  ch : char;

begin
	writeln('input the number of data and data character (readln)');
	readln(num);  readln(ch);
	writeln;
	while num > 0 do begin
		write(ch);
		num := num - 1
	end;
	writeln;
	writeln('input the number of data and data character (read)');
	read(num);  read(ch);
	writeln;
	while num > 0 do begin
		write(ch);
		num := num - 1
	end;
	writeln
end.

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing


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

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

(* type: gga_exc *)

$define xc_dimensions_2d

_2d_pbe_kappa := 0.4604:
_2d_pbe_mu    := 0.354546875:

_2d_pbe_f0 := s -> 1 + _2d_pbe_kappa*(1 - _2d_pbe_kappa/(_2d_pbe_kappa + _2d_pbe_mu*s^2)):
_2d_pbe_f  := x -> _2d_pbe_f0(X2S_2D*x):

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


text:
export Improve := module ()
    uses Domain_Type;
    export Simplifiers := table();
    export ModuleApply := proc(dom :: HDomain, $)::HDomain_mb;
        local es := sort( [indices(Domain:-Improve:-Simplifiers,nolist) ]
                          , key=(si->Simplifiers[si]:-SimplOrder))
             , bnd, sh;
        bnd, sh := op([1..2], dom);
        sh := foldr(((f,q)->proc() cmp_simp_sh(f,q,args) end proc)
                   ,(proc(_,b,$) b end proc), op(es))(bnd, sh);
        DOMAIN(bnd, sh);
    end proc;

    local ModuleLoad := proc($)
      LoadSimplifiers();
      LoadSimplifiers := (proc() end proc);
    end proc;#ModuleLoad

    local LoadSimplifiers := proc($)
      local lib_path, drs, improve_path, simpls_paths, simpl_path, names0, names1, new_names, nm;
      unprotect(Simplifiers);
      lib_path := LibraryTools:-FindLibrary(Hakaru);
      if lib_path=NULL then error "Hakaru library not found"; end if;
      lib_path := FileTools:-ParentDirectory(lib_path);
      drs := kernelopts(dirsep);
      improve_path := `cat`(lib_path,drs,"Domain",drs,"Improve");

      simpls_paths := FileTools:-ListDirectory(improve_path);
      for simpl_path in simpls_paths do
        names0 := {anames(user)};
        read(`cat`(improve_path,drs,simpl_path));
        names1 := {anames(user)};
        new_names := names1 minus names0;
        for nm in new_names minus {entries(Simplifiers,nolist)} do
          if nm :: `module`(SimplName,SimplOrder) then
            if assigned(Simplifiers[nm:-SimplName])
            and not Simplifiers[nm:-SimplName] = nm then
              WARNING("overwriting old simplifier! (%1)", Simplifiers[nm:-SimplName]):
            end if;
            Simplifiers[nm:-SimplName] := eval(nm);
          else
            WARNING("not recognized as a simplifier: %1", nm);
          end if;
        end do;
      end do;
      protect(Simplifiers);
    end proc;#LoadSimplifiers

    local combine_errs := proc(err::DomNoSol, mb, $)
        if mb :: DomNoSol then
            DNoSol(op(err),op(mb));
        else
            mb
        end if;
    end proc;

    # compose two simplifiers, combining errors if both fail
    local cmp_simp_sh := proc(simp0, simp1, bnd, sh :: {DomShape,DomNoSol}, $)::{DomShape,DomNoSol};
      local res, sh1; sh1 := Simplifiers[simp0](bnd, sh);
      if not sh1 :: DomNoSol then
          simp1(bnd, sh1);
      else
          combine_errs( sh1, simp1(bnd, sh) );
      end if;
    end proc;
end module;


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

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

(* type: lda_exc *)

$define xc_dimensions_2d

ax := -4/3*sqrt(2)/Pi:

f := (rs, z) -> ax*f_zeta_2d(z)/rs:


text:
program sample28p;
procedure p;begin writeln('Hello!') end;
procedure q;begin writeln('Everyone!') end;
begin call p; call q end.


text:
lprint("ConRec(P,Q,n,k,x,t,Ls,shift,Pindex,Qindex), ConRecq(P,Q,n,k,x,t,Ls,shift,Pindex,Qindex,q)");

## Compute the recurrence relations of the connection coefficients
## Input:
##   P,Q: two hypergeometric terms in k, degree index by n,
##        variable in x=x(t)
##
##   shift: belonds to {`shiftP`,`shiftQ`},
##            recurrence relation involving shifts of the index of P or Q.
##
##   Pindex, Qindex: the indices of P and Q
##
## Output:
##   ls: list of the coefficients of the recurrence relations
##   if length is 3, then
##   ls[1]c[m+1,n]+ls[2]c[m,n]+ls[3]c[m-1,n] = 0.
##   if length is 5, then
##   ls[1]c[m+2,n]+ls[2]c[m+1,n]+ls[3]c[m,n]+ls[4]c[m-1,n]+ls[5]c[m-2,n]=0.
##
## Example
##   P:=poch([la,2 la],n)/2^n*hyperterm([-n,n+2 la],[la+1/2],(1-s*x)/2,k);
##   Q:=2^(-n)*hyperterm([-n,n],[1/2],(1-x)/2,k);
##   L:=proc(f,x) diff(f,x); end proc:
##   ConRec(P,Q,n,k,x,x,L,`shiftP`,n,m);

ConRec:=proc(P,Q,n,k,x,t,Ls,shift,varm,varn)
local L,Lp,lsp,lsq,lspp,lsqp,lsps,lsqs,re;

  ## Ls is either an operator L,
  ## or it is the list of two operators [L, Lp]

  if type(Ls,list) then
    L:=Ls[1]; Lp:=Ls[2];
  else
    L:=Ls; Lp:=0;
  fi;

  ## the 3-term recurrences of P and Q
  lsp:=tr_coe(P,x,n,k,t);
  lsq:=tr_coe(Q,x,n,k,t);

  ## the 3-term recurrences of LP and LQ
  lspp:=tr_coe(L(P,t),x,n,k,t);
  lsqp:=tr_coe(L(Q,t),x,n,k,t);

  if shift=`shiftP` then
    if Lp<>0 then
      re:=Ext_Zeil([Lp(P,t),subs(n=n+1,P),P,subs(n=n-1,P)],k,[t]);
      lsps:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

      re:=Ext_Zeil([Lp(Q,t),subs(n=n+1,Q),Q,subs(n=n-1,Q)],k,[t]);
      lsqs:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

      return(recgenM(lsp,lsq,lspp,lsqp,n,varm,varn,lsps,lsqs));
    else
      return(recgenM(lsp,lsq,lspp,lsqp,n,varm,varn));
    fi;
  else
    if Lp<>0 then
      re:=Ext_Zeil([Lp(P,t),subs(n=n+1,P),P,subs(n=n-1,P)],k,[t]);
      lsps:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

      re:=Ext_Zeil([Lp(Q,t),subs(n=n+1,Q),Q,subs(n=n-1,Q)],k,[t]);
      lsqs:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

      return(recgenN(lsp,lsq,lspp,lsqp,n,varm,varn,lsps,lsqs));
    else
      return(recgenN(lsp,lsq,lspp,lsqp,n,varm,varn));
    fi;
  fi;

end proc:

## Compute the recurrence relations of the connection coefficients / q-case
## Input:
##   P,Q: two hypergeometric terms in k, degree index by n,
##        variable in x=x(t)
##
##   shift: belonds to {`shiftP`,`shiftQ`},
##            recurrence relation involving shifts of the index of P or Q.
##
##   Pindex, Qindex: the indices of P and Q
##
## Output:
##   ls: list of the coefficients of the recurrence relations
##   if length is 3, then
##   ls[1]c[m+1,n]+ls[2]c[m,n]+ls[3]c[m-1,n] = 0.
##   if length is 5, then
##   ls[1]c[m+2,n]+ls[2]c[m+1,n]+ls[3]c[m,n]+ls[4]c[m-1,n]+ls[5]c[m-2,n]=0.
##
## Example
##   P:=poch([la,2 la],n)/2^n*hyperterm([-n,n+2 la],[la+1/2],(1-s*x)/2,k);
##   Q:=2^(-n)*hyperterm([-n,n],[1/2],(1-x)/2,k);
##   L:=proc(f,x) diff(f,x); end proc:
##   ConRec(P,Q,n,k,x,x,L,`shiftP`,n,m);

ConRecq:=proc(P,Q,qn,k,x,t,Ls,shift,varm,varn,q)
local L,Lp,lsp,lsq,lspp,lsqp,lsps,lsqs,re,R;

  ## Ls is either an operator L,
  ## or it is the list of two operators [L, Lp]

  if type(Ls,list) then
    L:=Ls[1]; Lp:=Ls[2];
  else
    L:=Ls; Lp:=0;
  fi;

  ## the 3-term recurrences of P and Q
  re:=qExt_Zeil([x*P,subs(qn=qn*q,P),P,subs(qn=qn/q,P)],k,q,[t]);
  lsp:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

  re:=qExt_Zeil([x*Q,subs(qn=qn*q,Q),Q,subs(qn=qn/q,Q)],k,q,[t]);
  lsq:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

  ## the 3-term recurrences of LP and LQ
  R:=qhyper_simp(L(P,t));
  re:=qExt_Zeil([x*R,subs(qn=qn*q,R),R,subs(qn=qn/q,R)],k,q,[t]);
  lspp:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

  R:=qhyper_simp(L(Q,t));
  re:=qExt_Zeil([x*R,subs(qn=qn*q,R),R,subs(qn=qn/q,R)],k,q,[t]);
  lsqp:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

  if shift=`shiftP` then
    if Lp<>0 then
      re:=qExt_Zeil([qhyper_simp(Lp(P,t)),subs(qn=qn*q,P),P,subs(qn=qn/q,P)],k,q,[t]);
      lsps:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

      re:=qExt_Zeil([qhyper_simp(Lp(Q,t)),subs(qn=qn*q,Q),Q,subs(qn=qn/q,Q)],k,q,[t]);
      lsqs:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

      return(recgenMq(lsp,lsq,lspp,lsqp,qn,varm,varn,lsps,lsqs,q));
    else
      return(recgenMq(lsp,lsq,lspp,lsqp,qn,varm,varn,q));
    fi;
  else
    if Lp<>0 then
      re:=qExt_Zeil([qhyper_simp(Lp(P,t)),subs(qn=qn*q,P),P,subs(qn=qn/q,P)],k,q,[t]);
      lsps:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

      re:=qExt_Zeil([qhyper_simp(Lp(Q,t)),subs(qn=qn*q,Q),Q,subs(qn=qn/q,Q)],k,q,[t]);
      lsqs:=[-re[2]/re[1],-re[3]/re[1],-re[4]/re[1]];

      return(recgenNq(lsp,lsq,lspp,lsqp,qn,varm,varn,lsps,lsqs,q));
    else
      return(recgenNq(lsp,lsq,lspp,lsqp,qn,varm,varn,q));
    fi;
  fi;

end proc:

## compute the coefficients in the three-term relations
## The input:
## hypergeometric term 'P' on variable 'la' which is a polynomial in 'x'
## the order of the polynomial is 'n', the summation index is 'k'
## returen the recurrence
## la*P_n(la) = c[1]*P_{n+1}(la) + c[2]*P_n(la) + c[3]*P_{n-1}(la)
##
tr_coe:=proc(P,la,n,k,x)
local re;

     re:=Ext_Zeil([la*P, subs(n=n+1,P), P, subs(n=n-1, P)], k, [x]);
     return([-re[2]/re[1], -re[3]/re[1], -re[4]/re[1]]);
end proc:

recgenM:=proc(recp,recq,recpp,recqp,var,m,n)

  local c11,c12,b1,c21,c22,b2,res,u,v,
        delta,xi,eta,deltap,xip,etap,detc,
  recps,recqs,c31,c32,b3;

  c11:=subs(var=n-1,recq[1]): # u_{n-1}
  c12:=subs(var=n+1,recq[3]): # w_{n+1}
  b1:=[subs(var=m,recp[1]),
        subs(var=m,recp[2])-subs(var=n,recq[2]),
        subs(var=m,recp[3])]; # right hand side of (2.3)

  c21:=subs(var=n-1,recqp[1]): # u'_{n-1}
  c22:=subs(var=n+1,recqp[3]): # w'_{n+1}
  b2:=[subs(var=m,recpp[1]),
        subs(var=m,recpp[2])-subs(var=n,recqp[2]),
        subs(var=m,recpp[3])]; # right hand side of (2.4)

  detc:=factor(c11*c22-c12*c21);
  if detc=0 then
    return([seq(c11*b2[k]-c21*b1[k],k=1..3)]);
  else
    u:=map(factor, [seq((c22*b1[k]-c12*b2[k])/detc,k=1..3)]);  # (2.5)
    v:=map(factor, [seq((-c21*b1[k]+c11*b2[k])/detc,k=1..3)]); # (2.6)

    if nargs=7 then  # no sigma(x)

      delta:=factor(u[1]): xi:=factor(u[2]):
      eta:=factor(u[3]):

      deltap:=factor(v[1]): xip:=factor(v[2]):
      etap:=factor(v[3]):

      return(map(factor,[subs(n=n+1,delta)*subs(m=m+1,deltap),
          subs(n=n+1,delta)*subs(m=m+1,xip)+subs(n=n+1,xi)*deltap,
          subs(n=n+1,delta)*subs(m=m+1,etap)+subs(n=n+1,xi)*xip
             +subs(n=n+1,eta)*subs(m=m-1,deltap)-1,
          subs(n=n+1,xi)*etap+subs(n=n+1,eta)*subs(m=m-1,xip),
          subs(n=n+1,eta)*subs(m=m-1,etap)]));

    elif nargs=9 then # with sigma(x)

      recps:=args[8]: recqs:=args[9]:
      c31:=subs(var=n-1,recqs[1]): # u''_{n-1}
      c32:=subs(var=n+1,recqs[3]): # w''_{n+1}
      b3:=[subs(var=m,recps[1]),
             subs(var=m,recps[2])-subs(var=n,recqs[2]),
             subs(var=m,recps[3])]; # right hand side of (2.8)

      return(map(factor,[seq(c31*u[k]+c32*v[k]-b3[k],k=1..3)]));
    fi;
  fi;

end proc:

recgenN:=proc(recp,recq,recpp,recqp,var,m,n)

  local c11,c12,b1,c21,c22,b2,res,u,v,
        delta,xi,eta,deltap,xip,etap,detc,
  recps,recqs,c31,c32,b3;

  c11:=subs(var=m,recp[1]):
  c12:=subs(var=m,recp[3]):
  b1:=[subs(var=n-1,recq[1]),
        subs(var=n,recq[2])-subs(var=m,recp[2]),
        subs(var=n+1,recq[3])];

  c21:=subs(var=m,recpp[1]):
  c22:=subs(var=m,recpp[3]):
  b2:=[subs(var=n-1,recqp[1]),
        subs(var=n,recqp[2])-subs(var=m,recpp[2]),
        subs(var=n+1,recqp[3])];

  detc:=factor(c11*c22-c12*c21);
  if detc=0 then
    return([seq(c11*b2[k]-c21*b1[k],k=1..3)]);
  else
    u:=[seq((c22*b1[k]-c12*b2[k])/detc,k=1..3)];
    v:=[seq((-c21*b1[k]+c11*b2[k])/detc,k=1..3)];

    if nargs=7 then

      delta:=factor(u[1]): xi:=factor(u[2]):
      eta:=factor(u[3]):

      deltap:=factor(v[1]): xip:=factor(v[2]):
      etap:=factor(v[3]):

      return(map(factor,[subs(m=m-1,eta)*subs(n=n+1,etap),
          subs(m=m-1,xi)*etap+subs(m=m-1,eta)*subs(n=n+1,xip),
          subs(m=m-1,delta)*subs(n=n-1,etap)+subs(m=m-1,xi)*xip
             +subs(m=m-1,eta)*subs(n=n+1,deltap)-1,
          subs(m=m-1,delta)*subs(n=n-1,xip)+subs(m=m-1,xi)*deltap,
          subs(m=m-1,delta)*subs(n=n-1,deltap)]));

    elif nargs=9 then

      recps:=args[8]: recqs:=args[9]:

      c31:=subs(var=m,recps[1]):
      c32:=subs(var=m,recps[3]):
      b3:=[subs(var=n-1,recqs[1]),
             subs(var=n,recqs[2])-subs(var=m,recps[2]),
             subs(var=n+1,recqs[3])];

      return(map(factor,[seq(c31*u[4-k]+c32*v[4-k]-b3[4-k],k=1..3)]));
    fi;
  fi;

end proc:

recgenMq:=proc(recp,recq,recpp,recqp,var,qm,qn,q)

  local c11,c12,b1,c21,c22,b2,res,u,v,
        delta,xi,eta,deltap,xip,etap,detc,
  recps,recqs,c31,c32,b3;

  c11:=subs(var=qn/q,recq[1]): # u_{n-1}
  c12:=subs(var=qn*q,recq[3]): # w_{n+1}
  b1:=[subs(var=qm,recp[1]),
        subs(var=qm,recp[2])-subs(var=qn,recq[2]),
        subs(var=qm,recp[3])]; # right hand side of (2.3)

  c21:=subs(var=qn/q,recqp[1]): # u'_{n-1}
  c22:=subs(var=qn*q,recqp[3]): # w'_{n+1}
  b2:=[subs(var=qm,recpp[1]),
        subs(var=qm,recpp[2])-subs(var=qn,recqp[2]),
        subs(var=qm,recpp[3])]; # right hand side of (2.4)

  detc:=factor(c11*c22-c12*c21);
  if detc=0 then
    print("D=0");
    return([seq(c11*b2[k]-c21*b1[k],k=1..3)]);
  else
    u:=map(factor, [seq((c22*b1[k]-c12*b2[k])/detc,k=1..3)]);  # (2.5)
    v:=map(factor, [seq((-c21*b1[k]+c11*b2[k])/detc,k=1..3)]); # (2.6)

    if nargs=8 then  # no sigma(x)

      delta:=factor(u[1]): xi:=factor(u[2]):
      eta:=factor(u[3]):

      deltap:=factor(v[1]): xip:=factor(v[2]):
      etap:=factor(v[3]):

      return(map(factor,[subs(qn=qn*q,delta)*subs(qm=qm*q,deltap),
          subs(qn=qn*q,delta)*subs(qm=qm*q,xip)+subs(qn=qn*q,xi)*deltap,
          subs(qn=qn*q,delta)*subs(qm=qm*q,etap)+subs(qn=qn*q,xi)*xip
             +subs(qn=qn*q,eta)*subs(qm=qm/q,deltap)-1,
          subs(qn=qn*q,xi)*etap+subs(qn=qn*q,eta)*subs(qm=qm/q,xip),
          subs(qn=qn*q,eta)*subs(qm=qm/q,etap)]));

    elif nargs=10 then # with sigma(x)

      recps:=args[9]: recqs:=args[10]:
      c31:=subs(var=qn/q,recqs[1]): # u''_{n-1}
      c32:=subs(var=qn*q,recqs[3]): # w''_{n+1}
      b3:=[subs(var=qm,recps[1]),
             subs(var=qm,recps[2])-subs(var=qn,recqs[2]),
             subs(var=qm,recps[3])]; # right hand side of (2.8)

      return(map(factor,[seq(c31*u[k]+c32*v[k]-b3[k],k=1..3)]));
    fi;
  fi;

end proc:

recgenNq:=proc(recp,recq,recpp,recqp,var,qm,qn,q)

  local c11,c12,b1,c21,c22,b2,res,u,v,
        delta,xi,eta,deltap,xip,etap,detc,
  recps,recqs,c31,c32,b3;

  c11:=subs(var=qm,recp[1]):
  c12:=subs(var=qm,recp[3]):
  b1:=[subs(var=qn/q,recq[1]),
        subs(var=qn,recq[2])-subs(var=qm,recp[2]),
        subs(var=qn*q,recq[3])];

  c21:=subs(var=qm,recpp[1]):
  c22:=subs(var=qm,recpp[3]):
  b2:=[subs(var=qn/q,recqp[1]),
        subs(var=qn,recqp[2])-subs(var=qm,recpp[2]),
        subs(var=qn*q,recqp[3])];

  detc:=factor(c11*c22-c12*c21);
  if detc=0 then
    return([seq(c11*b2[k]-c21*b1[k],k=1..3)]);
  else
    u:=[seq((c22*b1[k]-c12*b2[k])/detc,k=1..3)];
    v:=[seq((-c21*b1[k]+c11*b2[k])/detc,k=1..3)];

    if nargs=8 then

      delta:=factor(u[1]): xi:=factor(u[2]):
      eta:=factor(u[3]):

      deltap:=factor(v[1]): xip:=factor(v[2]):
      etap:=factor(v[3]):

      return(map(factor,[subs(qm=qm/q,eta)*subs(qn=qn*q,etap),
          subs(qm=qm/q,xi)*etap+subs(qm=qm/q,eta)*subs(qn=qn*q,xip),
          subs(qm=qm/q,delta)*subs(qn=qn/q,etap)+subs(qm=qm/q,xi)*xip
             +subs(qm=qm/q,eta)*subs(qn=qn*q,deltap)-1,
          subs(qm=qm/q,delta)*subs(qn=qn/q,xip)+subs(qm=qm/q,xi)*deltap,
          subs(qm=qm/q,delta)*subs(qn=qn/q,deltap)]));

    elif nargs=10 then

      recps:=args[9]: recqs:=args[10]:

      c31:=subs(var=qm,recps[1]):
      c32:=subs(var=qm,recps[3]):
      b3:=[subs(var=qn/q,recqs[1]),
             subs(var=qn,recqs[2])-subs(var=qm,recps[2]),
             subs(var=qn*q,recqs[3])];

      return(map(factor,[seq(c31*u[4-k]+c32*v[4-k]-b3[4-k],k=1..3)]));
    fi;
  fi;

end proc:


text:
kernelopts(printbytes=false);
with(RigidMotionsParameterSpaceDecompostion);
LaunchComputeSamplePoints([a,b,c], "__REPLACE__", "N2");
done

max_stars_repo_name: copyme/MapleTools

max_issues_repo_name: copyme/MapleTools

max_forks_repo_name: copyme/MapleTools


text:
program gcd;
begin
var n1, n2: integer;
read(n1);
read(n2);
if n2 > n1 then
begin
var tmp: integer;
tmp := n1;
n1 := n2;
n2 := tmp;
end;
while n2 > 0 do
begin
var new: integer;
new := n1 % n2;
n1 := n2;
n2 := new;
end;
writeln(n1);
end.

max_stars_repo_name: sebschmi/codegen-course

max_issues_repo_name: sebschmi/codegen-course

max_forks_repo_name: sebschmi/codegen-course


text:
program assignboolean;
var x

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing


text:
# typedef struct ss {
#   int f1;
#   char f2:3;
#   char f3:5;
# } SS;
# SS foo(SS x) {
#   x.f2 = 32;
#   return x;
# }
type $SS <struct {@f1 i32, @f2:3 i8, @f3:5 i8}>
func $foo (
	var %x <$SS>) i32 {
  dassign %x 2 ( constval i32 32 )
  return ( dread i32 %x 2 ) }
 # 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:
# Recall that $\Lambda^*$ is defined to be the monoid of lists of
# finite sets of natural numbers.  Such a list can be represented
# in Maple like [{1,4,6},{2,3,4,5},{1}], for example.  Similarly,
# if $N$ is a finite set of natural numbers, then $\Lambda^*[N]$
# is the submonoid generated by the subses of $N$.

# The function `is_element/itloc`(N)(AA) will return true iff
# AA is an element of $\Lambda^*[N]$.

`is_element/itloc` := (N::set) -> proc(AA)
 return type(AA,list(set(nonnegint)));
end:

# The function `is_small_element/itloc`(N)(l)(AA) will return true
# iff AA is an element of $\Lambda^*[N]$ and has length $l$.

`is_small_element/itloc` := (N::set) -> (l::nonnegint) -> proc(AA)
 local A;

 if not(`is_element/itloc`(N)(AA)) then
  return false;
 fi;

 if nops(AA) > l then return false; fi;
 A := map(op,{op(AA)});
 if A minus N <> {} then return false; fi;
 return true;
end:

# The function `list_small_elements/itloc`(N)(l) will produce a
# list of all the elements of $\Lambda^*[N]$ that have length l.
# (This will obviously be unmanageably large unless N and l are
# very small.)

`list_small_elements/itloc` := (N::set) -> proc(l::nonnegint)
 option remember;
 local PN,L,i;

 PN := combinat[powerset](N);
 if l = 0 then
  return [[]];
 else
  L := `list_small_elements/itloc`(N)(l-1);
  return [seq(seq([op(AA),A],A in PN),AA in L)];
 fi;
end:

# The function `list_smaller_elements/itloc`(N)(l) lists
# elements of length <= l.

`list_smaller_elements/itloc` := (N::set) -> proc(l::nonnegint)
 [seq(op(`list_small_elements/itloc`(N)(k)),k=0..l)];
end:

# The function `count_small_elements/itloc`(N)(l) returns
# the number of elements of length l.

`count_small_elements/itloc` := (N::set) -> (l::nonnegint) ->
 2^(nops(N) * l);

# The function `random_small_element/itloc`(N)(l)() returns
# a randomly chosen element of $\Lambda^*[N]$ of length l.
# (Note the extra pair of brackets in the syntax.)

`random_small_element/itloc` := (N::set) -> (l::nonnegint) -> proc()
 local r;

 r := `random_element/subsets`(N);
 return [seq(r(),i=1..l)];
end:

# The function `o/itloc` is the multiplication law for the
# monoid $\Lambda^*$: so for example `o/itloc`(AA,BB,CC) is
# just the list obtained by concatenating AA, BB and CC.

`o/itloc` := () -> map(op,[args]);

# `id/itloc` is the empty list, which is the identity element
# of the monoid $\Lambda^*$.

`id/itloc` := [];

# `omega/itloc`(AA) is the union of all the sets in the list AA.
# `cap/itloc`(AA,B) is the sequence of intersections of sets in
# AA with the fixed set B. (This is defn-om in the LaTeX file.)

`omega/itloc` := (AA) -> map(op,{op(AA)});
`cap/itloc` := (AA,B) -> map(A -> `intersect`(A,B),AA);

# `sigma/itloc`(N)(AA) is the intersection of all the sets in the
# list AA, to be interpreted as N if the list AA is empty.
# This is essentially the same as in rem-sg-xi in the LaTeX file.

`sigma/itloc` := (N::set) -> (AA) -> `intersect`(N,op(AA));

# `prec/itloc`(A,B) returns true if every element of A is strictly
# less than every element of B.  The function `preceq/itloc`(A,B)
# returns true if every element of A is less than or equal to
# every element of B.  These are in defn-prec in the LaTeX file.
`prec/itloc` := proc(A::set,B::set)
 if A = {} or B = {} then return true; fi;
 return evalb(max(op(A)) < min(op(B)));
end:

`preceq/itloc` := proc(A::set,B::set)
 if A = {} or B = {} then return true; fi;
 return evalb(max(op(A)) <= min(op(B)));
end:

# The next five definitions correspond to defn-reduced in the LaTeX file.

`is_top_reduced/itloc` := proc(AA)
 local t;
 t := map(A -> max(-infinity,op(A)),AA);
 return evalb(t = sort(t));
end:

`is_bottom_reduced/itloc` := proc(AA)
 local b;
 b := map(A -> min(op(A),infinity),AA);
 return evalb(b = sort(b));
end:

`is_unnested/itloc` := proc(AA)
 local r,i;
 r := nops(AA);
 for i from 1 to r-1 do
  if AA[i] minus AA[i+1] = {} then
   return false;
  fi;
  if AA[i+1] minus AA[i] = {} then
   return false;
  fi;
 od;
 return true;
end:

`is_end_reduced/itloc` := proc(AA)
 `is_top_reduced/itloc`(AA) and
 `is_bottom_reduced/itloc`(AA);
end:

`is_reduced/itloc` := proc(AA)
 `is_top_reduced/itloc`(AA) and
 `is_bottom_reduced/itloc`(AA) and
 `is_unnested/itloc`(AA);
end:

# The next three definitions correspond to defn-reductions and
# lem-trimming-confluent
`beta/itloc` := proc(AA)
 local r,i,n,BB;

 r := nops(AA);
 if r = 0 then return []; fi;
 BB := table();
 BB[1] := AA[1];
 for i from 2 to r do
  if BB[i-1] = [] then
   return [[]];
  fi;
  n := min(BB[i-1]);
  BB[i] := select(j -> evalb(j >= n),AA[i]);
 od;
 return [seq(BB[i],i=1..r)];
end:

`tau/itloc` := proc(AA)
 local r,i,n,BB;

 r := nops(AA);
 if r = 0 then return []; fi;
 BB := table();
 BB[r] := AA[r];
 for i from r-1 to 1 by -1 do
  if BB[i+1] = [] then
   return [[]];
  fi;
  n := max(BB[i+1]);
  BB[i] := select(j -> evalb(j <= n),AA[i]);
 od;
 return [seq(BB[i],i=1..r)];
end:

`delta/itloc` := proc(AA)
 local r,m,i,BB,CC,Bi,Cm,U,V,ok;

 r := nops(AA);
 if r = 0 then return []; fi;

 ok := false;
 CC := AA;
 while not(ok) do
  BB := CC;
  r := nops(BB);
  CC := [BB[1]];
  m := 1;
  for i from 2 to r do
   Bi := {op(BB[i])};
   Cm := {op(CC[m])};
   U := Bi minus Cm;
   V := Cm minus Bi;
   if U = {} then
    CC := [seq(CC[j],j=1..m-1),BB[i]];
   elif V <> {} then
    CC := [op(CC),BB[i]];
    m := m+1;
   fi;
  od:
  ok := evalb(CC = BB);
 od;

 return BB;
end:

`end_reduce/itloc` := (AA) -> `tau/itloc`(`beta/itloc`(AA)):

`reduce/itloc` := (AA) -> `delta/itloc`(`tau/itloc`(`beta/itloc`(AA))):

`rho/itloc` := eval(`reduce/itloc`):

# `oo/itloc` concatenates its arguments, and then applies the
# reduction map.

`oo/itloc` := proc()
 `reduce/itloc`(`o/itloc`(args));
end:

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

# The set $\Theta(AA)$ of threads is defined in defn-threads.
# The function `is_element/threads/itloc`(AA)(a) returns true
# iff a is a thread for AA.

`is_element/threads/itloc` := (AA) -> proc(a)
 local r,i;

 r := nops(AA);
 if not(type(a,list) and nops(a) = r) then
  return false;
 fi;

 for i from 1 to r do
  if not(member(a[i],AA[i])) then
   return false;
  fi;
 od;

 if a <> sort(a) then
  return false;
 fi;

 return true;
end:

# The function `is_leq/threads/itloc`(AA)(a,b) returns true iff
# a[i] <= b[i] for all i.

`is_leq/threads/itloc` := (AA) -> proc(a,b)
 local i;
 for i from 1 to nops(AA) do
  if a[i] > b[i] then return false; fi;
 od;
 return true;
end:

# The function `list_elements/threads/itloc`(AA) returns a list
# of all the threads for AA.

`list_elements/threads/itloc` := proc(AA)
 local r,AA0,L0,L,U,a0;

 r := nops(AA);

 if r = 0 then
  return [[]];
 elif r = 1 then
  return [op(map(a -> [a],AA[1]))];
 fi;

 AA0 := [seq(AA[i],i=1..r-1)];
 L0 := `list_elements/threads/itloc`(AA0);
 L := NULL;
 for a0 in L0 do
  U := select(a -> a0[r-1] <= a,AA[r]);
  L := L,seq([op(a0),a],a in U);
 od:
 return [L];
end:

# The function `random_element/threads`(AA)() attempts to return
# a randomly chosen thread for AA.  It is not very intelligent,
# so it may fail to find a thread even if one exists.  It will
# return the symbol FAIL if it does not find a thread.

`random_element/threads/itloc` := (AA) -> proc(num_tries::posint := 10)
 local k,n,a,i,B;

 if min(op(map(nops,AA))) = 0 then return FAIL; fi;

 k := num_tries;
 n := nops(AA);

 while k > 0 do
  k := k-1;
  a := [random_element_of(AA[1])];
  for i from 2 to n do
   B := select(b -> b >= a[-1],AA[i]);
   if B = {} then
    break;
   else
    a := [op(a),random_element_of(B)];
   fi;
  od;
  if nops(a) = n then return a; fi;
 od;

 return FAIL;
end:

# We do not know an efficient method to compute the number of threads,
# so we set `count_elements/threads/itloc` to NULL.
`count_elements/threads/itloc` := NULL:

# If the poset of threads is nonempty then it will have a smallest
# element, which will be returned by the function
# `bot/threads/itloc`(AA).

`bot/threads/itloc` := proc(AA)
 local BB;
 BB := `beta/itloc`(AA);
 if member({},BB) then
  return FAIL;
 else
  return map(min,BB);
 fi;
end:

# If the poset of threads is nonempty then it will have a largest
# element, which will be returned by the function
# `top/threads/itloc`(AA).

`top/threads/itloc` := proc(AA)
 local CC;
 CC := `tau/itloc`(AA);
 if member({},CC) then
  return FAIL;
 else
  return map(max,BB);
 fi;
end:

# `is_threadable/itloc`(AA) returns true iff AA has a thread
`is_threadable/itloc` :=
 (AA) -> (`bot/threads/itloc`(AA) <> FAIL);

# `min_thread_sets/itloc`(AA) returns the minimal elements of the
# poset of thread sets for AA (as in defn-threads)
`min_thread_sets/itloc` := proc(AA)
 local U,V,W,a,b;
 U := {op(`list_elements/threads/itloc`(AA))};
 U := map(u -> {op(u)},U);
 U := [op(U)];
 U := sort(U,(a,b) -> evalb(nops(a) < nops(b)));
 V := [];
 for a in U do
  W := select(b -> b minus a = {},V);
  if W = [] then V := [op(V),a]; fi;
 od;
 return {op(V)};
end:

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

# $\Phi_0[N]$ is the submonoid of $\Lambda[N]$ generated by the
# functors $L_{K(n)}$.  It bijects naturally with the subset of
# $\Lambda^*[N]$ consisting of the zero element [{}], together
# with all elements [{a[1]},...,{a[r]}] with a[1] <= ... <= a[r].

`is_element/Phi0/itloc` := (N::set(nonnegint)) -> proc(AA)
 local P;

 if not(`is_element/itloc`(N)(AA)) then
  return false;
 fi;

 if AA = [{}] then return true; fi;

 if {1,op(map(nops,AA))} <> {1} then
  return false; # not all singletons
 fi;

 P := map(op,AA);
 if sort(P) <> P then
  return false;
 fi;

 return true;
end:

`is_leq/Phi0/itloc` := (N::set(nonnegint)) -> proc(AA,BB)
 local P,Q;

 if AA = [{}] then return true; fi;
 if BB = [{}] then return false; fi;
 P := {op(map(op,AA))};
 Q := {op(map(op,AA))};

 return evalb(P minus Q = {});
end:

`list_elements/Phi0/itloc` := proc(N::set(nonnegint))
 local L;
 L := combinat[choose](sort([op(N)]));
 L := map(P -> map(i -> {i},P),L);
 L := [[{}],op(L)];
 return L;
end:

`random_element/Phi0/itloc` := proc(N::set(nonnegint))
 local P,AA;

 if rand(1..6)() = 1 then
  return [{}];
 else
  P := sort([op(random_subset_of(N)())]);
  AA := map(u -> {u},P);
  return AA;
 fi;
end:

`count_elements/Phi0/itloc` := proc(N::set(nonnegint))
 1 + 2^nops(N);
end:

`o/Phi0/itloc` := proc()
 local AA,P;

 if member([{}],[args]) then
  return [{}];
 fi;

 AA := map(op,[args]);
 P := map(op,AA);
 if P <> sort(P) then
  return [{}];
 fi;

 return AA;
end:

`id/Phi0/itloc` := [];

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

# M11/itloc(N)(A) is the set $M'_{N,A}$ from defn-M-prime

`is_element/M11/itloc` := (N::set(nonnegint)) -> (A::set(nonnegint)) -> proc(T)
 evalb(type(T,set) and (T minus N = {}) and (T intersect A <> {}));
end:

`is_leq/M11/itloc` := (N) -> (A) -> (T0,T1) -> evalb(T0 minus T1 = {}):

`list_elements/M11/itloc` := (N::set(nonnegint)) -> proc(A::set(nonnegint))
 return [op(select(T -> T intersect A <> {},combinat[powerset](N)))];
end:

`random_element/M11/itloc` := (N::set(nonnegint)) -> (A::set(nonnegint)) -> proc()
 local U,V;

 U := `random_element/nonempty_subsets`(A)();
 V := `random_element/subsets`(N minus A)();
 return U union V;
end:

`count_elements/M11/itloc` := (N::set(nonnegint)) -> proc(A::set(nonnegint))
 2^nops(N) - 2^nops(N minus A);
end:

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

# M1/itloc(N)(AA) is the set $M'_{N,AA}$ from defn-M-prime-multi

`is_element/M1/itloc` := (N::set(nonnegint)) -> (AA) -> proc(TT)
 local r,i;

 r := nops(AA);

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

 if nops(TT) <> r then
  return false;
 fi;

 if map(op,{op(TT)}) minus N <> {} then
  return false;
 fi;

 for i from 1 to r do
  if TT[i] intersect AA[i] = {} then
   return false;
  fi;
 od;

 for i from 1 to r-1 do
  if not(`preceq/itloc`(TT[i],TT[i+1])) then
   return false;
  fi;
 od;

 return true;
end:

`is_leq/M1/itloc` := (N) -> (AA) -> proc(TT0,TT1)
 `and`(seq(TT0[i] minus TT1[i] = {},i=1..nops(TT0)));
end:

`list_elements/M1/itloc` := (N::set(nonnegint)) -> proc(AA)
 local r,m,L,L0,L1,L2,AA0,TT0;
 r := nops(AA);
 if r = 0 then
  return [[]];
 fi;

 L1 := `list_elements/M11/itloc`(N)(AA[r]);

 if r = 1 then
  return map(T -> [T],L1);
 fi;

 AA0 := [seq(AA[i],i=1..r-1)];
 L0 := `list_elements/M1/itloc`(N)(AA0);
 L := NULL;
 for TT0 in L0 do
  m := max(TT0[r-1]);
  L2 := select(T -> min(T) >= m,L1);
  L := L,seq([op(TT0),T],T in L2);
 od:
 L := [L];
 return L;
end:

`random_element/M1/itloc` := (N::set(nonnegint)) -> (AA) -> proc(num_tries := 10)
 local r,k,AA0,TT0,m,N0,A0,T;
 r := nops(AA);
 if r = 0 then
  return [];
 fi;

 k := num_tries;
 while k > 0 do
  k := k-1;
  AA0 := [seq(AA[i],i=1..r-1)];
  TT0 := `random_element/M1/itloc`(N)(AA0)(num_tries);
  if TT0 <> FAIL then
   N0 := N;
   if r > 1 then
    m := max(op(TT0[r-1]));
    N0 := select(i -> (i >= m),N);
   fi;
   A0 := N0 intersect AA[r];
   if A0 <> {} then
    T := `random_element/M11/itloc`(N0)(A0)();
    if T <> FAIL then
     return [op(TT0),T];
    fi;
   fi;
  fi;
 od;

 return FAIL;
end:

`count_elements/M1/itloc` := NULL:

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

`is_element/M2/itloc` := (N::set(nonnegint)) -> (AA) -> proc(T::set(nonnegint))
 `is_threadable/itloc`(`cap/itloc`(AA,T)) and
 (T minus N) = {};
end:

`is_leq/M2/itloc` := (N::set(nonnegint)) -> (AA) -> proc(T0,T1)
 evalb(T0 minus T1 = {});
end:

`list_elements/M2/itloc` := (N::set(nonnegint)) -> proc(AA)
 local L;

 L := `list_elements/threads/itloc`(AA);
 L := map(a -> {op(a)},L);
 return [seq(seq(T union U,U in combinat[powerset](N minus T)),T in L)];
end:

`random_element/M2/itloc` := (N::set(nonnegint)) -> (AA) -> proc()
 local a,T;
 a := `random_element/threads/itloc`(AA)();
 if a = FAIL then return FAIL; fi;
 T := {op(a)};
 T := {op(T),op(random_subset_of(N minus T))};
 return T;
end:

`count_elements/M2/itloc` := NULL:

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

`is_element/M3/itloc` := (N::set(nonnegint)) -> (AA) -> (T) -> proc(TT)
 `is_element/M1/itloc`(N)(AA)(TT) and
  `omega/itloc`(TT) minus T = {};
end:

`is_leq/M3/itloc` := (N::set(nonnegint)) -> (AA) -> (T) -> proc(TT0,TT1)
 `and`(seq(evalb(TT0[i] minus TT1[i] = {}),i=1..nops(TT0)));
end:

`list_elements/M3/itloc` := (N::set(nonnegint)) -> (AA) -> proc(T)
 local BB,L;

 BB := `cap/itloc`(AA,T);
 L := `list_elements/M1/itloc`(T)(BB);
 return L;
end:

`phi/itloc` := (N::set(nonnegint)) -> (AA) -> (T) -> (TT) -> proc(k)
 local r,i,BB,b;
 r := nops(AA);
 if k < 0 or k > 2*r then error("k is out of range"); fi;
 i := ceil(k/2);
 BB := `beta/itloc`(`cap/itloc`(AA,T));
 if member({},BB) then error("T is not a thread set for AA"); fi;
 b := map(min,BB);
 if k = 0 then
  return TT;
 elif k = 2*i then
  return [seq({b[j]},j=1..i),seq(TT[j],j=i+1..r)];
 else
  return [seq({b[j]},j=1..i-1),TT[i] union {b[i]},seq(TT[j],j=i+1..r)];
 fi;
end:


text:
program typeconv;
 var i : integer; b : boolean; c : char;
begin
  i := integer(false);
  while i <= integer(true) do begin
    writeln( boolean(i), ' : ', i);
    i := i + 1;
  end;
  writeln;
  i := integer(' ');
  while i < 127 do begin
    if i div 16 * 16 = i then writeln;
    write(char(i), ' ');
    i := i + 1;
  end;
  writeln
end.


text:
with(LinearAlgebra);
cdir := "/home/scott/wDocuments/research/software/hqca/hqca/maple/";
data := readdata(cat(cdir, "temp.rdm"), 6);
with(ArrayTools);
#
Rearrange := proc(x) local n, a, i, j, k, l; `local`(a, n, i, j, k, l); `description`("convert chemists to numpy and flatten an array to form a matrix"); n := Size(x); a := Array(1 .. n[1], 1 .. n[2], 1 .. n[3], 1 .. n[4], datatype = float[8]); for i to n[1] do for j to n[2] do for k to n[3] do for l to n[4] do a[i, k, j, l] := x[i, j, k, l]; end do; end do; end do; end do; return a; end proc;
Flatten := proc(x) local n, a, i, j, k, l; `local`(a, n, i, j, k, l); `description`("convert chemists to numpy and flatten an array to form a matrix"); n := Size(x); a := Array(1 .. n[1]*n[2], 1 .. n[3]*n[4], datatype = float[8]); for i to n[1] do for j to n[2] do for k to n[3] do for l to n[4] do a[(i - 1)*n[1] + j, (k - 1)*n[3] + l] := x[i, k, j, l]; end do; end do; end do; end do; return a; end proc;
New := Array(1 .. 3, 1 .. 3, 1 .. 3, 1 .. 3, datatype = float[8]);
newdata := data[3 .. ()];
for i in newdata do
    New[round(i[1]), round(i[2]), round(i[4]), round(i[3])] := i[5];
end do;
New[1, 2];
with(QuantumChemistry);

mol := [["H", 1.00000000, 0, 0], ["H", -1.00000000, 0, 0]];

h2 := Variational2RDM(mol, basis = "STO-3G", return_rdm = "rdm1_and_rdm2");

h2[rdm2];

pure := Purify2RDM(h2[rdm2], spin_free = true, electron_number = 2, conv_tol = 0.10000000, conditions = "DQG");
^2 D ^{i,j}_{k,l} = < i j l k  > -> A[1..r,1..r,1..r,1..r] <=> A[i,j,k,l]
;

G := Flatten(h2[rdm2]);

Eigenvalues(G);
F := Flatten(pure[rdm2]);

Eigenvalues(F);

max_stars_repo_name: damazz/HQCA

max_issues_repo_name: damazz/HQCA

max_forks_repo_name: damazz/HQCA


text:
var $g <*ptr>
var $h <*ptr> = addrof ptr $g
var $g <*ptr> = addrof ptr $h
var $h <*ptr>
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


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

`is_element/linord` := (A::set) -> proc(L)
 global reason;

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

 if nops(L) <> nops(A) then
  reason := ["is_element/linord","L does not the same length as A",L,A];
  return false;
 fi;

 if {op(L)} <> A then
  reason := ["is_element/linord","L is not an enumeration of A",L,A];
  return false;
 fi;

 return true;
end;

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

`is_equal/linord` := (A::set) -> proc(L0,L1)
 return evalb(L0 = L1);
end;

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

`op/linord` := (A::set) -> proc(L)
 local n,i;
 n := nops(A);
 return [seq(L[n-i],i=0..n-1)];
end;

`res/linord` := (A::set,B::set) -> proc(L)
 return select(a -> member(a,B),L);
end:

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

`flip/linord` := (A::set) -> (L) -> proc(a)
 local n,i;
 n := nops(L);
 if n <= 1 then return L; fi;
 i := 1;
 while i <= n and L[i] <> a do i := i+1; od;
 if i > n then error("a is not in L"); fi;
 return [op(i..n,L),op(1..i-1,L)];
end:

`loop/linord` := (A::set) -> (L) -> L;

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

`random_element/linord` := (A::set) -> proc()
 return combinat[randperm](A);
end:

`list_elements/linord` := proc(A::set)
 return combinat[permute](A);
end:

`count_elements/linord` := proc(A::set)
 return nops(A) !;
end:

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

`linord_is_leq` := (A::set) -> (L) -> proc(a,b)
 local x;
 for x in L do
  if x = a then return true; fi;
  if x = b then return false; fi;
 od;
 error("Neither a nor b is in L");
end:

`linord_is_less` := (A::set) -> (L) -> proc(a,b)
 return evalb((a <> b) and `linord_is_leq`(A)(L)(a,b));
end:


text:
TITLE

	CAP_extensive_form_cap_test_5;

STOCHASTIC

DATA

	I := DATAFILE("cap_test_5_SETS.dat");
	J := DATAFILE("cap_test_5_SETS.dat");
	K := DATAFILE("cap_test_5_SETS.dat");				! - num scenarios
	Q := I;							! - num 1st stage columns
	T := I*J;						! - num 2nd stage columns
	R := J+I;						! - num rows in T and W

SCENARIO

	k := 1..K;

INDEX
	FirstStCol := 1..Q;
	SecondStCol := 1..T;
	Row := 1..R;

PROBABILITIES

	p[k] :=  ALLEQUAL;

RANDOM DATA

	Vector_h[k, Row] := DATAFILE("cap_test_5_RANDOM.dat");	! - rows are related to scenarios

DATA

	Vector_c[FirstStCol] := DATAFILE("cap_test_5_DATA.dat");
	Vector_d[SecondStCol] := DATAFILE("cap_test_5_DATA.dat");
	Vector_b := DATAFILE("cap_test_5_DATA.dat");			! - max sum of h over all scenarios
	Matrix_A[FirstStCol] := DATAFILE("cap_test_5_DATA.dat");
	Matrix_T[Row, FirstStCol] := DATAFILE("cap_test_5_DATA.dat");
	Matrix_W[Row, SecondStCol] := DATAFILE("cap_test_5_DATA.dat");
	x_ub[FirstStCol] := DATAFILE("cap_test_5_DATA.dat");

STAGE1 BINARY VARIABLES

	x[FirstStCol];

STAGE2 VARIABLES

	y[SecondStCol];

MACROS

	StageFirst = sum(FirstStCol: Vector_c * x);
	StageSecond = sum(SecondStCol: Vector_d * y);

MIN

	Objective = StageFirst + StageSecond;

SUBJECT TO

	FirstStCnstr: sum(FirstStCol: Matrix_A * x) >= Vector_b;
    	SecondStCnstr[Row]: sum(FirstStCol: Matrix_T * x) + sum(SecondStCol: Matrix_W * y) >= Vector_h;

BOUNDS

	BinBounds[FirstStCol]: x <= x_ub;	! - why do we need it

END

max_stars_repo_name: pashtetGP/optconvert

max_issues_repo_name: pashtetGP/optconvert

max_forks_repo_name: pashtetGP/optconvert


text:
# Differential geometry of curves in R^3

analyse_curve := proc(x,t)
 local T;
 T := table():
 T["x"] := x;
 T["t"] := t;
 T["v"] := map(diff,x,t);
 T["speed"] := simplify(`norm_2/R`(3)(T["v"]));
 T["s"] := int(T["speed"],t);
 T["v_hat"] := simplify(T["v"] /~ T["speed"]);
 T["T"] := T["v_hat"];
 T["a"] := map(diff,T["v"],t);
 T["v_hat_dot"] := simplify(map(diff,T["v_hat"],t));
 T["curvature"] := simplify(`norm_2/R`(3)(T["v_hat_dot"]));
 T["N"] := simplify(T["v_hat_dot"] /~ T["curvature"]);
 T["radius_of_curvature"] := simplify(1/T["curvature"]);
 T["B"] := `cross_product`(T["T"],T["N"]);

 return eval(T);
end:


text:
######################################################################
# epi(A,B) is the set of surjective maps from A to B

`is_element/epi` := (A::set,B::set) -> proc(p)
 local a,b,hit;
 global reason;

 if not `is_element/maps`(A,B)(p) then
  reason := [convert(procname,string),"p is not a map from A to B",p,A,B,reason];
  return false;
 fi;

 hit := table();
 for b in B do hit[b] := false; od;
 for a in A do hit[p[a]] := true; od;
 for b in B do if not hit[b] then
  reason := [convert(procname,string),"b is not hit",b];
  return false;
 fi; od;
 return true;
end;

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

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

 return true;
end;

`is_leq/epi` := NULL;

`random_element/epi` := (A::set,B::set) -> proc()
 local p,m,n,A0,r,i;

 n := nops(A);
 m := nops(B);
 if m > n then return FAIL; fi;

 p := table();
 A0 := combinat[randperm]([op(A)]);

 for i from 1 to m do p[A0[i]] := B[i]; od;

 r := rand(1..m);

 for i from m+1 to n do p[A0[i]] := B[r()]; od;

 return eval(p);
end;

`list_elements/epi` := proc(A::set,B::set)
 local n,m,B0,b,E,E0,k,CC,C,c,p0,p;

 n := nops(A);
 m := nops(B);

 if n < m then return []; fi;
 if m = 0 then
  if n = 0 then
   return [table()];
  else
   return [];
  fi;
 fi;

 B0 := sort(B);
 b := B0[1];
 B0 := {op(2..-1,B0)};

 E := [];

 for k from 1 to n-m+1 do
  CC := combinat[choose](A,k);
  for C in CC do
   E0 := `list_elements/epi`(A minus C,B0);
   for p0 in E0 do
    p := copy(p0);
    for c in C do p[c] := b; od;
    E := [op(E),eval(p)];
   od;
  od;
 od;

 return E;
end;

`count_elements/epi` := (A::set,B::set) -> Stirling2(nops(A),nops(B)) * nops(B)!;


text:
var nTimes : int := 0;
print "How many times? ";
read nTimes;
var x : int;
for x in 0..nTimes-1 do
    print x;
    print " : Hello, World!\n";
end for;
assert (x = nTimes);

max_stars_repo_name: jgke/miniplc

max_issues_repo_name: jgke/miniplc

max_forks_repo_name: jgke/miniplc


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

(* Equation (7) from the paper *)

tpss_p := x -> X2S^2*x^2:
tpss_z := (x, t) -> x^2/(8*t):

tpss_alpha := (x, t) -> (t - x^2/8)/K_FACTOR_C:
tpss_qb    := (x, t) -> \
      9/20 * (tpss_alpha(x, t) - 1)/sqrt(1 + params_a_b*tpss_alpha(x, t)*(tpss_alpha(x, t) - 1)) \
      + 2*tpss_p(x)/3:

(* Equation (10) in all its glory *)
tpss_fxnum := (x, t) -> \
      + (MU_GE + params_a_c*tpss_z(x, t)^tpss_ff(tpss_z(x, t))/(1 + tpss_z(x, t)^2)^2)*tpss_p(x) \
      + 146/2025 * tpss_qb(x, t)^2 \
      - 73/405 * tpss_qb(x, t) * sqrt(1/2*(9/25*tpss_z(x, t)^2 + tpss_p(x)^2)) \
      + MU_GE^2/tpss_kappa(x, t) * tpss_p(x)^2 \
      + 2*sqrt(params_a_e)*MU_GE*9/25*tpss_z(x, t)^2 \
      + params_a_e*params_a_mu*tpss_p(x)^3:

tpss_fxden := x -> \
      (1 + sqrt(params_a_e)*tpss_p(x))^2:

tpss_fx    := (x, t) -> tpss_fxnum(x, t)/tpss_fxden(x):


text:
# float a[10];
# void init(void) {
#   int i;
#   for (i=0; i<10; i++)
#     a[i] = i*3;
# }

var $a <[10] f32>
var $f1 f32
func $init () void {
  var %i i32
  dassign %i ( constval i32 0 )
  while (
    lt i32 i32 ( dread i32 %i, constval i32 10) ) {
      iassign <* [10] f32> (
        array 1 a32 <* [10] f32> (addrof a32 $a, dread i32 %i),
	      mul f32 (dread f32 $f1, constval f32 3.33f) )
      dassign %i (
        add i32 (dread i32 %i, constval i32 1) )} }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
# This is a generic bracket operation.  With more than two arguments,
# it represents an iterated bracket associated like [a,[b,[c,[d,e]]]].
# Multilinearity rules are applied automatically, and the Jacobi
# rule is also used automatically to convert other kinds of iterated
# brackets into the fully right associated form.

bracket := proc()
 local n,m,i,j,u,a,b,c,v,ac,av,x,y;
 n := nargs;
 for i from 1 to n do
  u := seq(args[j],j=1..i-1);
  a := args[i];
  v := seq(args[j],j=i+1..n);
  if type(a,`+`) then
   return map(t -> bracket(u,t,v),a);
  elif type(a,`*`) then
   ac,av := selectremove(type,a,rational);
   if ac <> 1 then
    return ac * bracket(u,av,v);
   fi;
  elif type(a,specfunc(bracket)) then
   if i = n then
    return bracket(u,op(a));
   else
    b := op(1,a);
    c := bracket(op(2..-1,a));
    if nops(c) = 1 then c := op(c); fi;
    return bracket(u,b,c,v) - bracket(u,c,b,v);
   fi;
  fi;
 od;

 if type(args[n],specfunc(bracket)) then
  return bracket(args[1..n-1],op(args[n]));
 fi;

 return 'procname'(args);
end:


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

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

(* type: work_mgga_x *)

a := 0.5389:
b := 3:

p := x -> X2S^2*x^2:
q := u -> X2S^2*u:

(* Equation (15) *)
fab := z -> piecewise( \
    z<=0, 0, \
    z>=a, 1, \
    (1 + exp(a/(a-z)))^b/(exp(a/z) + exp(a/(a-z)))^b \
):

(* Equation (7) *)
mDelta := (x, u) -> 8*q(u)^2/81 - p(x)*q(u)/9 + 8*p(x)^2/243:

f_W    := x -> 5*p(x)/3:

(* Equation (8) *)
f_GE4  := (x, u) -> 1 + 5*p(x)/27 + 20*q(u)/9 + mDelta(x, u):

(* Equation (11) *)
f_GE4_M := (x, u) -> f_GE4(x, u)/sqrt(1 + mDelta(x, u)^2/(1 + f_W(x))^2):

(* Equation (17) *)
f   := (rs, x, t, u) -> f_W(x) + (f_GE4_M(x, u) - f_W(x))*fab(f_GE4_M(x, u) - f_W(x)):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
# Center of Mass Calculation for the Robot based on MDH frames
# Introduction
# Berechnung der Schwerpunkts-Kinematik (Positionen)
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# floatb -> floating base wird durch base twist (Geschwindigkeit der Basis) oder vollständige Orientierung (Euler-Winkel) berücksichtigt
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# kinematics_com -> Kinematik der Körperschwerpunkte
# worldframe -> Berechnung der Geschwindigkeit im Welt-KS (KSW)
# par1 -> Parametersatz 1 (Schwerpunkt als Parameter: SX,SY,SZ)
# Testergebnisse
# Test des exportierten Matlab-Codes erfolgreich (kinetische Energie)
# Authors
# Moritz Schappler, schappler@irt.uni-hannover.de, 2014-11
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Sources
# [GautierKhalil1990] Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots
# [KhalilDombre2002] Modeling, Identification and Control of Robots
# [Ortmaier2014] Vorlesungsskript Robotik I (WS 2014/15)
# [Ott2008] Cartesian Impedance Control of Redundant and Flexible-Joint Robots
# Initialization
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
codegen_act := true:
codegen_opt := 2:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotx":
read "../transformation/proc_troty":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../transformation/proc_trafo_mdh":
read "../robot_codegen_definitions/robot_env":
printf("%s. Generiere Schwerpunktskinematik für %s\n", FormatTime("%Y-%m-%d %H:%M:%S"), robot_name):
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
read sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
kin_constraints_exist := kin_constraints_exist: # nur zum Abschätzen der Komplexität
;
if not assigned(simplify_options) or simplify_options(3)=-1 then # Standard-Einstellungen:
  if not kin_constraints_exist then # normale serielle Ketten und Baumstrukturen
    use_simplify := 0: # Standardmäßig aus
  else # mit kinematischen Zwangsbedingungen
    use_simplify := 1: # standardmäßig simplify-Befehle anwenden
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(3): # dritter Eintrag ist für CoM-Kinematik
end if:

# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name, base_method_name):
Trf := Trf:
Trf_c := Trf_c:

# Positions of Center of Mass
r_W_i_Si := Matrix(3, NL):
r_W_W_Si := Matrix(3, NL):
mr_W_i_Si := Matrix(3, NL):
# Es gibt NL Körper (Basis und NJ bewegliche Körper).
# r_i_i_Si(1 .. 3, i) ist der Vektor vom Ursprung des Körper-KS i zum Schwerpunkt des Körpers i ausgedrückt im Körper-KS i
# r_W_i_Si(1 .. 3, i) ist der Vektor vom Ursprung des Körper-KS i zum Schwerpunkt des Körpers i ausgedrückt im Basis-KS
# mr_W_i_Si(1 .. 3, i) ist der Vektor vom Ursprung des Körper-KS i zum Schwerpunkt des Körpers i ausgedrückt im Basis-KS multipliziert mit der Masse des Segments i
# r_W_W_Si(1 .. 3, i) ist der Vektor vom Ursprung des Basis-KS zum Schwerpunkt des Körpers i ausgedrückt im Basis-KS
# Trf_c(1 .. 3, 1 .. 3, i) ist das Körperkoordinatensystem zu Körper i im Basis KS.
for i to NL do
  r_W_i_Si(1 .. 3, i) :=  Multiply(Matrix(Trf_c(1 .. 3, 1 .. 3, i)), r_i_i_Si(1 .. 3, i)):
  if use_simplify>=1 then
    tmp_t11:=time():
    tmp_l11 := length(r_W_i_Si(1 .. 3, i)):
    printf("%s. Vereinfache Schwerpunktskoordinaten. Länge: %d (r_i_i_Si).\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l11):
    r_W_i_Si(1 .. 3, i) := simplify2(r_W_i_Si(1 .. 3, i)):
    tmp_l21 := length(r_W_i_Si(1 .. 3, i)):
    tmp_t21:=time():
  end if:
  mr_W_i_Si(1 .. 3, i) := Multiply(Matrix(Trf_c(1 .. 3, 1 .. 3, i)), mr_i_i_Si(1 .. 3, i)):
  r_W_W_Si(1 .. 3, i) :=  Matrix(Trf_c(1 .. 3, 4, i)) + Matrix(r_W_i_Si(1 .. 3, i)):
  printf("%s. Schwerpunktsposition in Weltkoordinaten für Körper %d aufgestellt.\n", FormatTime("%Y-%m-%d %H:%M:%S"), i-1):#0=Basis
  if use_simplify>=1 then
    tmp_t12:=time():
    tmp_l12 := length(mr_W_i_Si(1 .. 3, i)):
    tmp_l13 := length(r_W_W_Si(1 .. 3, i)):
    printf("%s. Vereinfache Schwerpunktskoordinaten. Länge: %d/%d (r_W_W_Si/mr_W_i_Si).\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l13, tmp_l12):
    mr_W_i_Si(1 .. 3, i) := simplify2(mr_W_i_Si(1 .. 3, i)):
    r_W_W_Si(1 .. 3, i) := simplify2(r_W_W_Si(1 .. 3, i)):
    tmp_l22 := length(mr_W_i_Si(1 .. 3, i)):
    tmp_l23 := length(r_W_W_Si(1 .. 3, i)):
    tmp_t22:=time():
    printf("%s. Terme für Schwerpunktskoord. vereinfacht. Länge: %d->%d / %d->%d / %d->%d (r_W_i_Si / mr_W_i_Si / r_W_W_Si). Rechenzeit %1.1fs und %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l11, tmp_l21, tmp_l12, tmp_l22, tmp_l13, tmp_l23, tmp_t21-tmp_t11, tmp_t22-tmp_t12):
  end if:
end do:
# Maple Export
save mr_W_i_Si, r_W_W_Si, r_W_i_Si, sprintf("../codeexport/%s/tmp/kinematics_com_worldframe_floatb_%s_par1_maple.m", robot_name, base_method_name):
printf("Maple-Ausdrücke exportiert.\n"):
# Matlab-Export
if codegen_act then
  MatlabExport(convert_t_s(r_W_W_Si), sprintf("../codeexport/%s/tmp/fkine_com_floatb_%s_matlab.m", robot_name, base_method_name), codegen_opt):
end if:
# Calculate CoM of the whole robot
c:=Vector(3):
M_ges:=0:
for i from 1 to NL do
  c:=c + r_W_W_Si(1 .. 3,i)*M[i,1]:
  M_ges:=M_ges + M[i,1]:
end do:
c:=c/M_ges:
# Maple-Export
save c, sprintf("../codeexport/%s/tmp/kinematics_com_total_worldframe_floatb_%s_par1_maple.m", robot_name, base_method_name):
# Matlab-Export
if codegen_act then
  MatlabExport(convert_t_s(c), sprintf("../codeexport/%s/tmp/com_total_worldframe_floatb_%s_par1_matlab.m", robot_name, base_method_name), codegen_opt):
end if:

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


text:
# This defines the chain complex version of the Barratt-Eccles operad.
# Basis elements of E(A) are represented by expressions T(...), where
# there is at least one argument, each argument is a list containing
# each element of A precisely once, and no two adjacent arguments are
# the same.

`is_element/barratt_eccles` := (A::set) -> proc(x)
 local y,z,i;

 if x = 0 then return true; fi;

 if type(x,`+`) then
  return `and`(op(map(`is_element/barratt_eccles`(A),[op(x)])));
 fi;

 if type(x,`*`) then
  y,z := selectremove(type,[op(x)],integer);
  return (nops(z) = 1 and `is_element/barratt_eccles`(A)(z[1]));
 fi;

 if not type(x,specfunc(T)) then return false; fi;

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

 for i from 1 to nops(x) do
  if not(`is_element/ord`(A)(op(i,x))) then
   return false;
  fi;
 od;

 for i from 1 to nops(x) - 1 do
  if op(i,x) = op(i+1,x) then return false; fi;
 od;

 return true;
end:

# Auxiliary function feeding into `diff/barratt_eccles`

`diff0/barratt_eccles` := (A::set) -> proc(x)
 local y,n,i;

 y := 0;
 n := nops(x);

 for i from 1 to n do
  if 1 < i and i < n and op(i-1,x) = op(i,x) then
   continue;
  fi;

  y := y + (-1)^(i-1) * T(op(1..i-1,x),op(i+1..n,x));
 od;

 return y;
end:

# Differential on the chain complex
`diff/barratt_eccles` := (A::set) -> apply_linear(`diff0/barratt_eccles`(A));

# Auxiliary function feeding into `deg/barratt_eccles`
`deg0/barratt_eccles` := (A) -> (x) -> nops(x) - 1;

# Degree function
`deg/barratt_eccles` := (A) -> apply_deg(`deg0/barratt_eccles`(A));

# This is the circle product for the operad structure.
# It is assumed that B is a subset of A and u is in E(A/B) and v is in E(B),
# where A/B is implemented as A \ B u {B}.

`o0/barratt_eccles` := (A,B) -> proc(u,v)
 local d,e,SS,x,s,p,w,u0,v0,w0,i;

 d := nops(u) - 1;
 e := nops(v) - 1;
 SS := `list_elements/shuffles`(d,e);

 x := 0;

 for s in SS do
  p := `to_grid_path/shuffles`(d,e)(s);
  w := NULL;
  for i from 0 to d+e do
   u0 := op(p[i][1]+1,u);
   v0 := op(p[i][2]+1,v);
   w0 := subs(B = op(v0),u0);
   w := w,w0;
  od:
  x := x + `sgn/shuffles`(d,e)(s) * T(w);
 od:

 return x;
end:

`o/barratt_eccles` := (A,B) -> apply_bilinear(`o0/barratt_eccles`(A,B));

# We now have various functions related to an interesting filtration
# of the operad.
`flip_count/barratt_eccles` := (A::set) -> proc(x)
 local m,a,b,k,rr,i,r0,r1;

 m := table():
 for a in A do
  for b in A do
   m[a,b] := 0;
  od;
 od;

 k := nops(x);
 rr := map(`rank_table/ord`(A),[op(x)]);

 for i from 1 to k-1 do
  r0 := rr[i];
  r1 := rr[i+1];

  for a in A do
   for b in A do
    if (r1[b] - r1[a]) * (r0[b] - r0[a]) < 0 then
     m[a,b] := m[a,b] + 1;
    fi;
   od;
  od;
 od;

 return eval(m);
end:

`flip_count_matrix/barratt_eccles` := (A::set) -> proc(u,s_)
 local m;
 m := `flip_count/barratt_eccles`(A)(args);
 return Matrix([seq([seq(m[a,b],b in A)],a in A)]);
end:

`max_flip_count/barratt_eccles` := (A::set) -> proc(x)
 local m,mm,a,b;

 m := `flip_count/barratt_eccles`(A)(x);

 mm := 0;
 for a in A do
  for b in A do
   mm := max(mm,m[a,b]);
  od;
 od;

 return mm;
end:

`is_member/barratt_eccles_cells` := (A::set) -> proc(ms)
 local m,s,AA,a,b;

 if not(type(ms,list) and nops(ms) = 2) then return false; fi;

 m,s := op(ms);

 if not(`is_element/ord`(A)(s)) then return false; fi;
 if not(type(m,table)) then return false; fi;

 AA := {seq(seq([a,b],b in A),a in A)};
 if {indices(m)} <> AA then return false; fi;

 for a in A do
  for b in A do
   if not(type(m[a,b],nonnegint)) then return false; fi;
  od:
 od:

 for a in A do
  if m[a,a] <> 0 then return false; fi;
  for b in A do
   if m[a,b] <> m[b,a] then return false; fi;
  od;
 od;

 return true;
end:

`is_cell_member/barratt_eccles` := (A::set) -> (ms) -> proc(x)
 local m0,m1,s,a,b,y,z;

 if x = 0 then return true; fi;
 if type(x,`+`) then
  return `and`(seq(`is_cell_member/barratt_eccles`(A)(ms)(y),y in x));
 fi;
 if type(x,`*`) then
  y,z := selectremove(type,[op(x)],integer);
  if nops(z) = 1 then
   return `is_cell_member/barratt_eccles`(A)(ms)(z[1]);
  else
   return false;
  fi;
 fi;

 m0,s := op(ms);
 m1 := `flip_count/barratt_eccles`(A)(T(op(x),s));

 for a in A do
  for b in A do
   if m1[a,b] > m0[a,b] then return false; fi;
  od;
 od;

 return true;
end:

`filtration0/barratt_eccles` := (A::set) -> proc(x)
 local m,m_max,s,rr,k,i,r0,r1,a,b;

 if type(x,integer) then return 0; fi;

 return 1 + `max_flip_count/barratt_eccles`(A)(x);
end:

`filtration/barratt_eccles` := (A::set) ->
 apply_max_deg(`filtration0/barratt_eccles`(A));


text:
# 1 + 8 -> 9
func $testSum (var %i i32) i32 {
  return (
    add i32(constval i32 1, constval i32 8))
}

# sqrt(4.0) -> 2.0
func $testSqrt (var %i i32) f32 {
  return (
    sqrt f32(constval f32 4.0f))
}

# 5 + (5 + a) -> a + 10
func &testConstVarFold00 (var %a i32) i32 {
  return (add i32 (
    constval i32 5,
    add i32 (constval i32 5, dread i32 %a)))
}

# 5 + (5 - a) -> 10 - a
func &testConstVarFold01 (var %a i32) i32 {
  return (add i32 (
    constval i32 5,
    sub i32 (constval i32 5, dread i32 %a)))
}

# 5 + (a + 5) ->  a + 10
func &testConstVarFold02 (var %a i32) i32 {
  return (add i32 (
    constval i32 5,
    add i32 (dread i32 %a, constval i32 5)))
}

# 6 + (a - 5) ->  a + 1
func &testConstVarFold03 (var %a i32) i32 {
  return (add i32 (
    constval i32 6,
    sub i32 (dread i32 %a, constval i32 5)))
}

# 5 * (5 * a) -> 5 * (5 * a)
func &testConstVarFold04 (var %a i32) i32 {
  return (mul i32 (
    constval i32 5,
    mul i32 (constval i32 5, dread i32 %a)))
}

# 5 * (a * 5) -> 5 * (a * 5)
func &testConstVarFold05 (var %a i32) i32 {
  return (mul i32 (
    constval i32 5,
    mul i32 (dread i32 %a, constval i32 5)))
}

# 3 - (2 - a) -> a + 1
func &testConstVarFold06 (var %a i32) i32 {
  return (sub i32 (
    constval i32 3,
    sub i32 (constval i32 2, dread i32 %a)))
}

# 3 - (2 + a) -> 1 - a
func &testConstVarFold07 (var %a i32) i32 {
  return (sub i32 (
    constval i32 3,
    add i32 (constval i32 2, dread i32 %a)))
}

# 3 - (a + 2) -> 1 - a
func &testConstVarFold08 (var %a i32) i32 {
  return (sub i32 (
    constval i32 3,
    add i32 (dread i32 %a, constval i32 2)))
}

# 3 - (a - 2) -> 5 - a
func &testConstVarFold09 (var %a i32) i32 {
  return (sub i32 (
    constval i32 3,
    sub i32 (dread i32 %a, constval i32 2)))
}

# (2 + a) - 3 -> a - 1
func &testConstVarFold10 (var %a i32) i32 {
  return (sub i32 (
    add i32 (constval i32 2, dread i32 %a),
    constval i32 3))
}

# (2 - a) - 3 -> -a - 1
func &testConstVarFold11 (var %a i32) i32 {
  return (sub i32 (
    sub i32 (constval i32 2, dread i32 %a),
    constval i32 3))
}

# (a + 2) - 3 -> a - 1
func &testConstVarFold12 (var %a i32) i32 {
  return (sub i32 (
    add i32 (dread i32 %a, constval i32 2),
    constval i32 3))
}

#  (a - 2) - 3 -> a - 5
func &testConstVarFold13 (var %a i32) i32 {
  return (sub i32 (
    sub i32 (dread i32 %a, constval i32 2),
    constval i32 3))
}

# 1 + (5 + ((5 - a) - 3)) -> 8 - a
func &testConstVarFold15 (var %a i32) i32 {
  return (add i32 (
    constval i32 1,
    add i32 (
      constval i32 5,
      sub i32 (
        sub i32 (constval i32 5, dread i32 %a),
        constval i32 3))))
}

# 1 + (5 + ((5 - a) * 3)) -> ((5 - a) * 3) + 6
func &testConstVarFold16 (var %a i32) i32 {
  return (add i32 (
    constval i32 1,
    add i32 (
      constval i32 5,
      mul i32 (
        sub i32 (constval i32 5, dread i32 %a),
        constval i32 3))))
}

# (-a) + b --> (-a) + b
func &testConstantFoldAdd00 (var %a i32, var %b i32) i32 {
  return (add i32 (
    neg i32 (dread i32 %a),
    dread i32 %b))
}

#  a + (-b) --> a + (-b)
func &testConstantFoldAdd01 (var %a i32, var %b i32) i32 {
  return (add i32 (
    dread i32 %a,
    neg i32 (dread i32 %b)))
}

# 5.0f + (5 - a) -> same, no folding with floating point numbers
func &testConstVarFold20 (var %a i32) i32 {
  return (add f32 (
    constval f32 5.0f,
    sub i32 (constval i32 5, dread i32 %a)))
}

func &testConstFold0 (var %a i32, var %b i32, var %c i32) i32 {
  # ((a + 5) + ((2 + c) + (4 + 2))) -> ((a + c) + 13)
  return (add i32 (
    add i32 (dread i32 %a, constval i32 0x5),
    add i32 (
      add i32 (constval i32 0x2, dread i32 %c),
      add i32 (constval i32 0x4, constval i32 0x2))))
}

func &testConstFold1 (var %a i32, var %b i32, var %c i32) i32 {
  # ((a + 5) + (2 + c)) -> ((a + c) + 7)
  return (add i32 (
    add i32 (dread i32 %a, constval i32 0x5),
    add i32 (constval i32 0x2, dread i32 %c)))
}

func &testConstFold2 (var %a i32, var %b i32, var %c i32) i32 {
  # ((a - 5) - ((b - c) - (4 - 2))) -> ((a - (b - c)) - 3)
  return (sub i32 (
    sub i32 (dread i32 %a, constval i32 0x5),
    sub i32 (
      sub i32 (dread i32 %b, dread i32 %c),
      sub i32 (constval i32 0x4, constval i32 0x2))))
}

func &testConstFold3 (var %a i32, var %b i32, var %c i32) i32 {
  # ((a - 5) + ((b - c) - (4 - 2))) -> ((a + (b - c)) - 7)
  return (add i32 (
    sub i32 (dread i32 %a, constval i32 0x5),
    sub i32 (
      sub i32 (dread i32 %b, dread i32 %c),
      sub i32 (constval i32 0x4, constval i32 0x2))))
}

func &testConstFold4 (var %a i32, var %b i32, var %c i32) i32 {
  # ((a - 5) - ((4 - 2) - (4 - b))) -> ((a - b) - 3)
  return (sub i32 (
    sub i32 (dread i32 %a, constval i32 0x5),
    sub i32 (
      sub i32 (constval i32 0x4, constval i32 0x2),
      sub i32 (constval i32 0x4, dread i32 %b))))
}

func &testConstFold5 (var %a i32, var %b i32, var %c i32) i32 {
  # ((5 - a) - ((b - 4) - (4 - b))) -> (((-a) - (b - (-b))) + 13)
  return (sub i32 (
    sub i32 (constval i32 0x5, dread i32 %a),
    sub i32 (
      sub i32 (dread i32 %b, constval i32 0x4),
      sub i32 (constval i32 0x4, dread i32 %b))))
}

func &testConstFold6 (var %a i32, var %b i32, var %c i32) i32 {
  # ((5 - a) - ((4 - b) - (b - 4))) -> (((-a) - ((-b) - b)) - 3)
  return (sub i32 (
    sub i32 (constval i32 0x5, dread i32 %a),
    sub i32 (
      sub i32 (constval i32 0x4, dread i32 %b),
      sub i32 (dread i32 %b, constval i32 0x4))))
}

func &testConstFold7 (var %a i32, var %b i32, var %c i32) i32 {
  # ((a - 5) - ((4 + 2) - (4 - (b - c)))) -> ((a - (b - c)) - 7)
  return (sub i32 (
    sub i32 (dread i32 %a, constval i32 0x5),
    sub i32 (
      add i32 (constval i32 0x4, constval i32 0x2),
      sub i32 (
        constval i32 0x4,
        sub i32 (dread i32 %b, dread i32 %c)))))
}

func &testConstFold8 (var %a i32, var %b i32, var %c i32) i32 {
  # ((a - 5) - ((4 + 2) - (4 * (b - c)))) -> ((a - (-(4 * (b - c)))) - 11)
  return (sub i32 (
    sub i32 (dread i32 %a, constval i32 0x5),
    sub i32 (
      add i32 (constval i32 0x4, constval i32 0x2),
      mul i32 (
        constval i32 0x4,
        sub i32 (dread i32 %b, dread i32 %c)))))
}

func &testConstFold9 (var %a i32, var %b i32, var %c i32) i32 {
  # ((a - 5) - ((4 + 2) - (-(4 - (b - c))))) -> ((a - (-(b - c))) - 15)
  return (sub i32 (
    sub i32 (dread i32 %a, constval i32 0x5),
    sub i32 (
      add i32 (constval i32 0x4, constval i32 0x2),
      neg i32 (sub i32 (
        constval i32 0x4,
        sub i32 (dread i32 %b, dread i32 %c))))))
}

func &testConstFold10 (var %a i32, var %b i32, var %c i32) i32 {
  # ((4 + 2) * (-(4 - (3 - c)))) -> (6 * ((-c) - 1))
  return (mul i32 (
    add i32 (constval i32 0x4, constval i32 0x2),
    neg i32 (sub i32 (
      constval i32 0x4,
      sub i32 (constval i32 0x3, dread i32 %c)))))
}

func &testConstFold11 (var %a i32, var %b i32, var %c i32) i32 {
  # ((4 + 2) * (4 - (3 - c))) -> (6 * (c + 1))
  return (mul i32 (
    add i32 (constval i32 0x4, constval i32 0x2),
    sub i32 (
      constval i32 0x4,
      sub i32 (constval i32 0x3, dread i32 %c))))
}

func &testConstFold12 (var %a i32, var %b i32, var %c i32) i32 {
  # ((4 + 2) * (~(4 - (3 - c)))) -> (6 * (~(c + 1)))
  return (mul i32 (
    add i32 (constval i32 0x4, constval i32 0x2),
    bnot i32 (sub i32 (
      constval i32 0x4,
      sub i32 (constval i32 0x3, dread i32 %c)))))
}

func &testConstFold13 (var %a i32, var %b i32, var %c i32) i32 {
  # (5 * (((4 - 2) + (2 & 4)) - (4 - (c | 3)))) -> (5 * ((c | 3) - 2))
  return (mul i32 (
    constval i32 0x5,
    sub i32 (
      add i32 (
        sub i32 (constval i32 0x4, constval i32 0x2),
        band i32 (constval i32 0x2, constval i32 0x4)),
      sub i32 (
        constval i32 0x4,
        bior i32 (dread i32 %c, constval i32 0x3)))))
}

func &testConstFold14 (var %a i32, var %b i32, var %c i32) i32 {
  # (5 * (((4 - 2) + (2 & 4)) * (4 - (c | 3)))) -> (5 * (2 * (4 - (c | 3))))
  return (mul i32 (
    constval i32 0x5,
    mul i32 (
      add i32 (
        sub i32 (constval i32 0x4, constval i32 0x2),
        band i32 (constval i32 0x2, constval i32 0x4)),
      sub i32 (
        constval i32 0x4,
        bior i32 (dread i32 %c, constval i32 0x3)))))
}

func &testConstFold15 (var %a i32, var %b i32, var %c i32) i32 {
  # ((1 + 0) * (~(4 - (0 * c)))) -> (~4)
  return (mul i32 (
    add i32 (constval i32 0x1, constval i32 0x0),
    bnot i32 (sub i32 (
      constval i32 0x4,
      mul i32 (constval i32 0x0, dread i32 %c)))))
}

func &testConstFold16 (var %a i32, var %b i32, var %c i32) i32 {
  # ((~(4 - (c * 0))) * (1 + 0)) -> (~4)
  return (mul i32 (
    bnot i32 (sub i32 (
      constval i32 0x4,
      mul i32 (dread i32 %c, constval i32 0x0))),
    add i32 (constval i32 0x1, constval i32 0x0)))
}

func &testConstFoldDiv0 (var %a i32, var %b i32, var %c i32) i32 {
  # (5 / 0) * (1 + 1) -> 2
  return (mul i32 (
    div i32 (constval i32 0x5, constval i32 0x0),
    add i32 (constval i32 0x1, constval i32 0x1)))
}

func &testConstFoldFloor00 (var %a i32, var %b i32, var %c i32) i32 {
  # (floor(4.25) + c) * (1 + 1) -> (c + 4) * 2
  return (mul i32 (
    add i32 (floor i32 f32 (constval f32 4.25f), dread i32 %c),
    add i32 (constval i32 0x1, constval i32 0x1)))
}

func &testConstFoldFloor01 (var %a i32, var %b i32, var %c i32) i32 {
  # (floor(4.25) + 2) * (1 + 1) -> 12
  return (mul i32 (
    add i32 (floor i32 f32 (constval f32 4.25f), constval i32 0x2),
    add i32 (constval i32 0x1, constval i32 0x1)))
}

func &testConstFoldFloor02 (var %a i32, var %b i32, var %c i32) i32 {
  # floor((floor(4.25) + 1) + (1 + a)) -> floor(a + 6)
  return (floor i32 i32 (add i32 (
    add i32 (floor i32 f32 (constval f32 4.25f), constval i32 0x1),
    add i32 (constval i32 0x1, dread i32 %a))))
}

func &testConstFoldExtractbitsNode00 () i32 {
  # 487 = 00111100111
  # extractbits i32 6 5 (487) = 7 = 0111
  return (extractbits i32 6 5 (constval i32 487))
}

func &testConstFoldExtractbitsNode01 () i32 {
  # 487 = 00111100111
  # extractbits i32 4 5 (487) = -2 = 0xfffffffffffffffe
  return (extractbits i32 4 5 (constval i32 487))
}

func &testConstFoldExtractbitsNode02 (var %a i32, var %b i32, var %c i32) i32 {
  # 487 = 00111100111
  # extractbits i32 4 5 (487) = -2 = 0xfffffffffffffffe
  # extractbits((extractbits(487) + 1) + (1 + a)) -> extractbits(a)
  return (extractbits i32 4 5 (add i32 (
    add i32 (extractbits i32 4 5 (constval i32 487), constval i32 0x1),
    add i32 (constval i32 0x1, dread i32 %a))))
}

func &testConstFoldExtractbitsNode03 (var %a i32, var %b i32, var %c i32) i32 {
  # 487 = 00111100111
  # extractbits i32 6 5 (487) = 7 = 0111
  # extractbits((extractbits(487) + 1) + (1 + a)) -> extractbits(a + 9)
  return (extractbits i32 6 5 (add i32 (
    add i32 (extractbits i32 6 5 (constval i32 487), constval i32 0x1),
    add i32 (constval i32 0x1, dread i32 %a))))
}

func &testConstFoldCompareNode00 (var %a i32, var %b i32, var %c i32) i32 {
  # (ge(6, 2) + 1) + (1 + 2) -> 5
  return (add i32 (
    add i32 (
      ge i32 i32 (constval i32 6, constval i32 2),
      constval i32 1),
    add i32 (constval i32 1, constval i32 2)))
}

func &testConstFoldCompareNode01 (var %a i32, var %b i32, var %c i32) i32 {
  # ge((ge(6, 2) + 1), (1 + a)) -> ge(2, a + 1)
  return (ge i32 i32 (
    add i32 (
      ge i32 i32 (constval i32 6, constval i32 2),
      constval i32 1),
    add i32 (constval i32 1, dread i32 %a)))
}

func &testConstFoldTernaryNode00 (var %a i32, var %b i32, var %c i32) i32 {
  # (6 > 2 ? 5 : 0) + (1 + 2) -> 8
  return ( add i32 (
    select i32 (
      ge i32 i32 (constval i32 6, constval i32 2),
      constval i32 5,
      constval i32 0),
    add i32 (constval i32 1,  constval i32 2)))
}

func &testConstFoldTernaryNode01 (var %a i32, var %b i32, var %c i32) i32 {
  # (6 > 10 ? 5 : 0) + (1 + a) -> (a + 1)
  return ( add i32 (
    select i32 (
      ge i32 i32 (constval i32 6, constval i32 10),
      constval i32 5,
      constval i32 0),
    add i32 (constval i32 1, dread i32 %a)))
}

func &testConstFoldTernaryNode02 (var %a i32, var %b i32, var %c i32) i32 {
  # 6 > a ? (1 + 6) + (1 + a) : 0 -> 6 > a ? a + 8 : 0
  return ( select i32 (
    ge i32 i32 (constval i32 6, dread i32 %a),
    add i32 (
      add i32 (constval i32 1, constval i32 6),
      add i32 (constval i32 1, dread i32 %a)),
    constval i32 0))
}

type $Hello <struct {
    @h i32,
    @e i32,
    @l i32,
    @o i32
  }>

type $Person <struct {
  @age i32,
  @extra <$Hello>,
  @sex u1}>

func &testConstFoldIassignNode00 (var %a i32) void {
  # this is equivalent to dassign %a 0 (constval i32 0x6)
  iassign <* i32> 0 (addrof ptr %a 0, add i32 (constval i32 1, constval i32 5))
}

func &testConstFoldIassignNode01 (var %a i32) void {
  var %bob $Person

  # this is equivalent to dassign %bob 5 (constval i32 0x6)
  iassign <* $Hello> 3 (addrof ptr %bob 2, add i32 (constval i32 1, constval i32 5))
}

func &testConstFoldIassignNode02 (var %a i32) void {
  # this is equivalent to dassign %a 0 (constval i32 0x6)
  iassign <* i32> 0 (
    iaddrof ptr <* i32> 0 (addrof ptr %a 0),
    add i32 (constval i32 1, constval i32 5))
}

func &testConstFoldIassignNode03 (var %a i32) void {

  var %bob $Person
  # this is equivalent to dassign %bob 5 (constval i32 0x6)
  iassign <* <$Hello>> 3 (
    iaddrof ptr <* <$Person>> 2 (addrof ptr %bob 0),
    add i32 (constval i32 1, constval i32 5))
}

func &testConstFoldIassignNode04 (var %a i32) void {

  var %bob $Person
  # this is equivalent to dassign %bob 5 (constval i32 0x6)
  iassign <* <$Hello>> 3 (
    iaddrof ptr <* <$Hello>> 0 (addrof ptr %bob 2),
    add i32 (constval i32 1, constval i32 5))
}

func &testConstFoldIassignNode05 (var %a i32) void {
  var %b <* i32>

  dassign %b (addrof ptr %a 0)

  # this is equivalent to iassign <* i32> 0 (dread ptr %a, constval i32 0x6)
  iassign <* i32> 0 (
    iaddrof ptr <* i32> 0 (dread ptr %a),
    add i32 (constval i32 1, constval i32 5))
}

func &testConstFoldIassignNode06 (var %a i32) void {
  var %bob $Person
  # this is equivalent to dassign %a 0 (dread i32 %bob 1)
  dassign %a 0 (iread agg <* <$Person>> 1 (addrof ptr %bob 0))
}

func &testConstFoldSwitchNode00 () void {
  # switch (1+4) ... -> goto @lab1
  switch ( add i32 (constval i32 1, constval i32 4)) @labdft {
  -2: goto @lab0
  5: goto @lab1
  8: goto @lab9 }

  @lab0
   return (constval i32 1)
  @labdft
   return (constval i32 2)
  @lab9
   return (constval i32 3)
  @lab1
   return (constval i32 4)
}

func &testConstFoldSwitchNode01 () void {
  # switch (1 + 10) ... -> goto @labdft
  switch ( add i32 (constval i32 1, constval i32 10)) @labdft {
  -2: goto @lab0
  5: goto @lab1
  8: goto @lab9 }

  @lab0
   return (constval i32 1)
  @labdft
   return (constval i32 2)
  @lab9
   return (constval i32 3)
  @lab1
   return (constval i32 4)
}

func &testConstFoldSwitchNode02 (var %a i32) void {
  # switch (5 + (5 - a)) ... -> switch (10 - a) ...
  switch (add i32 (
    constval i32 5,
    sub i32 (constval i32 5, dread i32 %a))) @labdft {
  -2: goto @lab0
  5: goto @lab1
  8: goto @lab9 }

  @lab0
   return (constval i32 1)
  @labdft
   return (constval i32 2)
  @lab9
   return (constval i32 3)
  @lab1
   return (constval i32 4)
}

func &testConstFoldArrayNode00 (var %x i32) void {
  var %a <* [15][15] i32>
  # fold indexes -> array ptr <* <[15][15] i32>> (constval i32 0x6, constval i32 0x4)
  dassign %x (
    iread i32 <* i32> (
      array 1 ptr <* [15][15] i32> (addrof ptr %a,
                                add i32 (constval i32 1, constval i32 5),
                                sub i32 (constval i32 5, constval i32 1))))
}

func &testDepositbitsNodeNode00 (var %a i32, var %c i32) i32 {
  return ( depositbits i32 1 23 (
             add i32 (
               add i32 (dread i32 %a, constval i32 0x5),
               add i32 (
                 add i32 (constval i32 0x2, dread i32 %c),
                 add i32 (constval i32 0x4, constval i32 0x2))),
             add i32 (constval i32 0x4, constval i32 0x8)))
}

func &testDepositbitsNodeNode01 () i32 {
  return ( depositbits i32 1 23 (
             add i32 (
               add i32 (constval i32 0x5, constval i32 0x5),
               add i32 (
                 add i32 (constval i32 0x2, constval i32 0x2),
                 add i32 (constval i32 0x4, constval i32 0x2))),
             add i32 (constval i32 0x4, constval i32 0x8)))
}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
sgn := proc(s) local i,j; signum(mul(mul(s[j] - s[i],j=i+1..nops(s)),i=1..nops(s))); end:

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

is_table_on := (A::set) -> proc(x)
 if not(type(x,table)) then return false; fi;
 if {indices(x)} <> map(a -> [a],A) then
  return false;
 fi;

 return true;
end:

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

# Return a randomly selected element of a finite set A

random_element_of := proc(A)
 local n,i;

 if not type(A,list) and not type(A,set) then
  error("A is not a list or set");
 fi;

 n := nops(A);

 if n = 0 then
  error("A is empty");
 fi;

 i := rand(1..n)();
 return A[i];
end:

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

# Return a randomly selected element of a finite set A

random_subset_of := proc(A,size_::nonnegint)
 local n,i,U,V,r,x;

 if not type(A,list) and not type(A,set) then
  error("A is not a list or set");
 fi;

 n := nops(A);

 if nargs = 1 then
  r := rand(2);
  V := select(a -> (r() = 1),{op(A)});
  return V;
 else
  if size_ > n then
   return FAIL;
  fi;
  U := {op(A)};
  V := {};
  for i from 0 to size_-1 do
   x := U[rand(1..(n-i))()];
   U := U minus {x};
   V := V union {x};
  od:
  return V;
 fi;
end:

random_nonempty_subset_of := proc(A)
 local n,k,B;

 if not type(A,list) and not type(A,set) then
  error("A is not a list or set");
 fi;

 n := nops(A);
 if n = 0 then error("A is empty"); fi;
 k := rand(1..n)();
 B := combinat[randcomb](A,k);
 return B;
end:

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

# This function converts small numbers to single-character strings,
# which is convenient in some circumstances when we need to display
# a large number of small numbers.
#
#  0 ..  9 become "0" to "9"
# 10 .. 35 become "A" to "Z"
# 36 .. 61 become "a" to "z"

nat_code := table([
 seq(i = convert([48+i],bytes),i=0..9),
 seq(i = convert([65+i-10],bytes),i=10..35),
 seq(i = convert([97+i-36],bytes),i=36..61)
]);

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

# Return a list of pairs [p[i],v[i]] such that the p[i] are distinct
# primes in increasing order, the v[i] are nonzero integers, and
# n is the product of the p[i]^v[i].  The real work is done by Maple's
# ifactor() function, we just convert the result to a more convenient
# form.

ifactor_list := proc(n::rational)
 local u;
 if n = 0 then
  error("Cannot factor zero");
 fi;
 u := ifactor(abs(n));
 if type(u,`*`) then
  u := [op(u)];
 else
  u := [u];
 fi;
 u := map(x -> `if`(type(x,`^`),[op(op(1,x)),op(2,x)],[op(x),1]),u);
 return u;
end:

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

# Return true iff n = p^k for some prime p and integer k >= 0.

is_prime_power := proc(n::posint)
 not(type(ifactor(n),`*`));
end:

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

# Return the p-adic valuation of m.

padic_val := proc(m::rational,p)
 local k,l;
 if m = 0 then return infinity; fi;
 if denom(m) <> 1 then
  return padic_val(numer(m),p) - padic_val(denom(m),p);
 fi;
 k := 0;
 l := m;
 while irem(l,p,'l') = 0 do k := k+1; od:
 return k;
end:

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

# Return the list of the base p digits of n, least significant first.

digit_list := proc(n::nonnegint,p::posint)
 local L,m;
 L := NULL;
 m := n;
 while m > 0 do
  L := L,irem(m,p);
  m := iquo(m,p);
 od:
 return [L];
end:

# Return the sum of the base p digits of n.

digit_sum := (n::nonnegint,p::posint) -> `+`(0,op(digit_list(n,p)));

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

# If u is an algebraic expression with an overall minus sign, then
# remove the mius sign.

strip_sign := proc(u)
 if type(u,list) or type(u,set) then
  return map(strip_sign,u);
 elif type(u,`*`) then
  return map(x -> `if`(type(x,numeric),abs(x),x),u);
 else
  return u;
 fi;
end:

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

# Given a Grobner basis such that the corresponding quotient ring is
# finite-dimensional, return the usual basis for that quotient.

cobasis := proc(basis,vars)
 local lm_basis,B,B0,V,b,i,v;

 lm_basis := map(u -> LeadingMonomial(u,vars), basis);
 B := [1];
 V := [op(vars)];
 for v in V do
  B0 := B;
  B := NULL;
  for b in B0 do
   i := 0;
   while NormalForm(b*v^i,lm_basis,vars) <> 0 do
    B := B,b*v^i;
    i := i+1;
   od:
  od:
  B := [B];
 od:
 return B;
end:

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

# Returns the set of monomials occuring in u

monomials := proc(u)
 if type(u,`+`) or type(u,list) or type(u,set) then
  return map(op,map(monomials,{op(u)}));
 elif type(u,`*`) then
  return {remove(type,u,constant)};
 elif type(u,constant) then
  return {1};
 else
  return {u};
 fi;
end:

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

# Returns the coefficient in u of the monomial m.
# If u involves variables that do not occur in m, then these will
# be included in the coefficient.

multi_coeff := proc(u,m)
 local c,x;
 c := u;
 for x in indets(m) do
  c := coeff(c,x,degree(m,x));
 od;
 return c;
end:

coeff_list := (u,mm) -> map2(multi_coeff,u,mm);

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

# This function returns a pair [d,q], where d is the gcd of the integers
# in u, and q is a list of integers such that the sum of the terms
# q[i] * u[i] is equal to d.  Moreover, q is the unique such list with
# 0 <= q[i] < d[i]/d[i-1] for i = 2,...,n, where d[i] is the gcd of
# u[1],...,u[i].  Note that there is no bound on q[1], which will
# typically be large and negative.

igcd_alt := proc(u::list(posint))
 local n,d,e,v,s,t,p,q,m,i;
 n := nops(u);
 d := table():
 e := table():
 v := table():
 d[1] := u[1];
 for i from 2 to nops(u) do
  d[i] := igcdex(d[i-1],u[i],s[i],t[i]);
  e[i-1] := d[i-1]/d[i];
  v[i] := u[i]/d[i];
 od:
 e[n] := d[n];

 p[n] := s[n];
 q[n] := t[n];
 m := iquo(q[n],e[n-1]);
 q[n] := q[n] - m*e[n-1];
 p[n] := p[n] + m*v[n];

 for i from n-1 to 2 by -1 do
  p[i] := s[i]*p[i+1];
  q[i] := t[i]*p[i+1];
  m := floor(q[i]/e[i-1]);
  q[i] := q[i] - m*e[i-1];
  p[i] := p[i] + m*v[i];
 od;

 q[1] := p[2];
 return [d[n],[seq(q[i],i=1..n)]];
end:

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

# Permutation matrix
# (duplicated in groups/general_linear.mpl)

perm_matrix := proc(s)
 local n,M,i;
 n := nops(s);
 M := Matrix(n,n);
 for i from 1 to n do
  M[s[i],i] := 1;
 od:
 return M;
end:

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

# Return u, but with all tiny numerical constants converted to 0.

trim := proc(u,epsilon := 10.^(-90))
 local c,v;
 if type(u,`+`) or type(u,list) or type(u,set) or type(u,Matrix) or type(u,Vector) then
  return map(trim,u,epsilon);
 elif type (u,`*`) then
  c,v := selectremove(type,u,complex(numeric));
  if abs(Re(c)) < epsilon then c := Im(c) * I; fi;
  if abs(Im(c)) < epsilon then c := Re(c); fi;
  return c * v;
 elif type(u,complex(numeric)) then
  c := u;
  if abs(Re(c)) < epsilon then c := Im(c) * I; fi;
  if abs(Im(c)) < epsilon then c := Re(c); fi;
  return c;
 else
  return u;
 fi;
end:

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

symmetric_difference := (A,B) ->
 {op(A),op(B)} minus ({op(A)} intersect {op(B)});

is_all_zero := proc(x)
 if type(x,{list,set}) then
  return `and`(op(map(is_all_zero,x)));
 elif type(x,Vector) then
  return is_all_zero(convert(x,list));
 elif type(x,Matrix) then
  return is_all_zero(convert(x,listlist));
 else
  return evalb(x = 0);
 fi;
end:

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

is_all_zero := proc(x)
 if type(x,{list,set}) then
  return `and`(op(map(is_all_zero,x)));
 elif type(x,Vector) then
  return is_all_zero(convert(x,list));
 elif type(x,Matrix) then
  return is_all_zero(convert(x,listlist));
 else
  return evalb(x = 0);
 fi;
end:

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

make_index := proc(L::{list,Vector})
 local L0,T,i;

 L0 := L;
 if type(L0,Vector) then L0 := convert(L0,list); fi;

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

 return eval(T);
end:

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

sum_terms := proc(x)
 if x = 0 then
  return [];
 elif type(x,`+`) then
  return [op(x)];
 else
  return [x];
 fi;
end;

coeff_split := proc(x,T := integer)
 if type(x,T) then
  return [x,1];
 elif type(x,`*`) then
  return [selectremove(type,x,T)];
 else
  return [1,x];
 fi;
end;


text:
# # Weyl group related operations on weyl word
# #       by Shizuo Kaji
# requires coxeter and weyl package by J.Stembridge
# http://www.math.lsa.umich.edu/~jrs/maple.html
#
# global variables:
# R: Lie type
# B: the set of simple roots

WeylOps := module() option package;
  export `.`,inv,descent,minrep,PReducedWd,WeylFunc,e,c,symmetric,e2w,e2a,a2e,a2w,w2a,deg_X,dimr,relative_roots,e2ref,action,oneline;
  # inverse of weyl word
  inv := proc (weyl::list);
	return ListTools[Reverse](weyl);
  end proc;

  # operator '.' as product
  `.` := overload(
  [
    proc(a::list, b::list) option overload(callseq_only); global R;
       return(reduce([op(a),op(b)],R));
    end,
    proc(a::indexed, b::indexed) option overload(callseq_only);
        if op(0,a)=`s` and op(0,b)=`s` then
            return(s[op(1,a).op(1,b)]);
        else
            return(a*b);
        end if;
    end,
    proc(a::`+`, b::anything)  option overload(callseq_only);
        return map(`.`,a,b);
    end,
    proc(a::anything, b::`+`)  option overload(callseq_only);
        return map2(`.`,a,b);
    end,
    proc(a::`*`, b::anything)  option overload(callseq_only);
        return op(1,a)*((a/op(1,a)).b);
    end,
    proc(a::anything, b::`*`)  option overload(callseq_only);
        return op(1,b)*(a.(b/op(1,b)));
    end
  ]
):

# descent set of the weyl word "weyl"
descent := proc(weyl::list) local i, V; global R;
	V := {};
	for i from 1 to rank(R) do
	   if nops(weyl.[i]) < nops(weyl) then
    		V := {op(V), i};
	   end if;
    end do:
	return V;
end proc;

# minimal coset representative of "weyl" according to "roots"
minrep := proc(weyl::list, roots::set) local i,temp; global R;
  for i in roots do
      temp:=weyl.[i]:
      if nops(temp) < nops(weyl) then return(procname(temp, roots)); end if:
  end do:
  return(weyl):
end proc:

# enumerate all reduced words (minimal coset representatives) of "R" upto length "len"
PReducedWd := proc (R, roots::set:={}, maxlen::integer:=nops(longest_elt(R)))
  local B, r, v, i, j, len, newvects, newwords, words, oldvects;
    B := base(R); r := rank(B);
	newvects := [add(weights(R)[i], i={$1..r} minus roots)];
    words := Array(0..maxlen); words[0]:=[[]];
    for len from 1 to maxlen do;
        oldvects := newvects;
        newvects := []; newwords := NULL;
        for j to nops(oldvects) do
            for i to r do
                if 0 < iprod(B[i], oldvects[j]) then
                    v := reflect(B[i], oldvects[j]);
                    if not member(v, newvects) then
                        newvects := [op(newvects), v];
                        newwords := newwords, reduce([i,op(words[len-1][j])],R);
                    end if
                end if
            end do
        end do;
        if newwords=NULL then break; end if;
        words[len] := [newwords];
    end do;
    return convert(words[1..len-1],list);
end proc;

# iterative application of a function of simple roots for a weyl word poly
WeylFunc := proc(weyl, f::polynom, func::procedure) local i,temp;
	temp := f;
    if type(weyl,'numeric') then
        return weyl*f;
    elif type(weyl,'list') then
        for i from nops(weyl) to 1 by -1 do
            temp := func(weyl[i], temp);
        end do;
        return temp;
    elif type(weyl, 'indexed') then
        if op(0,weyl)=`s` then
            return procname(op(1,weyl), f, func);
        else
            return weyl*f;
        end if;
    elif type(weyl, linear) then
	  return add(coeff(weyl,ind)*procname(ind, f, func), ind=indets(weyl));
    else
        print("error in WeylFunc:",weyl);
    end if;
end proc:

## Utility Functions
# coordinate variables
e := proc (j::integer);
    return cat('e', j);
end proc;

# elementary symmetric function "c(i)" of "t[1..r]"
c := proc (deg::integer,sym::symbol := `t`) global R;
    return symmetric(deg, [seq(sym[i], i = 1 .. rank(R))])
end proc;

# elementary symmetric function of degree "deg" in "vars"
symmetric := proc (deg::integer, vars::list) local T;
    return sort(expand(coeff(mul(T+i, i=vars), T, nops(vars)-deg)));
end proc;

# convert vector in "e" into weights "w[1..r]"
e2w := proc (el::linear) global R;
    return add(weight_coords(el, R)[i]*w[i], i=1..rank(R));
end proc;

# convert vector in "e" into simple roots "a[1..r]"
e2a := proc (el::{polynom, list}) global R;
    if type(el, rational) then
	  return el;
    elif type(el, linear) then
      return add(root_coords(el, R)[i]*a[i], i=1..rank(R));
    elif type(el, `*`) or type(el, list) then
	  return map(procname, el);
    elif type(el, `^`) then
	  return procname(op(1,el))^op(2,el);
    else
        print("error in e2a:",el);
        return infinity;
    end if;
end proc;

# convert vector in "e" into simple roots "a[1..r]"
a2e := proc (apol::{polynom, list}) global B;
    return eval(apol,[seq(a[i]=B[i],i=1..nops(B))]);
end proc:

# convert a-poly into w-poly
a2w := proc (apol::polynom) global R,B;
    return eval(apol, [seq(a[j]=add(weight_coords(B[j], R)[i]*w[i], i=1..nops(B)), j=1..nops(B))]);
end proc;

# convert w-poly into a-poly
w2a := proc (wpol::polynom) global R,B;
    return eval(wpol, [seq(w[j]=add(root_coords(weights(R)[j], R)[i]*a[i], i=1..nops(B)), j=1..nops(B))]);
end proc;

# degree of a polynomial in X[]
deg_X := proc(f::polynom) local d,max_d,term,temp,x,y;
    if type(f,`+`) then
        max_d:=0;
        for term in f do
            d := procname(term);
            if d>max_d then max_d:=d; end if;
        end do;
        return max_d;
    elif type(f,`^`) then
        procname(op(1,f))*op(2,f);
    else
        temp := f;
        for x in indets(f) do
            temp := eval(temp, x=y^(nops(op(1,x))));
        end do;
        return degree(temp);
    end if;
end proc:

# dimension of the standard representation of R
dimr := proc(R) local lie_type; option remember;
	lie_type := name_of(R);
    if lie_type = cat('A',rank(R)) then
        return rank(R)+1;
    elif lie_type = G2 then
        return 3;
    elif lie_type = E6 or lie_type = E7 then
        return 8;
    else
        return rank(R);
    end if;
end proc:

# positive roots relative to a parabolic subgroup
relative_roots := proc(R,roots::set) local L;
    L:={seq(base(R)[i],i=roots)};
    L:=orbit(L,R);
    return {op(pos_roots(R))} minus {op(L)};
end proc:

# vector to reflection
e2ref:=proc(el) local i,w; global R;
   vec2fc(reflect(el,interior_pt(R)),R,'w');
   return w;
end proc:

# action of a reduced word on a linear sum of ei
action:=proc(w::list,el::polynom) global B;
    return reflect(seq(B[i],i=w),el);
end proc:

# one line notation of elements in classical groups
oneline:=proc(weyl::list) local deg, lie_type, L, i; global R,pg;
	lie_type := name_of(R);
    if lie_type = cat('A',rank(R)) then
        deg := rank(R)+1;
    elif lie_type = cat('B',rank(R)) or lie_type = cat('C',rank(R)) or lie_type = cat('D',rank(R)) then
        deg := 2*rank(R);
    else
        return -infinity;
    end if;
    L:=convert(multperm(weyl,pg), 'permlist', deg);
    if lie_type = cat('B',rank(R)) or lie_type = cat('C',rank(R)) or lie_type = cat('D',rank(R)) then
        L:=L[-rank(R)..-1];
        for i from 1 to nops(L) do;
            if L[i]<=rank(R) then
                L[i]:=L[i]-rank(R)-1;
            else
                L[i]:=L[i]-rank(R);
            end if;
        end do;
    end if;
    return convert(L,array);
end proc:

end module:

max_stars_repo_name: shizuo-kaji/WeylGroup

max_issues_repo_name: shizuo-kaji/WeylGroup

max_forks_repo_name: shizuo-kaji/WeylGroup