text:
# Berechnungen für Systeme mit kinematischen Zwangsbedingungen
# Beschreibung
# Dieses Arbeitsblatt ist für Systeme mit kinematischen Zwangsbedingungen relevant.
# Es werden die Zusammenhänge zwischen den Minimalkoordinaten des Systems mit Zwangsbedingungen mit dem offenen System ohne Zwangsbedingungen hergestellt (nach [NakamuraGho1989], [ParkChoPlo1999]).
# Das System ohne Zwangsbedingungen kann durch die Erstellung einer neuen robot_env definiert werden (Weglassen der Zwangsbedingungen, offene Baumstruktur).
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# kinematic_constraints_calculations -> Berechnungen bezüglich kinematischer Zwangsbedingungen
# Quellen
# [KhalilKle1986] Khalil, W. & Kleinfinger, J.: A new geometric notation for open and closed-loop robots (1986)
# [NakamuraGho1989] Nakamura, Yoshihiko and Ghodoussi, Modjtaba: Dynamics computation of closed-link robot mechanisms with nonredundant and redundant actuators (1989)
# [ParkChoPlo1999] Park, FC and Choi, Jihyeon and Ploen, SR: Symbolic formulation of closed chain dynamics in independent coordinates (1999)
# Autor
# Moritz Schappler, schappler@irt.uni-hannover.de, 2017-12
# Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
interface(rtablesize=30):
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 2: # Hoher Optimierungsgrad.
;
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
# Definition und Zwangsbedingeun
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
if not assigned(simplify_options) or simplify_options(1)=-1 then
  if NJ <= 5 then
    use_simplify := 1: # Nur bei sehr einfachen Systemen (Viergelenkkette). Sonst hängt es sich auf.
  else
    use_simplify := 0: # Sonst keine Vereinfachung.
  end if:
else
  use_simplify := simplify_options(1): # erster Eintrag ist für Zwangsbedingungen
end if:
# Lesen der Zwangsbedingungen
kin_constraints_exist := false:
constrfile := sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
if FileTools[Exists](constrfile) then
  read constrfile:
end if:
if kin_constraints_exist = true then:
  kintmp_qs := kintmp_qs: # gelesene Variable sonst nicht sichtbar
  kintmp_qt := kintmp_qt: # gelesene Variable sonst nicht sichtbar
  kintmp_subsexp := kintmp_subsexp: # gelesene Variable sonst nicht sichtbar
  printf("%s. Kinematische Zwangsbedingungen gelesen.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
else
  printf("%s. Es gibt keine Zwangsbedingungen. Offene Struktur. Keine weiteren Berechnungen notwendig.", FormatTime("%Y-%m-%d %H:%M:%S")):
  restart: # Funktioniert nicht. Gedacht, damit im worksheet-Modus nichts mehr passieren kann.
  robot_name := "": # Damit werden die Export-Befehle ungültig
  codegen_act := false:
  quit:
end if:
# Abhängige Koordinaten definieren
# Nicht alle Ersetzungsausdrücke sind auch Koordinaten
NSE := RowDimension(kintmp_s):
# temporäre kinematische Ausdrücke (die auch die Zwangsbedingungen enthalten)
kintmp_s:
# Gelenkvariable für Drehgelenke (nach [KhalilKle1986])
theta_s := convert_t_s(theta):
# Gelenkvariable für Schubgelenke (nach [KhalilKle1986])
d_s := convert_t_s(d):
# Finde und zähle die Vorkommnisse von kintmp_s in theta_s und d_s. Die Anzahl der abhängigen Gelenkvariablen muss mit den definierten Minimalkoordinaten übereinstimmen, sonst ist das benutzerspezifizierte Modell falsch.
# Zähle, welche Zwangsbedingungen tatsächlich vorkommen
kintmp_zaehler := Matrix(NSE,1):
for i from 1 to NSE do
  if has(has~(<theta_s(1..NJ,1);d_s(1..NJ,1)>, kintmp_s(i,1)), true) then
    kintmp_zaehler(i,1) := 1:
  end if:
end do:
# Zähle, welche Gelenke abhängig sind (durch ZB beschrieben).
# (in der Gelenkvariablen stehen die Inhalte von kintmp)
jv_zaehler := Matrix(NJ,1):
for i from 1 to NJ do
  if has(has~(theta_s(i,1), kintmp_s), true) or has(has~(d_s(i,1), kintmp_s), true) then
    jv_zaehler(i,1) := 1:
  end if:
end do:
# Zählen, welche Gelenke statische Transformationen sind
sigma2_zaehler := Matrix(NJ,1):
for i from 1 to NJ do
  if sigma(i,1)=2 then
    sigma2_zaehler(i,1) := 1:
  end if:
end do:
# Anzahl der abhängigen Gelenkvariablen (Zählmethode: Gelenkkoordinaten durchgehen)
NC1 := add(jv_zaehler(n), n = 1 .. NJ):
# Anzahl der abhängigen Gelenkvariablen (Zählmethode: Zwangsbedingungen durchgehen)
NC2 := add(kintmp_zaehler(n), n = 1 .. NSE):
# Anzahl der konstanten Transformationen
NS := add(sigma2_zaehler(n), n = 1 .. NJ):
# Anzahl der "realen" Gelenke (statische Transformationen herausrechnen)
NRJ := NJ - NS:
# Es gibt (NL - 1) durch Gelenke bewegte Körper, NVJ/2 Koordinatensysteme (hier als "Gelenk" gezählt) sind statisch aufgrund der Schleifenschlussbedingungen
if NQJ+NC1 <> NRJ then
  printf("%s. Fehler: Anzahl der Abhängigen Gelenke %d (aus Gelenkvariablen) passt nicht zur Anzahl der Gelenke %d und der Minimalkoordinaten %d\n", FormatTime("%Y-%m-%d %H:%M:%S"), NC1, NRJ, NQJ):
end if:
if NC1 <> NC2 then
  printf("%s. Fehler: Anzahl der Abhängigen Gelenke %d (aus Gelenkvariablen) passt nicht zur Anzahl der Zwangsbedingungen %d (aus Verwendung der Zwangsbedingungs-Variablen)\n", FormatTime("%Y-%m-%d %H:%M:%S"), NC1, NC2):
end if:
# Indizes der Zwangsbedingungen (in kintmp_s), die wirklich verwendet werden. Wird später nicht mehr benötigt.
Ind_ZB := Matrix(NC,1):
j := 1:
for i from 1 to NSE do
  if kintmp_zaehler(i,1) = 1 then
    Ind_ZB(j,1) := i:
    j := j + 1:
  end if:
end do:
# Koordinaten der offenen Struktur (MDH-Gelenkvariablen)
# Darstellung der Gelenkvariablen in abhängigkeit der Minimalkoordinaten. Diese Gelenkvariablen sind die Minimalkoordinaten des offenen Systems (ohne Zwangsbedingungen)
jv_s := theta_s*~(Matrix(NJ,1,1)-sigma) + d_s*~sigma - qoffset:
jv_qs := jv_s: # Initialisierung
jv_qt := Matrix(NJ,1): # Initialisierung
;
# Ersetzen der abhängigen Gelenkkoordinaten
for i from 1 to NJ do # Index über alle Gelenkvariablen
  # Substituiere sin und cos der Winkel (einfachere Ausdrücke. Sollten aber hier nicht drin sein.)
  jv_qs(i,1) := subs_kintmp_exp(jv_qs(i,1)):
  # Substituiere die verbleibenden Winkel direkt
  for jj from 1 to NSE do # Index über zu ersetzende Winkel
    jv_qs(i,1) := subs( { kintmp_s(jj, 1) = kintmp_qs(jj, 1) }, jv_qs(i,1) ):
  end do:
  jv_qt(i,1) := convert_s_t( jv_qs(i,1) ):
end do:
# Jacobi-Matrix der abhängigen Winkel
# Ableitung der Gelenkkoordinaten nach den Minimalkoordinaten.
# Wird benötigt, um die Geschwindigkeit der Gelenkkoordinaten der offenen Struktur zu berechnen und um die Momente der offenen Struktur auf die geschlossene Struktur umzurechnen.
# Ist Matrix Phi aus [ParkChoPlo1999]. qJ hier ist qa dort. qr dort entspricht den Koordinaten desselben Systems ohne Zwangsbedingungen.
# Entspricht der Matrix W aus [NakamuraGho1989] (nur dass hier die Aufteilung q1,q2 nicht gestapelt in q auftreten muss).
# Wenn statische Transformationen definiert sind (mit sigma=2), dann werden die entsprechenden Spalten hier trotzdem angegeben (das ist nicht unbedingt intuitiv, aber von der Programmierung bedeutend einfach)
Phi_s := Matrix(NJ, NQJ):
# Jacobi-Matrix berechnen
for i from 1 to NJ do
  for j from 1 to NQJ do
    Phi_s(i,j) := diff(jv_qs(i,1), qJ_s(j,1)):
  end do:
end do:
# Term vereinfachen
if use_simplify >= 1 then
  tmp_t1 := time(): tmp_l1 := length(Phi_s):
  printf("%s. Vereinfache Passiv-Gelenk-Jacobi (Dimension %d x %d). Länge: %d.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), NJ, NQJ, tmp_l1):
  for i from 1 to NJ do
    for j from 1 to NQJ do
      tmp_t1ij := time(): tmp_l1ij := length(Phi_s(i,j)):
      printf("%s. Vereinfache Passiv-Gelenk-Jacobi Eintrag %d,%d. Länge: %d.\n", \
        FormatTime("%Y-%m-%d %H:%M:%S"), i, j, tmp_l1ij):
      Phi_s(i,j) := simplify2(Phi_s(i,j)):
      tmp_t2ij := time(): tmp_l2ij := length(Phi_s(i,j)):
      printf("%s. Passiv-Gelenk-Jacobi der ZB Eintrag %d,%d vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
        FormatTime("%Y-%m-%d %H:%M:%S"), i, j, tmp_l1ij, tmp_l2ij, tmp_t2ij-tmp_t1ij):
    end do:
  end do:
  tmp_t2 := time(): tmp_l2 := length(Phi_s):
  printf("%s. Passiv-Gelenk-Jacobi der ZB vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t1):
end if:
# Zeitableitung der Jacobi-Matrix
# Wird benötigt, um die Beschleunigung der Gelenkkoordinaten der offenen Struktur zu berechnen.
Phi_t := convert_s_t(Phi_s):
PhiD_t := diff~(Phi_t, t):
PhiD_s := convert_t_s(PhiD_t):
# Exportiere Code für folgende Skripte
# Speichere Maple-Ausdruck (Eingabe-Format und internes Format)
save jv_s, jv_qs, jv_qt, sprintf("../codeexport/%s/tmp/kinematic_constraints_explicit_maple.m", robot_name):
save Phi_s, sprintf("../codeexport/%s/tmp/kinematic_constraints_explicit_jacobian_maple.m", robot_name):
save PhiD_s, sprintf("../codeexport/%s/tmp/kinematic_constraints_explicit_jacobian_time_derivative_maple.m", robot_name):
printf("%s. Ausdrücke für Kinematische ZB gespeichert (Maple)\n", FormatTime("%Y-%m-%d %H:%M:%S")):
if codegen_act = true then
  MatlabExport(jv_zaehler, sprintf("../codeexport/%s/tmp/kinconstr_index_dependant_joints_matlab.m", robot_name), 2):
  MatlabExport(jv_qs, sprintf("../codeexport/%s/tmp/kinconstr_expl_matlab.m", robot_name), 2):
  MatlabExport(Phi_s, sprintf("../codeexport/%s/tmp/kinconstr_expl_jacobian_matlab.m", robot_name), 2):
  MatlabExport(PhiD_s, sprintf("../codeexport/%s/tmp/kinconstr_expl_jacobianD_matlab.m", robot_name), 2):
  printf("%s. Ausdrücke für Kinematische ZB gespeichert (Matlab)\n", FormatTime("%Y-%m-%d %H:%M:%S")):
end if:

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


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

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

`gamma/stasheff_star` := (A::set,B::set) -> (p) -> proc(U,V)
 local R,S,T,r,s,t,b,n,JJ,J,pJ,J1,i,j,k;

 S,s := op(U);
 R := table();
 r := table();
 for b in B do
  R[b],r[b] := op(V[b]);
 od;

 T := [seq(op(R[b]),b in S)];

 t := table();

 n := nops(T);
 JJ := {seq(seq({seq(T[k],k=i..j)},j=i..n),i=1..n)};
 for J in JJ do
  pJ := map(j -> p[j],J);
  J1 := select(j -> member(p[j],pJ),A);
  if J = J1 then
   t[J] := s[pJ];
  elif nops(pJ) = 1 then
   t[J] := r[op(pJ)][J];
  else
   t[J] := 0;
  fi;
 od:

 return [T,eval(t)];
end;

ext: mpl


text:
gfun:=proc(x)
    uses Grid;
    global y,z,w,u,v;
    print(x,y,z,w,u,v);
    Barrier();
end proc:
Foo:=module()
    option package;
    local x:=1,z:=3;
    export fun;
    fun:=proc()
        local y:=2,w:=4,u:=5,v:=6;
        Grid:-Launch(gfun,x,imports=[':-y'=y,':-z'=z,'w'=w,'u',"v"]);
    end proc:
end module:
Foo:-fun();
y;z;


text:
#@ Not autoload

MU_n_max := 10;

make_fgl_Lazard := proc(d)
 global a,m,ma,am;
 local e,T,u,x,y,q,i,j,k,Fm,Fmp,err,sol;

 T := table();
 T["degree"] := d;
 Order := d+2;

 T["log_m"] := unapply(x + add(m[i]*x^(i+1),i=1..d-1),x);
 T["exp_m"] :=
   unapply(
   convert(
    solve(x = series(T["log_m"](y) + sin(y)^(d+1),y=0,d+1),y),
    polynom,x
   ),
   x
  );

 Fm := x+y;
 for k from 2 to d do
  Fmp[0] := 1;
  Fmp[1] := subs({x=e*x,y=e*y},Fm);
  for j from 2 to k do
   Fmp[j] := rem(expand(Fmp[1] * Fmp[j-1]),e^(k+1),e);
  od;
  err := expand(Fmp[1] - e*x - e*y +
   add(m[i]*(Fmp[i+1] - e^(i+1)*(x^(i+1)+y^(i+1))),i=1..k-1));
  Fm := Fm - subs(e=1,err);
 od:

 T["sum_m"] := unapply(Fm,x,y):

 am[1] := m[1];
 for k from 2 to d-1 do
  q := igcd_alt([seq(binomial(k+1,i),i=1..k)])[2];
  am[k] := -add(q[i]*coeff(coeff(Fm,x,i),y,k+1-i),i=1..k);
 od:
 sol := solve({seq(a[i]=am[i],i=1..d-1)},{seq(m[i],i=1..d-1)}):
 for i from 1 to d-1 do ma[i] := subs(sol,m[i]); od:

 T["log_a"] := unapply(collect(expand(eval(subs(m = ma,T["log_m"](x)))),x),x);
 T["exp_a"] := unapply(collect(expand(eval(subs(m = ma,T["exp_m"](x)))),x),x);
 T["sum_a"] := unapply(collect(expand(eval(subs(m = ma,T["sum_m"](x,y)))),{x,y}),x,y);

 return eval(T);
end:

MU_vars := plex(seq(a[i],i=1..MU_n_max));
MUMU_vars := plex(seq(a[i],i=1..MU_n_max),seq(b[i],i=1..MU_n_max));


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

hcth_a_beta  :=  0.0042:
hcth_a_gamma :=  6:
hcth_a_c0    :=  1.09878:
hcth_a_c1    := -2.51173:
hcth_a_c2    :=  0.0156233:

hcth_a_aux := x -> 1 + hcth_a_gamma*hcth_a_beta*x*arcsinh(x):

hcth_a_f := x -> hcth_a_c0 + hcth_a_beta/X_FACTOR_C*x^2*(hcth_a_c1/hcth_a_aux(x)
     + hcth_a_c2/(hcth_a_beta*hcth_a_aux(x)^2)):

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


text:
function compare(a, b)
    if islist(a) and islist(b) then
        if size(a) != size(b) then
	    return 0
	end
	for i in 0:size(a)-1
	    if compare(a[i], b[i]) == 0 then
	        return 0
	    end
	    i = i + 1
	end
	return 1
    end
    return a == b
end

function error(a, b, s)
    "error " + s
    a
    b
    a-b
    exit(1)
end

function test(a, b, s)
    if not compare(a, b) then
        error(a, b, s)
    end
end

function test_addition()
    a = { 5, 1, (2,3), (1,2), [1,2,3], [1,2], [(1,2),(3,4)], \
          [(1,2),(3,4)], \
	  [1,2|3,4], [1,2|3,4], [(1,2),(3,4)|(5,6),(7,8)], \
	  [(1,2),(3,4)|(5,6),(7,8)], "abc" }
    b = { 3, (2,3), 1, (3,4), [4,5,6], [(1,2),(3,4)], [1,2], [(5,6),(7,8)], \
          [5,6|7,8], [(1,2),(3,4)|(5,6),(7,8)], [1,2|3,4], \
	  [(8,7),(6,5)|(4,3),(2,1)], "def" }
    q = { 8, (3,3), (3,3), (4,6), [5,7,9], [(2,2),(5,4)], [(2,2),(5,4)], [(6,8),(10,12)], \
          [6,8|10,12], [(2,2),(5,4)|(8,6),(11,8)], [(2,2),(5,4)|(8,6),(11,8)], \
	  [(9,9),(9,9)|(9,9),(9,9)], "abcdef" }
    s = { "real + real", "real + complex", "complex + real", \
          "complex + complex", "vector + vector", "vector + cvector", \
	  "cvector + vector", "cvector + cvector", \
	  "matrix + matrix", "matrix + cmatrix", "cmatrix + matrix", \
	  "cmatrix + cmatrix", "string + string" }
    c = a + b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
        i = i + 1
    end
end

function test_subtraction()
    a = { 5, 1, (2,3), (1,2), [4,5,6], [1,2], [(1,2),(3,4)], [(1,2),(3,4)], \
          [5,6|7,8], [1,2|3,4], [(1,2),(3,4)|(5,6),(7,8)], \
	  [(8,7),(6,5)|(4,3),(2,1)], "abc" }
    b = { 3, (2,3), 1, (3,4), [1,2,3], [(1,2),(3,4)], [1,2], [(5,6),(7,8)], \
          [1,2|3,4], [(1,2),(3,4)|(5,6),(7,8)], [1,2|3,4], \
	  [(1,2),(3,4)|(5,6),(7,8)], "def" }
    x = { 2, (-1,-3), (1,3), (-2,-2), [3,3,3], [(0,-2),(-1,-4)], \
          [(0,2),(1,4)], [(-4,-4),(-4,-4)], \
	  [4,4|4,4], [(0,-2),(-1,-4)|(-2,-6),(-3,-8)], \
	  [(0,2),(1,4)|(2,6),(3,8)], [(7,5),(3,1)|(-1,-3),(-5,-7)], "abcfed" }
    s = { "real - real", "real - complex", "complex - real", \
          "complex - complex", "vector - vector", "vector - cvector", \
	  "cvector - vector", "cvector - cvector", \
	  "matrix - matrix", "matrix - cmatrix", "cmatrix - matrix", \
	  "cmatrix - cmatrix", "string - string" }
    c = a - b
    for i in 0:size(a)-1
        test(x[i], c[i], s[i])
        i = i + 1
    end
end

function test_multiplication()
    a = { 5,                  2,                  (2,3),            (1,2), \
          2,                  [1,2,3],            [1,2,3], \
	  2,                  (1,2),              (1,2), \
	  [(1,2),(3,4)],      [1,2],              [(1,2),(3,4)], \
	  [1,2],              [(1,2),(3,4)],      [(1,2),(3,4)], \
	  2,                  [1,2|3,4],          [1,2|3,4], \
	  [1,2|3,4],          (1,2),              2, \
	  (1,2),              [1,2|3,4],          [(1,2),(3,4)|(5,6),(7,8)], \
	  [(1,2),(3,4)|(5,6),(7,8)],              [(1,2),(3,4)|(5,6),(7,8)], \
	  [1,2|3,4],          [(1,2),(3,4)|(5,6),(7,8)], \
	  [(1,2),(3,4)|(5,6),(7,8)],              [1,2|3,4], \
	  [(1,2),(3,4)|(5,6),(7,8)], \
	  2,                  "abc" }

    b = { 3,                  (2,3),              2,                (3,4), \
          [1,2,3],            2,                  [2,3,4], \
          [(1,2),(3,4)],      [1,2],              [(1,2),(3,4)], \
	  2,                  (1,2),              (1,2), \
	  [(1,2),(3,4)],      [1,2],              [(4,3),(2,1)], \
	  [1,2|3,4],          2,                  [1,2], \
	  [1,3|2,4],          [1,2|3,4],          [(1,2),(3,4)|(5,6),(7,8)], \
	  [(1,2),(3,4)|(5,6),(7,8)], (1,2),       [1,2], \
	  2,                  (1,2), \
	  [(1,2),(3,4)],      [(1,2),(3,4)], \
	  [1,2|3,4],          [(1,2),(3,4)|(5,6),(7,8)], \
	  [(8,7),(6,5)|(4,3),(2,1)], \
	  "abc",              2 }

    q = { 15,                 (4,6),              (4,6),            (-5,10), \
          [2,4,6],            [2,4,6],            20, \
          [(2,4),(6,8)],      [(1,2),(2,4)],      [(-3,4),(-5,10)], \
	  [(2,4),(6,8)],      [(1,2),(2,4)],      [(-3,4),(-5,10)], \
	  (7,-10),            (7,10),             (20,10), \
	  [2,4|6,8],          [2,4|6,8],          [5,11], \
	  [5,11|11,25],       [(1,2),(2,4)|(3,6),(4,8)], [(2,4),(6,8)|(10,12),(14,16)], \
	  [(-3,4),(-5,10)|(-7,16),(-9,22)], [(1,2),(2,4)|(3,6),(4,8)], [(7,10),(19,22)], \
	  [(2,4),(6,8)|(10,12),(14,16)], [(-3,4),(-5,10)|(-7,16),(-9,22)], \
	  [(7,10),(15,22)],   [(-10,28),(-18,68)], \
	  [(10,14),(14,20)|(26,30),(38,44)], [(11,14),(17,20)|(23,30),(37,44)], \
	  [(-6,48),(-2,28)|(2,136),(6,84)], \
	  "abcabc",           "abcabc" }

    s = { "real * real",      "real * complex",   "complex * real", "complex * complex", \
          "real * vector",    "vector * real",    "vector * vector", \
	  "real * cvector",   "complex * vector", "complex * cvector", \
	  "cvector * real",   "vector * complex", "cvector * complex", \
	  "vector * cvector", "cvector * vector", "cvector * cvector", \
	  "real * matrix",    "matrix * real",    "matrix * vector", \
	  "matrix * matrix",  "complex * matrix", "real * cmatrix", \
	  "complex * cmatrix", "matrix * complex", "matrix * complex", \
	  "cmatrix * real",   "cmatrix * complex", \
	  "matrix * cvector", "cmatrix * cvector", \
	  "cmatrix * matrix", "matrix * cmatrix", \
	  "cmatrix * cmatrix", \
	  "real * string",    "string * real" }
    c = a*b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
        i = i + 1
    end
end

function test_pointwise_multiplication()
    a = { [1,2,3], [1,2,3], [(1,2),(3,4),(5,6)], [(1,2),(3,4)], \
          [1,2|3,4], [1,2|3,4], [(1,2),(3,4)|(5,6),(7,8)], \
	  [(1,2),(3,4)|(5,6),(7,8)] }

    b = { [4,5,6], [(4,5),(6,7),(8,9)],[7,8,9], [(5,6),(7,8)], \
          [5,6|7,8], [(1,2),(3,4)|(5,6),(7,8)], [1,2|3,4], \
	  [(1,2),(3,4)|(5,6),(7,8)] }

    q = { [4,10,18], [(4,5),(12,14),(24,27)], \
          [(7,14),(24,32),(45,54)], [(-7,16),(-11,52)], \
	  [5,12|21,32], [(1,2),(6,8)|(15,18),(28,32)], \
	  [(1,2),(6,8)|(15,18),(28,32)], \
	  [(-3,4),(-7,24)|(-11,60),(-15,112)] }

    s = { "vector .* vector", "vector .* cvector", "cvector .* vector", \
          "cvector .* cvector", "matrix .* matrix", "matrix .* cmatrix", \
	  "cmatrix .* matrix", "cmatrix .* cmatrix" }

    c = a.*b

    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
	i = i + 1
    end
end

function test_division()
    a = { 6, 13, (4,6), (-5,10), [2,4,6], [1,2], [(1,2),(3,4)], \
          [(1,2),(3,4)], [2,4|6,8], [1,2|3,4], \
	  [(1,2),(3,4)|(5,6),(7,8)], [(1,2),(3,4)|(5,6),(7,8)] }
    b = { 2, (2,-3), 2, (3,4), 2, (1,2), 2, (1,2), 2, (1,2), 2, (1,2) }
    q = { 3, (2,3), (2,3), (1,2), [1,2,3], [(0.2,-0.4),(0.4,-0.8)], \
          [(0.5,1),(1.5,2)], [(1,0),(2.2,-0.4)], [1,2|3,4], \
	  [(0.2,-0.4),(0.4,-0.8)|(0.6,-1.2),(0.8,-1.6)], \
	  [(0.5,1),(1.5,2)|(2.5,3),(3.5,4)], \
	  [(1,0),(2.2,-0.4)|(3.4,-0.8),(4.6,-1.2)] }
    s = { "real / real", "real / complex", "complex / real", \
          "complex / complex", "vector / real", "vector / complex", \
	  "cvector / real", "cvector / complex", "matrix / real", \
	  "matrix / complex", "cmatrix / real", "cmatrix / complex" }
    c = a/b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
        i = i + 1
    end
end

function test_pointwise_division()
    a = { [4,10,18], [5,25,100], \
          [(7,14),(24,32),(45,54)], [(11,52),(-11,52)], \
	  [5,12|21,32], [5,25|100,113], \
	  [(1,2),(6,8)|(15,18),(28,32)], \
	  [(-3,4),(-7,24)|(11,52),(-15,112)] }

    b = { [4,5,6], [(1,2),(3,4),(6,8)], [7,8,9], [(8,7),(7,8)], \
          [5,6|7,8], [(1,2),(3,4)|(6,8),(7,8)], [1,2|3,4], \
	  [(1,2),(3,4)|(8,7),(7,8)] }

    q = { [1,2,3], [(1,-2),(3,-4),(6,-8)], [(1,2),(3,4),(5,6)], [(4,3),(3,4)], \
          [1,2|3,4], [(1,-2),(3,-4)|(6,-8),(7,-8)], [(1,2),(3,4)|(5,6),(7,8)], \
	  [(1,2),(3,4)|(4,3),(7,8)] }

    c = a./b

     s = { "vector ./ vector", "vector ./ cvector", "cvector ./ vector", \
          "cvector ./ cvector", "matrix ./ matrix", "matrix ./ cmatrix", \
	  "cmatrix ./ matrix", "cmatrix ./ cmatrix" }

    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
	i = i + 1
    end
end

function test_floored_division()
    test(2, 8 div 3, "real div real")
end

function test_modulo()
    test(2, 8 mod 3, "real mod real")
end

function test_exponent()
    a = { 3, 3, (3,0), (3,0) }
    b = { 2, (2,0), 2, (2,0) }
    q = { 9, (9,0), (9,0), (9.00000000000000177636,0) }
    s = { "real ^ real", "real ^ complex", "complex ^ real", \
          "complex ^ complex" }
    c = a^b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
        i = i + 1
    end
end

function test_negation()
    a = {3, (1,2), [1,2,3], [(1,2),(3,4)], [1,2|3,4], \
         [(1,2),(3,4)|(5,6),(7,8)], "abc"}
    q = {-3, (-1,-2), [-1,-2,-3], [(-1,-2),(-3,-4)], [-1,-2|-3,-4], \
         [(-1,-2),(-3,-4)|(-5,-6),(-7,-8)], "cba" }
    s = { "- real", "- complex", "- vector", "- cvector", "- matrix", \
          "- cmatrix", "- string" }
    c = -a
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
	i = i + 1
    end
end

function test_arithmetic()
    test_addition()
    test_subtraction()
    test_multiplication()
    test_pointwise_multiplication()
    test_division()
    test_pointwise_division()
    test_floored_division()
    test_modulo()
    test_exponent()
    test_negation()
end

function test_or()
    a = {1, 1, 0, 0}
    b = {1, 0, 1, 0}
    q = {1, 1, 1, 0}
    c = a or b
    for i in 0:size(a)-1
        test(q[i], c[i], "or")
	i = i + 1
    end
end

function test_and()
    a = {1, 1, 0, 0}
    b = {1, 0, 1, 0}
    q = {1, 0, 0, 0}
    c = a and b
    for i in 0:size(a)-1
        test(q[i], c[i], "and")
	i = i + 1
    end
end

function test_not()
    a = { 0, 1 }
    q = { 1, 0 }
    c = not a
    for i in 0:size(a)-1
        test(q[i], c[i], "not")
	i = i + 1
    end
end

function test_bor()
    test(14, 10 bor 12, "bor")
end

function test_bxor()
    test(6, 10 bxor 12, "bxor")
end

function test_band()
    test(8, 10 band 12, "band")
end

function test_bnot()
    test(-2, bnot 1, "bnot")
end

function test_logic()
    test_or()
    test_and()
    test_not()
    test_bor()
    test_bxor()
    test_band()
    test_bnot()
end

function test_eq()
    a = { 3, 1, (1,2), (1,2), [1,2,3], [1,2], [(1,0),(2,0)], \
          [(1,2),(3,4)], [1,2|3,4], [1,2|3,4], [(1,0),(2,0)|(3,0),(4,0)], \
	  [(1,2),(3,4)|(5,6),(7,8)], "abc" }
    b = { 3, (2,3), 3, (3,4), [1,2,3], [(1,0),(2,0)], [1,2], \
          [(1,2),(3,4)], [1,2|3,4], [(1,0),(2,0)|(3,0),(4,0)], [1,2|3,4], \
	  [(1,2),(3,4)|(5,6),(7,8)], "def" }
    q = { 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0 }
    s = { "real == real", "real == complex", "complex == real", \
           "complex == complex", "vector == vector", "vector == cvector", \
	   "cvector == vector", "cvector == cvector", "matrix == matrix", \
	   "matrix == cmatrix", "cmatrix == matrix", "cmatrix == cmatrix", "string == string" }
    c = a == b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
	i = i + 1
    end
end

function test_ne()
    a = { 3, 1, (1,2), (1,2), [1,2,3], [1,2], [(1,0),(2,0)], \
          [(1,2),(3,4)], [1,2|3,4], [1,2|3,4], [(1,0),(2,0)|(3,0),(4,0)], \
	  [(1,2),(3,4)|(5,6),(7,8)], "abc" }
    b = { 3, (2,3), 3, (3,4), [1,2,3], [(1,0),(2,0)], [1,2], \
          [(1,2),(3,4)], [1,2|3,4], [(1,0),(2,0)|(3,0),(4,0)], [1,2|3,4], \
	  [(1,2),(3,4)|(5,6),(7,8)], "def" }
    q = { 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1 }
    s = { "real != real", "real != complex", "complex != real", \
           "complex != complex", "vector != vector", "vector != cvector", \
	   "cvector != vector", "cvector != cvector", "matrix != matrix", \
	   "matrix == cmatrix", "cmatrix == matrix", "cmatrix == cmatrix", "string != string" }
    c = a != b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
	i = i + 1
    end
end

function test_lt()
    a = { 3, "abc" }
    b = { 3, "def" }
    q = { 0, 1 }
    s = { "real < real", "string < string" }
    c = a < b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
	i = i + 1
    end
end

function test_le()
    a = { 3, "abc" }
    b = { 3, "def" }
    q = { 1, 1 }
    s = { "real <= real", "string <= string" }
    c = a <= b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
	i = i + 1
    end
end

function test_gt()
    a = { 3, "abc" }
    b = { 3, "def" }
    q = { 0, 0 }
    s = { "real > real", "string > string" }
    c = a > b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
	i = i + 1
    end
end

function test_ge()
    a = { 3, "abc" }
    b = { 3, "def" }
    q = { 1, 0 }
    s = { "real >= real", "string >= string" }
    c = a >= b
    for i in 0:size(a)-1
        test(q[i], c[i], s[i])
	i = i + 1
    end
end

function test_relational_operators()
    test_eq()
    test_ne()
    test_lt()
    test_le()
    test_gt()
    test_ge()
end

function test_vector_index()
    a = [1, 2, 3, 4]
    test(2, a[1], "vector[i]")
    test([2, 3], a[1:2], "vector[i:i]")
    test([2, 4], a[1:2:2], "vector[i:i:i]")
    a[1] = 20
    test([1, 20, 3, 4], a, "vector[i] = real")
    a[1:2] = [200, 300]
    test([1, 200, 300, 4], a, "vector[i:i] = vector")
    a[1:2:2] = [2000, 4000]
    test([1, 2000, 300, 4000], a, "vector[i:i:i] = vector")
end

function test_cvector_index()
    a = [(1,2),(3,4),(5,6),(7,8)]
    test((3,4), a[1], "cvector[i]")
    test([(3,4),(5,6)], a[1:2], "cvector[i:i]")
    test([(3,4),(7,8)], a[1:2:2], "cvector[i:i:i]")
    a[1] = (30,40)
    test([(1,2),(30,40),(5,6),(7,8)], a, "cvector[i] = real")
    a[1:2] = [(300,400),(500,600)]
    test([(1,2),(300,400),(500,600),(7,8)], a, "cvector[i:i] = cvector")
    a[1:2:2] = [(3000,4000), (7000,8000)]
    test([(1,2),(3000,4000), (500,600), (7000,8000)], a, "cvector[i:i:i] = vector")
end

function test_matrix_index()

    ~ simple index

    a = [1,2,3|4,5,6|7,8,9]
    test([4,5,6], a[1], "matrix[i]")
    test([4,5,6|7,8,9], a[1:2], "matrix[i:i]")
    test([1,2,3|7,8,9], a[0:2:2], "matrix[i:i:i]")
    a[1] = [40,50,60]
    test([1,2,3|40,50,60|7,8,9], a, "matrix[i] = vector")
    a[1:2] = [44,55,66|77,88,99]
    test([1,2,3|44,55,66|77,88,99], a, "matrix[i:i] = matrix")
    a[0:2:2] = [10,20,30|70,80,90]
    test([10,20,30|44,55,66|70,80,90], a, "matrix[i:i:i] = matrix")

    ~ row index

    a = [1,2,3|4,5,6|7,8,9]
    test([4,5,6], a[1,], "matrix[i,]")
    test([4,5,6|7,8,9], a[1:2,], "matrix[i:i,]")
    test([1,2,3|7,8,9], a[0:2:2,], "matrix[i:i:i,]")
    a[1,] = [40,50,60]
    test([1,2,3|40,50,60|7,8,9], a, "matrix[i,] = vector")
    a[1:2,] = [44,55,66|77,88,99]
    test([1,2,3|44,55,66|77,88,99], a, "matrix[i:i,] = matrix")
    a[0:2:2,] = [10,20,30|70,80,90]
    test([10,20,30|44,55,66|70,80,90], a, "matrix[i:i:i,] = matrix")

    ~ column index

    a = [1,2,3|4,5,6|7,8,9]
    test([2, 5, 8], a[,1], "matrix[,i]")
    test([2,3|5,6|8,9], a[,1:2], "matrix[,i:i]")
    test([1,3|4,6|7,9], a[,0:2:2], "matrix[,i:i:i]")
    a[,1] = [20,50,80]
    test([1,20,3|4,50,6|7,80,9], a, "matrix[,i] = vector")
    a[,1:2] = [22,33|55,66|88,99]
    test([1,22,33|4,55,66|7,88,99], a, "matrix[,i:i] = matrix")
    a[,0:2:2] = [10,30|40,60|70,90]
    test([10,22,30|40,55,60|70,88,90], a, "matrix[,i:i:i] = matrix")

    ~ matrix index

    a = [1,2,3|4,5,6|7,8,9]
    test(6, a[1,2], "matrix[i,i]")
    test([5,6], a[1,1:2], "matrix[i,i:i]")
    test([4,6], a[1,0:2:2], "matrix[i,i:i:i]")
    test([5,8], a[1:2,1], "matrix[i:i,i]")
    test([5,6|8,9], a[1:2,1:2], "matrix[i:i,i:i]")
    test([4,6|7,9], a[1:2,0:2:2], "matrix[i:i,i:i:i]")
    test([2,8], a[0:2:2,1], "matrix[i:i:i,i]")
    test([2,3|8,9], a[0:2:2,1:2], "matrix[i:i:i,i:i]")
    test([1,3|7,9], a[0:2:2,0:2:2], "matrix[i:i:i,i:i:i]")

    a[1,2] = 60
    test([1,2,3|4,5,60|7,8,9], a, "matrix[i,i] = real")
    a[1,1:2] = [55,66]
    test([1,2,3|4,55,66|7,8,9], a, "matrix[i,i:i] = vector")
    a[1,0:2:2] = [400,600]
    test([1,2,3|400,55,600|7,8,9], a, "matrix[i,i:i:i] = vector")
    a[1:2,1] = [50,80]
    test([1,2,3|400,50,600|7,80,9], a, "matrix[i:i,i] = vector")
    a[1:2,1:2] = [555,666|888,999]
    test([1,2,3|400,555,666|7,888,999], a, "matrix[i:i,i:i] = matrix")
    a[1:2,0:2:2] = [404,606|707,909]
    test([1,2,3|404,555,606|707,888,909], a, "matrix[i:i,i:i:i] = matrix")
    a[0:2:2,1] = [20,80]
    test([1,20,3|404,555,606|707,80,909], a, "matrix[i:i:i,i] = vector")
    a[0:2:2,1:2] = [22,33|88,99]
    test([1,22,33|404,555,606|707,88,99], a, "matrix[i:i:i,i:i] = matrix")
    a[0:2:2,0:2:2] = [10,30|70,90]
    test([10,22,30|404,555,606|70,88,90], a, "matrix[i:i:i,i:i:i] = matrix")
end

function test_cmatrix_index()

    ~ simple index

    a = [(1,2),(3,4),(5,6)
         (7,8),(9,10),(11,12)
	 (13,14),(15,16),(17,18)]

    test([(7,8),(9,10),(11,12)], a[1], "cmatrix[i]")
    test([(7,8),(9,10),(11,12)|(13,14),(15,16),(17,18)], a[1:2], "cmatrix[i:i]")
    test([(1,2),(3,4),(5,6)|(13,14),(15,16),(17,18)], a[0:2:2], "cmatrix[i:i:i]")
    a[1] = [(70,80),(90,100),(110,120)]
    test([(1,2),(3,4),(5,6)|(70,80),(90,100),(110,120)|(13,14),(15,16),(17,18)], a, "cmatrix[i] = cvector")
    a[1:2] = [(77,88),(99,101),(111,122)|(133,144),(155,166),(177,188)]
    test([(1,2),(3,4),(5,6)|(77,88),(99,101),(111,122)|(133,144),(155,166),(177,188)], a, "cmatrix[i:i] = cmatrix")
    a[0:2:2] = [(11,22),(33,44),(55,66)|(31,41),(51,61),(71,81)]
    test([(11,22),(33,44),(55,66)|(77,88),(99,101),(111,122)|(31,41),(51,61),(71,81)], a, "cmatrix[i:i:i] = cmatrix")

    ~ row index

    a = [(1,2),(3,4),(5,6)
         (7,8),(9,10),(11,12)
	 (13,14),(15,16),(17,18)]

    test([(7,8),(9,10),(11,12)], a[1,], "cmatrix[i,]")
    test([(7,8),(9,10),(11,12)|(13,14),(15,16),(17,18)], a[1:2,], "cmatrix[i:i,]")
    test([(1,2),(3,4),(5,6)|(13,14),(15,16),(17,18)], a[0:2:2,], "cmatrix[i:i:i,]")
    a[1,] = [(70,80),(90,100),(110,120)]
    test([(1,2),(3,4),(5,6)|(70,80),(90,100),(110,120)|(13,14),(15,16),(17,18)], a, "cmatrix[i,] = cvector")
    a[1:2,] = [(77,88),(99,101),(111,122)|(133,144),(155,166),(177,188)]
    test([(1,2),(3,4),(5,6)|(77,88),(99,101),(111,122)|(133,144),(155,166),(177,188)], a, "cmatrix[i:i,] = cmatrix")
    a[0:2:2,] = [(11,22),(33,44),(55,66)|(31,41),(51,61),(71,81)]
    test([(11,22),(33,44),(55,66)|(77,88),(99,101),(111,122)|(31,41),(51,61),(71,81)], a, "cmatrix[i:i:i,] = cmatrix")

    ~ column index

    a = [(1,2),(3,4),(5,6)
         (7,8),(9,10),(11,12)
	 (13,14),(15,16),(17,18)]

    test([(3,4),(9,10),(15,16)], a[,1], "cmatrix[,i]")
    test([(3,4),(5,6)|(9,10),(11,12)|(15,16),(17,18)], a[,1:2], "cmatrix[,i]")
    test([(1,2),(5,6)|(7,8),(11,12)|(13,14),(17,18)], a[,0:2:2], "cmatrix[,i:i:i]")
    a[,1] = [(30,40),(90,100),(150,160)]
    test([(1,2),(30,40),(5,6)|(7,8),(90,100),(11,12)|(13,14),(150,160),(17,18)], a, "cmatrix[,i] = cvector")
    a[,1:2] = [(33,44),(55,66)|(99,101),(111,122)|(155,166),(177,188)]
    test([(1,2),(33,44),(55,66)|(7,8),(99,101),(111,122)|(13,14),(155,166),(177,188)], a, "cmatrix[,i:i] = cmatrix")
    a[,0:2:2] = [(11,22),(50,60)|(77,88),(11,21)|(31,41),(71,81)]
    test([(11,22),(33,44),(50,60)|(77,88),(99,101),(11,21)|(31,41),(155,166),(71,81)], a, "cmatrix[,i:i:i] = cmatrix")

    ~ matrix index

    a = [(1,2),(3,4),(5,6)
         (7,8),(9,10),(11,12)
	 (13,14),(15,16),(17,18)]

    test((11,12), a[1,2], "cmatrix[i,i]")
    test([(9,10),(11,12)], a[1,1:2], "cmatrix[i,i:i]")
    test([(7,8),(11,12)], a[1,0:2:2], "cmatrix[i,i:i:i]")
    test([(9,10),(15,16)], a[1:2,1], "cmatrix[i:i,i]")
    test([(9,10),(11,12)|(15,16),(17,18)], a[1:2,1:2], "cmatrix[i:i,i:i]")
    test([(7,8),(11,12)|(13,14),(17,18)], a[1:2,0:2:2], "cmatrix[i:i,i:i:i]")
    test([(3,4),(15,16)], a[0:2:2,1], "cmatrix[i:i:i,i]")
    test([(3,4),(5,6)|(15,16),(17,18)], a[0:2:2,1:2], "cmatrix[i:i:i,i:i]")
    test([(1,2),(5,6)|(13,14),(17,18)], a[0:2:2,0:2:2], "cmatrix[i:i:i,i:i:i]")

    a[1,2] = (111,112)
    test([(1,2),(3,4),(5,6)|(7,8),(9,10),(111,112)|(13,14),(15,16),(17,18)], a, "cmatrix[i,i] = complex")
    a[1,1:2] = [(209,210),(211,212)]
    test([(1,2),(3,4),(5,6)|(7,8),(209,210),(211,212)|(13,14),(15,16),(17,18)], a, "cmatrix[i,i:i] = cvector")
    a[1,0:2:2] = [(307,308),(311,312)]
    test([(1,2),(3,4),(5,6)|(307,308),(209,210),(311,312)|(13,14),(15,16),(17,18)], a, "cmatrix[i,i:i:i] = cvector")
    a[1:2,1] = [(409,410),(415,416)]
    test([(1,2),(3,4),(5,6)|(307,308),(409,410),(311,312)|(13,14),(415,416),(17,18)], a, "cmatrix[i:i,i] = cvector")
    a[1:2,1:2] = [(509,510),(511,512)|(515,516),(517,518)]
    test([(1,2),(3,4),(5,6)|(307,308),(509,510),(511,512)|(13,14),(515,516),(517,518)], a, "cmatrix[i:i,i:i] = cmatrix")
    a[1:2,0:2:2] = [(607,608),(611,612)|(613,614),(617,618)]
    test([(1,2),(3,4),(5,6)|(607,608),(509,510),(611,612)|(613,614),(515,516),(617,618)], a, "cmatrix[i:i,i:i:i] = cmatrix")
    a[0:2:2,1] = [(703,704),(715,716)]
    test([(1,2),(703,704),(5,6)|(607,608),(509,510),(611,612)|(613,614),(715,716),(617,618)], a, "cmatrix[i:i:i,i] = cvector")
    a[0:2:2,1:2] = [(803,804),(805,806)|(815,816),(817,818)]
    test([(1,2),(803,804),(805,806)|(607,608),(509,510),(611,612)|(613,614),(815,816),(817,818)], a, "cmatrix[i:i:i,i:i] = cmatrix")
    a[0:2:2,0:2:2] = [(901,902),(905,906)|(913,914),(917,918)]
    test([(901,902),(803,804),(905,906)|(607,608),(509,510),(611,612)|(913,914),(815,816),(917,918)], a, "cmatrix[i:i:i,i:i:i] = cmatrix")

end

function test_string_index()
    a = "abcdefghijklmnopqrstuvwxyz"
    test("m", a[12], "string[i]")
    test("bcd", a[1:3], "string[i:i]")
    test("bdf", a[1:3:2], "string[i:i:i]")
    b = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
    b[12] = "mark"
    test("ABCDEFGHIJKLmNOPQRSTUVWXYZ", b, "string[i] = s")
    b[0:13] = a[13:13]
    test("nopqrstuvwxyzNOPQRSTUVWXYZ", b, "string[i:i] = s")
    b = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
    c = string(26)
    c[1:13:2] = a
    c[0:13:2] = b
    test("AaBbCcDdEeFfGgHhIiJjKkLlMm", c, "string[i:i:i]")
end

function test_list_index()
    a = {1, 2, 3, 4}
    test(2, a[1], "list[i]")
    test({2, 3}, a[1:2], "list[i:i]")
    test({2, 4}, a[1:2:2], "list[i:i:i]")
    a[1] = {2, 2.5}
    test({1, {2, 2.5}, 3, 4}, a, "list[i] = value")
    a[1:2] = {20, 30}
    test({1, 20, 30,4}, a, "list[i:i] = list")
    a[1:2:2] = {200, 400}
    test({1, 200, 30, 400}, a, "list[i:i:i] = list")
end

function test_index()
    test_vector_index()
    test_cvector_index()
    test_matrix_index()
    test_cmatrix_index()
    test_string_index()
    test_list_index()
end

function test_standard_functions()
    ~ not testing matvec(), exit(), read(), eof()
    ~ some tests may fail due to machine rounding
    test([0, 0, 0], vector(3), "vector()")
    test("   ", string(3), "string()")
    test({0, 0, 0}, list(3), "list()")
    test(3, size([0, 0, 0]), "size(vector)")
    test(3, size("abc"), "size(string)")
    test(3, size({1, 2, 3}), "size(list)")
    test({1,2,3,4,5,6}, concat({1,2,3},{4,5,6}), "concat(list,list)")
    test(REAL_TYPE, type(0), "type(e)")
    test(65, ascii("A"), "ascii(string)")
    test(1, real((1,2)), "real(complex)")
    test(2, imag((1,2)), "imag(complex)")
    test(2, abs(-2), "abs(real)")
    test(5, abs((3,4)), "abs(complex)")
    test(0, arg((1,0)), "arg(complex)")
    test(25, norm((3,4)), "norm(complex)")
    test((1,-2), conj((1,2)), "conj(complex)")
    test((1,0), polar(1,0), "polar(real,real)")
    test(5, sqrt(25), "sqrt(real)")
    test((1,2), sqrt((-3,4)), "sqrt(complex)")
    test(1, exp(0), "exp(real)")
    test((1,0), exp((0,0)), "exp(complex)")
    test(0, expm1(0), "exp(real)")
    test(0, log(1), "log(real)")
    test((0,0), log((1,0)), "log(complex)")
    test(0, log1p(0), "log1p(real)")
    test(1, log10(10), "log10(real)")
    test(0, sin(0), "sin(real)")
    test((0,0), sin((0,0)), "sin(complex)")
    test(1, cos(0), "cos(real)")
    test((1,0), cos((0,0)), "cos(complex)")
    test(0, tan(0), "tan(real)")
    test((0,0), tan((0,0)), "tan(complex)")
    test(1, csc(PI/2), "csc(real)")
    test((1,0), csc((PI/2,0)), "csc(complex)")
    test(1, sec(0), "sec(real)")
    test((1,0), sec((0,0)), "sec(complex")
    test(1+EPS, cot(PI/4), "cot(real)")
    test((1+EPS,0), cot((PI/4,0)), "cot(complex)")
    test(0, asin(0), "asin(real)")
    test((0,0), asin((0,0)), "asin(complex)")
    test(0, acos(1), "acos(real)")
    test((0,0), acos((1,0)), "acos(complex)")
    test(0, atan(0), "atan(real)")
    test((0,0), atan((0,0)), "atan(complex)")
    test(PI/2, acsc(1), "acsc(real)")
    test((PI/2,0), acsc((1,0)), "acsc(comple)")
    test(0, asec(1), "sec(real)")
    test((0,0), asec((1,0)), "sec(complex)")
    test(PI/4, acot(1), "acot(real)")
    test((PI/4,0), acot((1,0)), "acot(complex)")
    test(0, sinh(0), "sinh(real)")
    test((0,0), sinh((0,0)), "sinh(complex)")
    test(1, cosh(0), "cosh(real)")
    test((1,0), cosh((0,0)), "cosh(complex)")
    test(0, tanh(0), "tanh(real)")
    test((0,0), tanh((0,0)), "tanh(complex)")
    test(0.850918128239321521, csch(1), "csch(real)")
    test((0.850918128239321521,0), csch((1,0)), "csch(complex)")
    test(1, sech(0), "sech(real)")
    test((1,0), sech((0,0)), "sech(complex)")
    test(0.642092615934330624, coth(1), "coth(real)")
    test((0.642092615934330735,0), coth((1,0)), "coth(real)")
    test(0, asinh(0), "asinh(real)")
    test((0,0), asinh((0,0)), "asinh(complex)")
    test(0, acosh(1), "acosh(real)")
    test((0,0), acosh((1,0)), "acosh(complex)")
    test(0, atanh(0), "atanh(real)")
    test((0,0), atanh((0,0)), "atanh(complex)")
    test(0.88137358701954303, acsch(1), "acsch(real)")
    test((0.881373587019542899,0), acsch((1,0)), "acsch(complex)")
    test(0, asech(1), "asech(real)")
    test((0,0), asech((1,0)), "asech(complex)")
    test(0.549306144334054734, acoth(2), "acoth(real)")
    test((0.549306144334054845,0), acoth((2,0)), "acoth(complex)")
    test(120, gamma(6), "gamma(real)")
end

function mpl_test()
    test_arithmetic()
    test_logic()
    test_relational_operators()
    test_index()
    test_standard_functions()
end

mpl_test()

"Done testing"

max_stars_repo_name: yugyfoog/MPL

max_issues_repo_name: yugyfoog/MPL

max_forks_repo_name: yugyfoog/MPL


text:
######################################################################
# An ordering of a set A is encoded as a list containing each element
# of A precisely once.

`is_element/ord` := (A::set) -> proc(R)
 type(R,list) and nops(R) = nops(A) and {op(R)} = A;
end;

`is_equal/ord` := (A::set) -> (R,S) -> evalb(R = S);

`is_leq/ord` := NULL;

`random_element/ord` := (A::set) -> () -> combinat[randperm]([op(A)]);

`list_elements/ord` := (A::set) -> combinat[permute](sort([op(A)]));

`count_elements/ord` := (A::set) -> nops(A)!;

`rank_table/ord` := (A::set) -> proc(R)
 local i;
 table([seq(R[i] = i, i = 1 .. nops(R))]):
end:


text:
######################################################################
# pairs(A) is the set of all pairs (a,b) in A^2 with a <> b

`is_element/pairs` := (A::set) -> proc(ab)
 type(ab,list) and nops(ab) = 2 and
 member(ab[1],A) and member(ab[2],A) and ab[1] <> ab[2];
end;

`is_equal/pairs` := (A::set) -> (ab,cd) -> evalb(ab = cd):

`is_leq/pairs` := NULL;

`random_element/pairs` := (A::set) -> proc()
 local n,r,i,j;
 n := nops(A);

 if n < 2 then return FAIL; fi;

 i := rand(1..n)();
 j := rand(1..n-1)();
 if j >= i then j := j+1; fi;

 return [A[i],A[j]];
end;

`list_elements/pairs` := proc(A::set)
 local n,i,j;

 n := nops(A);

 [seq(op([seq([A[i],A[j]],j=1..i-1),seq([A[i],A[j]],j=i+1..n)]),i=1..n)];
end:

`count_elements/pairs` := (A::set) -> nops(A) * (nops(A)-1);


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

g96_c1 := 1/137:

f_g96  := x->
  1 + g96_c1/X_FACTOR_C * x^(3/2):

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


text:
$ifndef _CLASSIFYER_
$define _CLASSIFYER_

$include "Basic.mpl"
$include "Condition.mpl"
$include "Fetch.mpl"
$include "InvOrder.mpl"
$include "InvSimplify.mpl"
$include "Logout.mpl"
$include "Utils.mpl"
$include "InvSol.mpl"
$include "Closure.mpl"
$include "ClassifyHolder.mpl"
$include "GenSol.mpl"
$include "TeqSol.mpl"
$include "TsolsList.mpl"

# 绑定函数
reset:=ClassifyHolder:-reset;
addSol:=ClassifyHolder:-addSol;
getSols:=ClassifyHolder:-getSols;
getIeqCode:=ClassifyHolder:-getIeqCode;
addUnsolvedSol:=ClassifyHolder:-addUnsolvedSol;
getUnsolvedSols:=ClassifyHolder:-getUnsolvedSols;

# 重新求解的入口
classify:=proc(As,A,eqs)
    local sol,n;
    sol:=Object(InvSol);
    sol:-As:=As;
    sol:-A:=A;
    n:=LinearAlgebra:-RowDimension(A);
    sol:-nvars:=n;
    sol:-vars:={seq(a[i],i=1..n)};
    sol:-oeq:=eqs;
    sol:-state:=0;
    reset();
    resolve(sol);
end proc:

# 可重用的求解入口
resolve:=proc(s::InvSol)
    flogf[0](convert(procname,string));
    if      (s:-state=0) then
        return solveOeq(s);     # 针对偏微分方程组进行求解
    elif    (s:-state=1) then
        return solveIeq(s);     # 针对不变量方程进行求解
    elif    (s:-state=2) then
        return checkNewInvariants(s); # 检查是否产生了新的不变量
    elif    (s:-state=3) then
        return solveTeq(s);     # 求解变换方程
    elif    (s:-state=4) then
        return selectRep(s);    # 选择最优代表元
    elif    (s:-state=5) then
        addSol(s);              # 求解成功
        return;
    else
        error "unkown state";
    end if;
end proc:

# 针对偏微分方程组进行求解
solveOeq:=proc(s::InvSol)
    local deltas;
    flogf[0](convert(procname,string));
    deltas:=getNewInvariants(s);
    if deltas=[] or convert(deltas,set) subset convert(s:-Deltas,set) then
        flogf[1]("没有新的不变量");
        solveByClosure(s);
    else
        flogf[1]("解得新的不变量");
        flog[1]~(deltas);
        genIeq(s,deltas);
    end if;
end proc:

# 求解新的不变量
# 两处判断用的都是这个函数，因此要考虑通用性，代入的条件到底是什么，
# 该如何定义，该如何和别的地方进行协调？
getNewInvariants:=proc(s::InvSol)
    local oeq,deltas;
    flogf[0](convert(procname,string));
    flogf[0]("-------------------------------------------------");
    flogf[0]("求解新的不变量");
    flogf[0]("附加约束");
    flog[0](s:-addcons);

    oeq:=getRealOeq(s);
    flogf[0]("偏微分方程组");
    flog[0](oeq);

    # 求解新不变量（已化简）
    deltas:=getInvariants(oeq);

    # 检查是否有解
    if not type(indets(deltas,name),set(specindex(a))) then
        deltas:=[];
    end if;

    return deltas;
end proc:

# 按封闭进行求解
solveByClosure:=proc(s::InvSol)
    local c,s1,s2;
    flogf[0](convert(procname,string));
    c:=getMinClosure(getClosure(s:-A,getZeroCons(s)));
    c:={seq(a[x],x in c)};
    # 封闭不全为零
    flogf[1]("封闭不全为零");
    flog[1](c);
    s1:=Object(s);
    # 添加展示性约束
    s1:-discons:=s1:-discons union {add(x^2,x in c)<>0};
    s1:-rep:=v[op([1,1],c)];
    s1:-state:=5;
    resolve(s1);
    # 封闭全为零
    s2:=Object(s);
    # 这里选择了在addcons中加入信息，而不是加入到解中
    solveClosureAllZero(s2,c);
end proc:

# 处理封闭全为零的情况
solveClosureAllZero:=proc(s::InvSol,c::set(specindex(a)))
    local s1;
    flogf[0](convert(procname,string));
    flogf[1]("-------------------------------------------------");
    flogf[1]("封闭全为零");
    flog[1](c);
    addZeroCons(s,{seq(x=0,x in c)});
    updateVars(s);
    if numelems(s:-vars)>1 then
        solveOeq(s);
    elif numelems(s:-vars)=1 then
        s1:=Object(s);
        s1:-rep:=v[op([1,1],s:-vars)];
        s1:-state:=5;
        resolve(s1);
    else
        return;
    end if;
end proc:

# 生成不变量方程组
genIeq:=proc(s::InvSol,deltas)
    local spos,pos,n;
    flogf[0](convert(procname,string));
    if andmap(type,rhs~(s:-ieq),0) then
        spos:=numelems(s:-Deltas)+1;
        s:-Deltas:=[s:-Deltas[],deltas[]];
        s:-orders:=findInvariantsOrder~(s:-Deltas);
        n:=numelems(s:-Deltas);
        for pos from spos to n do
            buildIeq(s,pos);
        end do;
    else
        appendIeq(s,deltas);
    end if;
end proc:

# 按照拓展的方式生成不变量方程
# 对于 Delta[1]=c[1],...,Delta[n]=c[n]的方程，若生成了新的不变量Delta[n+1]
# 新的方程为Delta[1]=c[1],...,Delta[n]=c[n],Delta[n+1]=c[n+1]
appendIeq:=proc(_s::InvSol,deltas)
    local s,cid,getCname;
    flogf[0](convert(procname,string));
    cid:=findCid(_s);
    getCname:=proc()
        cid:=cid+1;
        return c[cid];
    end proc:
    s:=Object(_s);
    s:-Deltas:=[s:-Deltas[],deltas[]];
    s:-orders:=findInvariantsOrder~(s:-Deltas);
    s:-ieq:=[s:-ieq[],seq(deltas[i]=getCname(),i=1..numelems(deltas))];
    s:-state:=1;
    resolve(s);
end proc:

# 按位置进行讨论，建立不变量方程
# 讨论 Delta[k] 的取值
# 此时 Delta[1..(k-1)]=0，Delta[(k+1)..n]=c[j]
buildIeq:=proc(_s::InvSol,pos::posint)
    local s,n,xpos,cid,getCname,x,rs;
    flogf[0](convert(procname,string));
    n:=numelems(_s:-Deltas);
    if type(_s:-orders[pos],even) then
        xpos:=[1,-1,0];
    else
        xpos:=[1,0];
    end if;
    cid:=0;# 因为前面全是0
    getCname:=proc()
        cid:=cid+1;
        return c[cid];
    end proc:
    rs:=Array(1..n,x->`if`(x>pos,getCname(),0));
    for x in xpos do
        # 只求解全零方程，其它情况将包含在下一个pos的情况中
        if x=0 and pos<> n then
            next;
        end if;
        rs[pos]:=x;
        s:=Object(_s);
        s:-ieqCode:=getIeqCode();
        s:-ieq:=[seq(s:-Deltas[i]=rs[i],i=1..n)];
        s:-state:=1;
        resolve(s);
    end do;
end proc:

# 获取当前c的最后一个下标
# 本做法保证一个不变量方程中的任意常数c，从c[1]开始依次编号
# 而不与其它不变量方程中的任意常数c进行比较
findCid:=proc(s::InvSol)
    local rs:=select(type,rhs~(s:-ieq),specindex(c));
    flogf[0](convert(procname,string));
    if rs=[] then
        return 0;
    else
        return op(1,rs[-1]);
    end if;
end proc:

# 是否等价于封闭全为零的解
isClosureSol:=proc(s::list)
    return andmap(x->evalb(rhs(x)=0 or lhs(x)=rhs(x)),s);
end proc:

# 求解不变量方程
solveIeq:=proc(s::InvSol)
    local isols,icons,acons,tc,cc,discons,subcons,ind,i,isol;
    flogf[0](convert(procname,string));
    flogf[1]("-------------------------------------------------");
    flogf[1]("求解不变量方程");
    displayIeq(s);

    # 针对全部变量求解
    isols,icons:=ieqsolve(s:-ieq,{seq(a[i],i=1..s:-nvars)});

    flog[1](isols);
    flogf[1]("约束条件");
    flog[1](icons);

    # 等价于封闭全为零进行求解
    if andmap(isClosureSol,isols) then
        flogf[1]("等价于封闭为零");
        for isol in isols do
            solveClosureAllZero(Object(s),
            indets(select(x->rhs(x)=0,isol),name));
        end do;
        return;
    end if;

    # 进入封闭的求解条件
    # 若只存在单变量大于零或小于零的的约束，则等价成非零约束进行求解。
    # 只存在常数约束和非零约束
    # 只存在一个具有单个单变量非零约束的一般解
    # 若存在多个具有单个单变量非零约束的一般解，
    #     若两个解不能比较，则选择其中一个进行求解。
    #     若两个解等价，则分别进行封闭求解。
    # 若不能进入封闭进行求解，则一个解对应一个特解分别进行求解

    # 首先满足只存在单变量约束
    acons:=`union`(icons[]);
    if not andmap(x->evalb(numelems(indets(x,name))=1),acons) then
        return branchSolve(s,isols,icons);
    end if;

    # 尝试将大于零或小于零的约束等价为非零约束
    # 单变量约束只有3种可能：= <> <
    acons:=select(type,acons,`<`);
    if acons<>{} then
        tc:=table();
        for cc in acons do
            if lhs(cc)=0 or rhs(cc)=0 then
                tappend(tc,indets(cc,name)[],cc);
            else
                # 存在不是大于零或者小于零的约束，则不进入封闭进行计算
                return branchSolve(s,isols,icons);
            end if;
        end do;
        # 对于同一个变量存在大于零和小于零的约束
        # 一般这个约束会体现在值域的约束中，但是现在还没有计算值域
        # 而且从其次性进行考虑，出现这种值域约束的同时，肯定会出现根号约束，
        # 大部分情况下根号约束是不能处理的，不进行封闭计算。
        # 可以给出一个例子，a[1]^2-a[2]^2=0，可以解出 a[1]=±a[2]
        # 反正太复杂了我不走封闭了
        if ormap(x->evalb(numelems(x)>1),[entries(tc,nolist)]) then
            return branchSolve(s,isols,icons);
        end if;
        discons:={entries(tc,nolist)};
        discons:=map(x->x[],discons);
        subcons:=seq((x)=(indets(x)[]<>0),x in discons);
        # 处理掉<约束之后，只剩下单变量等式约束和非零约束，必然能够进行封闭运算
        s:-discons:=s:-discons union discons;
        icons:=subs(subcons,icons);
    end if;

    # 然后至少要有一个非零约束
    if not ormap(isNonZeroCon,`union`(icons[])) then
        return branchSolve(s,isols,icons);
    end if;

    # 尝试进入封闭
    # TODO 还没加入一般解只有一个非零约束的判定
    ind:=findGenSolInd(icons,s:-nvars);
    if type(ind,posint) then
        # 只有一个极大元
        closureRefine(s,isols,icons,ind);
    else
        # 有多个极大元，在不能比较的部分中选择包含元素个数最少的那组
        ind:=MinSelect(ind,numelems)[];
        for i in ind do
            closureRefine(s,isols,icons,i);
        end do;
    end if;

end proc:

# 按照分支方法进行求解
branchSolve:=proc(s::InvSol,isols,icons)
    local i,n,_s;
    flogf[0](convert(procname,string));
    flogf[1]("分支求解");
    n:=numelems(isols);
    for i from 1 to n do
        _s:=Object(s);
        _s:-isols:=[isols[i]];
        _s:-icons:=[icons[i]];
        _s:-state:=2;
        _s:-isolInd:=1;
        _s:-useBranch:=true;
        resolve(_s);
    end do;
end proc:

# 解的精简，按照封闭求解
# 方程的解至多只含等式约束和非零的约束
closureRefine:=proc(_s::InvSol,isols,icons,ind::posint)
    local CL,zcons,s,s0,s1;
    flogf[0](convert(procname,string));
    flogf[1]("一般解为");
    flog[1](isols[ind]);
    s:=Object(_s);
    # 取封闭
    CL:=getClosure(s:-A,getZeroCons(s));
    zcons:=select(isNonZeroCon,icons[ind]);# 提取约束
    zcons:=indets(zcons,name);# 提取变量
    zcons:=map(x->op(1,x),zcons);# 提取下标
    zcons:=[zcons[]];
    CL:=`union`(CL[zcons][]);
    CL:={seq(a[x],x in CL)};# 系数表示
    flogf[1]("封闭为");
    flog[1](CL);
    # 封闭不全为零
    s1:=Object(s);
    s1:-state:=2;
    s1:-isols:=isols;
    s1:-icons:=icons;
    s1:-isolInd:=ind;
    s1:-discons:=s1:-discons union {add(x^2,x in CL)<>0};
    resolve(s1);
    # 封闭全为零
    if checkIeq(s,CL) then
        flogf[1]("封闭全为零有解");
        s0:=Object(s);
        solveClosureAllZero(s0,CL);
    else
        flogf[1]("封闭全为零无解");
    end if;
end proc:

# 检查封闭全为零是否满足方程
checkIeq:=proc(s::InvSol,c::set(specindex(a)))
    local eq,e;
    try
        eq:=subs(seq(x=0,x in c),s:-ieq);
        eq:=remove(x->lhs(x)=rhs(x),eq);
        for e in eq do
            if not isVar(lhs(e)) then
                if isVar(rhs(e)) then
                    return false;
                else
                    if lhs(x)<>rhs(x) then
                        return false;
                    end if;
                end if;
            end if;
        end do;
        return true;
    catch :
        return false;
    end try;
end proc:

isVar:=proc(x)
    return indets(x,name)<>{};
end proc:

# 取特解
checkNewInvariants:=proc(s::InvSol)
    local deltas;
    flogf[0](convert(procname,string));
    deltas:=getNewInvariants(s);
    if deltas=[] or convert(deltas,set) subset convert(s:-Deltas,set) then
        flogf[1]("没有新的不变量");
        solveRep(s);
    else
        flogf[1]("解得新的不变量");
        flog[1]~(deltas);
        genIeq(s,deltas);
    end if;
end proc:

# 取特解
solveRep:=proc(s::InvSol)
    local rsols,_s,rsol;
    flogf[0](convert(procname,string));
    rsols:=fetchSpecSol~(s:-isols,s:-icons,nonzero);# 针对所有解取特解
    flogf[1]("一般解为");
    flog[1](s:-isols[s:-isolInd]);
    flogf[1]("取特解");
    rsols:=`union`(rsols[]);
    rsols:=convert(rsols,list);
    flog[1](rsols);
    if not s:-useBranch then
        s:-rsols:=rsols;
        s:-state:=3;
        resolve(s);
    else
        for rsol in rsols do
            _s:=Object(s);
            _s:-rsols:=[rsol];
            _s:-state:=3;
            resolve(_s);
        end do;
    end if;
    return;
end proc:

# 建立不变量方程并求解
solveTeq:=proc(s::InvSol)
    local ESC;
    ESC:=map[2](specTeqSolve,s,s:-rsols);
    if not ormap(TeqSol:-hasSol,ESC) then
        WARNING("-----------------------------变换方程均无解\n");
        addUnsolvedSol(s);
        return;
    end if;
    s:-tsols:=ESC;
    s:-state:=4;
    return resolve(s);
end proc:

specTeqSolve:=proc(sol::InvSol,spec::list)
    local ax,_ax,res;
    ax:=Matrix([seq(a[i],i=1..sol:-nvars)]);
    _ax:=Matrix(spec);
    res:=TeqSol(spec,[solveSpecTeq(ax,_ax,sol),
                      solveSpecTeq(_ax,ax,sol)]);
    return res;
end proc:

solveSpecTeq:=proc(va,vb,s::InvSol)
    local teq,tsol,tcon,var,n,eqs,eq,_eq,_con,_sol,needResolve;
    n:=numelems(va);
    teq:=subs(s:-isols[s:-isolInd][],convert(va-vb.s:-A,list));
    var:={seq(epsilon[i],i=1..n)};
    flogf[0]("求解变换方程");
    flog[0](teq);
    needResolve:=false;
    try # 设置最长求解时间
        tsol:=timelimit(3,teqsolve(teq,var));
    catch:
        flogf[1]("直接求解超时");
        flog[1](teq);
        needResolve:=true;
    end try;
    if needResolve or tsol=[] then
        # 求解失败，尝试二次求解法方法
        # 首次求解
        eqs:=teqsolve(teq);
        # 二次求解
        tsol:=[];
        tcon:=[];
        for eq in eqs do
            if not checkTeqSol(teq,eq) then
                _sol:=[[]];
                _con:=[{}];
            else
                _eq,_con:=selectremove(has,eq,epsilon);
                _con:=remove(x->type(x,`=`) and (lhs(x)=rhs(x)),_con);
                _con:=convert(_con,set);
                _sol:=teqsolve(_eq,var,_explicit);
                _con:=map(x->getTsolCons(x,s) union _con,_sol);
            end if;
            tsol:=[tsol[],_sol[]];
            tcon:=[tcon[],_con[]];
        end do;
    else
        # 求解成功，直接计算约束
        tcon:=map(x->getTsolCons(x,s),tsol);
    end if;
    # 清理矛盾解
    tsol:=zip((s,c)->if (undefined in rhs~(c)) then NULL else s end if,
             tsol,tcon);
    tcon:=remove(c->(undefined in rhs~(c)),tcon);
    return [teq,tsol,tcon];
end proc:

# 检查解是否满足条件
checkTeqSol:=proc(teq,sol)
    return andmap(type,RealDomain:-simplify(subs(sol[],teq)),0);
end proc:

# 自定义变换方程求解函数
teqsolve:=proc({_explicit::boolean:=false})
    local res;
    if _explicit then
        res:=[RealDomain:-solve(_rest,explicit)];
    else
        res:=convert~([RealDomain:-solve(_rest)],radical);
    end if;
    return convert~(res,list);
end proc:

# 获取变换方程的解的约束条件
# 删除只含 epsilon 的约束
# 删除 s:-discons 主要是为了删除只能 a[k]>0 的约束，这依赖于 solveIeq 的相关操作。
# 删除通解中存在的约束
getTsolCons:=proc(tsol,s::InvSol)
    return (select(has,findSolutionDomain(tsol),{a,c}) minus getIsolCons(s));
end proc:

# 检查是否是非零约束
# 只考虑单变量非零约束
# 这里只考虑 x<>0 的形式，不考虑 0<>x 的情况
# 因为在求解约束条件时，能够保证把 0<>x 变成 x<>0
isNonZeroCon:=proc(x)
    return op(0,x)=`<>` and rhs(x)=0 and numelems(indets(x,name))=1;
end proc:

# 不变量方程的求解函数
ieqsolve:=proc(eq::list,vars::set)
    local isols,icons,zcons,ind;
    isols:=convert~([RealDomain:-solve(eq,vars,explicit)],list);
    icons:=findSolutionDomain~(isols);
    # 对于单变量约束的情况，尝试进行补全
    zcons:=select(isNonZeroCon,`union`(icons[]));
    if numelems(zcons)=1 then
        isols:={isols[],convert~([RealDomain:-solve([eq[],indets(zcons,name)[]],vars,explicit)],list)[]};
        isols:=[isols[]];
        isols:=SolveTools[SortByComplexity](isols);
        icons:=findSolutionDomain~(isols);
    end if;
    # 删除对c有约束的解
    ind:=map(k->`if`(chkCCons(icons[k]),NULL,k),[seq(1..numelems(isols))]);
    isols:=isols[ind];
    icons:=icons[ind];
    return isols,icons;
end proc:

# 检查是存在仅和c有关的约束
chkCCons:=proc(con)
    return ormap(x->type(indets(x,name),set(specindex(c))),con);
end proc:

# 选择最优代表元
selectRep:=proc(s::InvSol)
    local tlist;
    tlist:=TsolsList(s:-tsols);
    printTsolsList(tlist);
    s:-tsolsList:=tlist;
    s:-rep:=tlist:-reps[1];
    s:-tInd:=tlist:-torders[1];
    s:-tsol:=tlist:-tsols[1];
    s:-tcon:=tlist:-tcons[1];
    s:-state:=5;
    resolve(s);
end proc:

$endif


text:
(* 2020 Susi Lehtola and Miguel A. L. Marques

   This script is used to check the convergence of the series
   expansions of the various attenuation functions, and to
   establish cutoff values for their series expansions.
*)

$include "attenuation.mpl"

(* For exact results we use 1000 digit *)
exact_digits := 1000:
(* Double precision has 15 digits *)
double_digits := 15:
(* An acceptable relative error at double precision is 1e-13 *)
error_thr := 1e-13:

check_asymptotics := proc(f, a);
    (* Grid spacing *)
    da := 0.01:
    (* Value for a to start with *)
    amax := 5:

    (* First, we compare we find a point where to make the cutoff *)
    for acut from amax by -da to da do
        (* Exact result *)
        Digits := exact_digits:
        exact := evalf(f(acut)):
        (* Double precision *)
        Digits := double_digits:
        local doubleprec := evalf(f(acut)):
        (* Error *)
        local err := doubleprec/exact-1:
        #printf("Cutoff %e exact % e double % e error % e\n", acut, exact, doubleprec, err):
        if (abs(err) < error_thr) then
            break:
        end if:
    end do:

    (* Now we find the series expansion that has the same level of agreement *)
    for expord from 4 to 1000 by 2 do
        f_series := a -> eval(throw_out_large_n(convert(series(f(b), b=infinity, expord+padding_order), polynom), expord), b=a):
        local ser := evalf(f_series(acut)):
        local err := ser/exact-1:
        #printf("Expansion order %3d original % e asymptotic % e error % "
        #       "e\n", expord, exact, ser, err);
        if (abs(err) < error_thr) then
            (* Check if the expansion is accurate everywhere *)
            accurate := true:
            for aval from acut by da to amax do
                (* Exact result *)
                Digits := exact_digits:
                lexact := evalf(f(aval)):
                (* Double precision *)
                Digits := double_digits:
                lser := evalf(f_series(aval)):
                (* Error *)
                local lerr := lser/lexact-1:
                #printf("a= %e exact % e double % e error % e\n", aval, lexact, lser, lerr):
                if (abs(lerr) > error_thr) then
                    accurate := false:
                    break:
                end if:
            end do:
            if accurate then
                break:
            end if:
        end if:
    end do:

    printf("Cutoff should be at a= %.2f and use a series expansion of "
           "order %d\n", acut, expord):
end proc:

printf("attenuation_erf\n"):
check_asymptotics(attenuation_erf0, a):

printf("\nattenuation_erf_f2\n"):
check_asymptotics(attenuation_erf_f20, a):

printf("\nattenuation_erf_f3\n"):
check_asymptotics(attenuation_erf_f30, a):

printf("\nattenuation_gau\n"):
check_asymptotics(attenuation_gau0, a):

printf("\nattenuation_yukawa\n"):
check_asymptotics(attenuation_yukawa0, a):


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

# Autor
# Moritz Schappler, schappler@imes.uni-hannover.de, 2018-02
# Institut fuer mechatronische Systeme, Leibniz Universitaet Hannover
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(LinearAlgebra):
with(StringTools): # Für Zeitausgabe
;
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
codegen_act := true:
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
kin_constraints_exist := true:
# Zwangsbedingungen auslesen und Nachverarbeiten
# Lese Ergebnisse der Kinematik von KAS5m7 aus.
# Entferne dort die Ersetzungsausdrücke
read  sprintf("../codeexport/KAS5m7TE/tmp/kinematic_constraints_maple_inert.m", robot_name):
kintmp_qs := kintmp_qs:
kintmp_qt := kintmp_qt:
kintmp_subsexp := kintmp_subsexp:
# Ausdrücke entfernen
kintmp_subsexp := kintmp_subsexp*0:
#
# Speichern der Ergebnisse
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert", robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
# Liste der Parameter speichern
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_qs )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
kintmp_s:


text:
input := StringTools:-Trim(FileTools:-Text:-ReadFile("AoC-2021-15-input.txt"));
grid := Matrix(map(s->parse~(StringTools:-Explode(s)),StringTools:-Split(input,"\n"))):
(gridw,gridl) := upperbound(grid);

inc := a -> (a mod 9) +~ 1;
gb:= Matrix([seq([seq((inc@@(i+j-2))(grid),j=1..5)],i=1..5)] );
(gbw,gbl) := upperbound(gb);

nbs := proc(i,j,gridsize)
local out := NULL;
    if i < gridsize then out := out, [i+1, j]; end if;
    if j < gridsize then out := out, [i, j+1]; end if;
    if i > 1 then out := out, [i-1,j]; end if;
    if j > 1 then out := out, [i, j-1]; end if;
    return [out];
end proc:

with(GraphTheory):
G := MakeWeighted(MakeDirected(SpecialGraphs:-GridGraph(gridl,gridw)));

for e in Edges(G) do
    SetEdgeWeight(G, e, grid[parse(e[2])]);
end do:

DijkstrasAlgorithm(G, "1,1", cat("",gridl,",",gridw))[2];
BG := MakeWeighted(MakeDirected(SpecialGraphs:-GridGraph(gbl,gbw)));

W := WeightMatrix(BG, copy=false);

for i from 1 to gbl do for j from 1 to gbw do
       a := j + (i-1)*(gbw);
       for n in nbs(i,j,gbw) do
           b := n[2] + (n[1]-1)*gbw; W[b,a] := gb[i,j];
       end do;
end do; end do:

answer2 := DijkstrasAlgorithm(BG, "1,1", cat("",gbl,",",gbw))[2];


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

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

(* type: mgga_exc *)

(* prefix:
  mgga_x_mbr_params *params;

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

(* replace: "br89_x\(" -> "xc_mgga_x_br89_get_x(" *)

$include "mgga_x_br89.mpl"

params_a_at := 0:

(*Equation 15. The three equations below are from mgga_x_tm.mpl; note
that the numerical factors in the denominator in eqn 15 just
correspond to going from 'x' to 's'.*)
tm_p  := x -> (X2S*x)^2:
tm_y  := x -> (2*params_a_lambda - 1)^2 * tm_p(x):
tm_f0 := x -> (1 + 10*(70/27)*tm_y(x) + params_a_beta*tm_y(x)^2)^(1/10):

(* definition below equation 16 *)
mbr_D := (ts, xs) -> 2*ts - 1/4 * (2*params_a_lambda - 1)^2 * xs^2:

(* k_\sigma is not defined in the paper, but Subrata Jana
said on GitLab that k_\sigma = (6\pi^2\rho_\sigma)^{1/3} *)
k_sigma := (6*Pi^2)^(1/3):

(* Equation 18. Note that there's a typo in the paper: since Becke's
tau is two times the kinetic energy density, there should be a factor
of two in front of tau uniform as well. *)
br89_Q := (x, u, t) ->
  1/6*(
  + 6*(params_a_lambda^2 - params_a_lambda + 1/2)*(2*t - 2*K_FACTOR_C - 1/36*x^2)
  + 6/5*k_sigma^2*(tm_f0(x)^2 - 1)
  - 2*params_a_gamma*mbr_D(t, x)
  ):


text:
# Step 2 of 3: computer algebra

improve := proc(lo :: LO(name, anything), {_ctx :: t_kb := empty}, opts := [], $)
local r, `&context`;
  userinfo(5, improve, "input: ", print(lo &context _ctx));
  _Env_HakaruSolve := true;
  r:= LO(op(1,lo), reduce(op(2,lo), op(1,lo), _ctx, opts));
  userinfo(5, improve, "output: ", print(r));
  r
end proc;

# Walk through integrals and simplify, recursing through grammar
# h - name of the linear operator above us
# kb - domain information
reduce := proc(ee, h :: name, kb :: t_kb, opts := [], $)
  local e, elim, subintegral, w, ww, x, c, kb1, with_kb1, dom_specw, dom_specb
       , body, dom_spec, ed, mkDom, vars, rr
       , do_domain := evalb( not ( "no_domain" in {op(opts)} ) ) ;
  e := ee;

  if do_domain then
    rr := reduce_Integrals(e, h, kb, opts);
    if rr <> FAIL then return rr end if;
  end if;
  if e :: 'applyintegrand(anything, anything)' then
    map(simplify_assuming, e, kb)
  elif e :: `+` then
    map(reduce, e, h, kb, opts)
  elif e :: `*` then
    (subintegral, w) := selectremove(depends, e, h);
    if subintegral :: `*` then error "Nonlinear integral %1", e end if;
    subintegral := convert(reduce(subintegral, h, kb, opts), 'list', `*`);
    (subintegral, ww) := selectremove(depends, subintegral, h);
    simplify_factor_assuming(`*`(w, op(ww)), kb)
      * `*`(op(subintegral));
  elif e :: Or(Partition,t_pw) then
    if e :: t_pw then e := PWToPartition(e); end if;
    e := Partition:-Simpl(e);
    e := kb_Partition(e, kb, simplify_assuming,
                      ((rhs, kb) -> %reduce(rhs, h, kb, opts)));
    e := eval(e, %reduce=reduce);
    # big hammer: simplify knows about bound variables, amongst many
    # other things
    Testzero := x -> evalb(simplify(x) = 0);
    e := Partition:-Simpl(e);
    if ee::t_pw and e :: Partition then
      e := Partition:-PartitionToPW(e);
    end if;
    e;
  elif e :: t_case then
    subsop(2=map(proc(b :: Branch(anything, anything))
                   eval(subsop(2='reduce'(op(2,b),x,c,opts),b),
                        {x=h, c=kb})
                 end proc,
                 op(2,e)),
           e);
  elif e :: 'Context(anything, anything)' then
    kb1 := assert(op(1,e), kb);
    # A contradictory `Context' implies anything, so produce 'anything'
    # In particular, 42 :: t_Hakaru = false, so a term under a false
    # assumption should never be inspected in any way.
    if kb1 :: t_not_a_kb then
        return 42
    end if;
    applyop(reduce, 2, e, h, kb1, opts);
  elif e :: 'toLO:-integrate(anything, Integrand(name, anything), list)' then
    x := gensym(op([2,1],e));
    # If we had HType information for op(1,e),
    # then we could use it to tell kb about x.
    subsop(2=Integrand(x, reduce(subs(op([2,1],e)=x, op([2,2],e)), h, kb, opts)), e)
  else
    simplify_assuming(e, kb)
  end if;
end proc;

# "Integrals" refers to any types of "integrals" understood by domain (Int,
# Sum currently)
reduce_Integrals := module()
  export ModuleApply;
  local
  # The callbacks passed by reduce_Integrals to Domain:-Reduce
    reduce_Integrals_body, reduce_Integrals_into
  # tries to evaluate a RootOf
  , try_eval_Root
  # tries to evaluate Int/Sum/Ints/Sums
  , elim_intsum;

  reduce_Integrals_body := proc(h,opts,x,kb1) reduce(x,h,kb1,opts) end proc;
  reduce_Integrals_into := proc(h,opts,kind,e,vn,vt,kb,$)
    local rr;
    rr := elim_intsum(Domain:-Apply:-do_mk(args[3..-1]), h, kb,opts);
    rr := subsindets(rr, specfunc(RootOf), x->try_eval_Root(x,a->a));
    return rr;
  end proc;

  ModuleApply := proc(expr, h, kb, opts, $)
    local rr;
    rr := Domain:-Reduce(expr, kb
      ,curry(reduce_Integrals_into,h,opts)
      ,curry(reduce_Integrals_body,h,opts)
      ,(_->:-DOM_FAIL));
    rr := kb_assuming_mb(Partition:-Simpl)(rr, kb, x->x);
    if has(rr, :-DOM_FAIL) then
      return FAIL;
    elif has(rr, FAIL) then
      error "Something strange happened in reduce_Integral(%a, %a, %a, %a)\n%a"
           , expr, kb, kb, opts, rr;
    end if;
    rr;
  end proc;

  try_eval_Root := proc(e0::specfunc(`RootOf`),on_fail := (_->FAIL), $)
    local ix,e := e0;
    try
      if nops(e)=2 or nops(e)=3
      and op(-1,e) :: `=`(identical(index),{specindex(real),nonnegint}) then
        ix := op([2,-1],e);
        if ix :: specindex(real) then ix := op(ix); end if;
        e := op(0,e)(op(1,e));
      else
        ix := NULL;
      end if;
      e := convert(e, 'radical', ix);
      if e :: specfunc(RootOf) then return on_fail(e) end if;
      return e;
    catch: return on_fail(e0); end try;
  end proc;

  # Try to find an eliminate (by evaluation, or simplification) integrals which
  # are free of `applyintegrand`s.
  elim_intsum := module ()
    export ModuleApply := proc(inert0, h :: name, kb :: t_kb, opts, $)
       local ex, e, inert := inert0;
       ex := extract_elim(inert, h, kb);
       e[0] := apply_elim(h, kb, ex);
       e[1] := check_elim(inert, e[0]);
       if e[1] = FAIL then inert
       else
         e[2] := reduce(e[1],h,kb,opts);
         if has(e[2], {csgn}) then
           WARNING("Throwing away an eliminated result result containing csgn (this "
                   "could be a bug!):\n%1\n(while running %2)", e[2], ex);
           inert;
         else e[2] end if;
       end if
    end proc;

    local known_tys := table([Int=int_assuming,Sum=sum_assuming,Ints=ints,Sums=sums]);
    local extract_elim := proc(e, h::name, kb::t_kb,$)
      local t, intapps, var, f, e_k, e_args, vs, blo, bhi;
      vs := {op(KB:-kb_to_variables(kb))};
      t := 'applyintegrand'('identical'(h), 'anything');
      intapps := indets(op(1,e), t);
      if intapps = {} then
        return FAIL;
      end if;
      e_k := op(0,e); e_args := op([2..-1],e);
      if Domain:-Has:-Bound(e) and assigned(known_tys[e_k]) then
        var := Domain:-ExtBound[e_k]:-ExtractVar(e_args);
        ASSERT(var::DomBoundVar);
        blo, bhi := Domain:-ExtBound[e_k]:-SplitRange
                    (Domain:-ExtBound[e_k]:-ExtractRange(e_args));
        if ormap(b->op(1,b) in map((q-> (q,-q)), vs) and op(2,b)::SymbolicInfinity
                ,[[blo,bhi],[bhi,blo]]) then
          return FAIL end if;
        if var :: list then var := op(1,var) end if;
        if not depends(intapps, var) then
          f := known_tys[e_k];
        else
          return FAIL;
        end if;
      end if;
      [ op(1,e), f, var, [e_args] ];
    end proc;

    local apply_elim := proc(h::name,kb::t_kb,todo::{list,identical(FAIL)})
      local body, f, var, rrest;
      if todo = FAIL then return FAIL; end if;
      body, f, var, rrest := op(todo);
      banish(body, h, kb, infinity, var,
             proc (kb1,g,$) do_elim_intsum(kb1, f, g, op(rrest)) end proc);
    end proc;

    local check_elim := proc(e, elim,$)
      if has(elim, {MeijerG, undefined, FAIL}) or e = elim or elim :: SymbolicInfinity then
        return FAIL;
      end if;
      return elim;
    end proc;

    local do_elim_intsum := proc(kb, f, ee, v::{name,name=anything})
      local w, e, x, g, t, r;
      w, e := op(Domain:-Extract:-Shape(ee));
      w := Domain:-Shape:-toConstraints(w);
      e := piecewise_And(w, e, 0);
      e := f(e,v,_rest,kb);
      x := `if`(v::name, v, lhs(v));
      g := '{sum, sum_assuming, sums}';
      if f in g then
        t := {'identical'(x),
              'identical'(x) = 'Not(range(Not({SymbolicInfinity, undefined})))'};
      else
        g := '{int, int_assuming, ints}';
        t := {'identical'(x),
              'identical'(x) = 'anything'};
        if not f in g then g := {f} end if;
      end if;
      for r in indets(e, 'specfunc'(g)) do
        if 1<nops(r) and op(2,r)::t then return FAIL end if
      end do;
      e
    end proc;
  end module; # elim
end module; # reduce_Integrals

int_assuming := proc(e, v::name=anything, kb::t_kb, $)
  simplify_factor_assuming('int'(e, v), kb);
end proc;

sum_assuming := proc(e, v::name=anything, kb::t_kb)
  simplify_factor_assuming('sum'(e, v), kb);
end proc;

# Int( .., var=var_ty ) == var &X var_ty
isBound_IntSum := kind -> module()
  option record;
  export MakeKB := (`if`(kind=Sum,KB:-genSummation,KB:-genLebesgue));
  export ExtractVar := (e->op(1,e));
  export ExtractRange := (e->op(2,e));
  export MakeRange := `..`;
  export SplitRange := (e->op(e));
  export Constrain := `if`(kind=Sum,`<=`,`<`);
  export DoMk := ((e,v,t)->kind(e,v=t));
  export Min := `min`; export Max := `max`;
  export VarType := 'name';
  export RangeType := 'range';
  export MapleType := 'And'('specfunc'(kind), 'anyfunc(anything,name=range)');
  export BoundType := `if`(kind=Sum,'integer','real');
  export RecogBound := `if`(kind=Sum,
            (proc(k,b)
               if   k = `<=` then (x->subsop(2=b,x))
               elif k = `>=` then (x->subsop(1=b,x))
               elif k = `<`  then (x->subsop(2=(b-1),x))
               elif k = `>`  then (x->subsop(1=b+1,x))
               end if;
             end proc),
            (proc(k,b)
               if   k in {`<=`,`<`} then (x->subsop(2=b,x))
               elif k in {`>=`,`>`} then (x->subsop(1=b,x))
               end if;
             end proc));
end module;

# Ints( .., var::name, var_ty::range, dims::list(name=range) ) ==
#        [ var   , map(lhs,dims) ] :: list(name)  &X
#        [ var_ty, map(rhs,dims) ] :: list(range)
isBound_IntsSums := kind -> module()
  option record;
  export MakeKB := proc(vars, lo, hi, kb, $)
    local var, dims, ty, rngs, x, kb1;
    var  := op(1, vars);
    rngs := zip(`..`,lo,hi);
    ty   := op(1, rngs);
    dims := subsop(1=NULL,zip(`=`,vars,rngs));
   x, kb1 := genType(var,
                      mk_HArray(`if`(kind=Ints,
                                     HReal(open_bounds(ty)),
                                     HInt(closed_bounds(ty))),
                                dims),kb);
    if nops(dims) > 0 then
      kb1 := assert(size(x)=op([-1,2,2],dims)-op([-1,2,1],dims)+1, kb1);
    end if;
    x, kb1;
  end proc;
  export ExtractVar   := ((v,t,d)->[v,map(lhs,d)[]]);
  export ExtractRange := ((v,t,d)->[t,map(rhs,d)[]]);
  export MakeRange    := ((a,b)->zip(`..`,a,b));
  export SplitRange   := (rs->(map(x->op(1,x),rs), map(x->op(2,x),rs)));
  export Constrain    := ((a,b)->zip(`if`(kind=Ints, `<`, `<=`),a,b)[]);
  export DoMk         := ((e,v,t)->kind( e,op(1,v),op(1,t), subsop(1=NULL,zip(`=`,v,t)) ));
  export Min          := ((a,b)->zip(`min`,a,b));
  export Max          := ((a,b)->zip(`max`,a,b));
  export VarType      := 'And(list(name),satisfies(x->x<>[]))';
  export RangeType    := 'And(list(range),satisfies(x->x<>[]))';
  export MapleType    := 'And'('specfunc'(kind),'anyfunc'('anything', 'name', 'range', 'list(name=range)'));
  export BoundType    := TopProp;
  export RecogBound   := (_->NULL);
end module;

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
with(Threads:-Task):
fun:=proc()
    Continue(passed,x,Task=[`+`,1,2,3],y,Tasks=[`*`,[1,2],[2,3,3]],z);
end proc:
Start(passed,x,Task=[fun],Tasks=[`+`,[2,3,3],[6,6,6]]);


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

$include "mgga_c_kcis.mpl"

kcisk_p := 0.193:
kcisk_gcnst := 20/(3*Pi*(6*Pi^2)^(1/3)):

(* Eqs. (27) *)
(* There was a "bug" in the reference implementation of Stefan Kurth. The
   f_pw was calculated for an unpolarized density instead of a polarized one.
   Here, we follow the original code. *)
kcisk_gamma0 := rs -> m_max(0, kcis_beta + 2^(1/3)*kcisk_gcnst*f_pw(2^(1/3)*rs, 0)/n_total(rs)^(1/3)):

kcisk_gamma1 := rs -> m_max(0, kcis_beta + kcisk_gcnst*f_pw(rs, 0)/n_total(rs)^(1/3)):

(* Eq. (19) and (22) *)
(* There seems to be a misspel in Eq. (22) regarding the "1 +". This follows the
   reference implementation of Stefan Kurth *)
kcis_gga0 := (rs, xt) -> f_pw(rs, 0) /
  (1 + kcisk_p*log(1 + kcisk_gamma0(rs)*kcis_t(rs, xt)^2/(kcisk_p*m_abs(f_pw(rs, 0))) )):

kcis_gga1 := (rs, xt) -> f_pw(rs, 1) /
  (1 + kcisk_p*log(1 + 2^(-1/3)*kcisk_gamma1(rs)*kcis_t(rs, xt)^2/(kcisk_p*m_abs(f_pw(rs, 1))) )):


text:
######################################################################
# single_cubes(k) is the usual set of injective maps [0,1]^k -> [0,1]^k

`is_element/single_cubes` := (k::posint) -> proc(xy)
 local i,x,y;
 global reason;

 if not(type(xy,[[realcons$k],[realcons$k]])) then
  reason := [convert(procname,string),"xy is not in R^k x R^k",xy];
  return false;
 fi;

 x,y := op(xy);
 for i from 1 to k do
  if not(0 <= x[i] and x[i] < y[i] and y[i] <= 1) then
   reason := [convert(procname,string),"not (0 <= x[i] < y[i] <= 1)",i,x[i],y[i]];
   return false;
  fi;
 od;

 return true;
end;

`is_equal/single_cubes` := (k::posint) -> proc(xy1,xy2)
 global reason;

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

 return true;
end:

`is_leq/single_cubes` := NULL;

`random_element/single_cubes` := (k::posint) -> proc(d := 6)
 local r,x,y,i,s,t;

 r := rand(0..d);
 x := NULL;
 y := NULL;
 for i from 1 to k do
  s := 0;
  t := 0;
  while s = t do s := r(); t := r(); od;
  x := x,min(s,t)/d;
  y := y,max(s,t)/d;
 od;
 return [[x],[y]];
end;

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

`act/single_cubes` := (k::posint) -> (xy) -> (t) ->
 xy[1] +~ (t *~ (xy[2] -~ xy[1]));

`id/single_cubes` := (k::posint) -> [[0$k],[1$k]];

`compose/single_cubes` := (k::posint) -> (xy,uv) ->
  map(`act/single_cubes`(k)(xy),uv);

`overlap/single_cubes` := (k::posint) -> proc(xy,uv)
 local x,y,u,v,i;

 x,y := op(xy);
 u,v := op(uv);
 for i from 1 to k do
  if max(x[i],u[i]) >= min(y[i],v[i]) then
   return false;
  fi;
 od;
 return true;
end;

`centre/single_cubes` := (k::posint) -> proc(xy)
 (xy[1] +~ xy[2]) /~ 2;
end;


text:
program change;  /* sample2a.mpl */ { compute change program. What is change? }
             { no 5, 50, 500, & 5000 yen used }
  var n, count : integer;
begin
  writeln('please input change');  readln(n);    count := 0;
while n>0 do begin
    while n>=10 do begin
        while n>=  100 do begin
           while n>=1000 do begin
             while n  >=10000 do begin
               count := count +1; n := n-10000;
             end;
             if count >0 then writeln('10000 yen : ',  count);
             count := 0;
             if n < 1000 then break;
             count := count +1; n := n-1000;
           end;
             if count >0 then writeln(' 1000 yen : ',  count);
             count := 0;
             if n < 100 then break;
              count := count +1; n := n-100;
         end;
         if count >0 then     writeln('  100 yen : ',  count);
             count := 0;
             if n < 10 then break;
              count := count +1; n := n-10;
      end;
         if count >0 then     writeln('   10 yen : ',  count);
             count := 0;
             if n < 1 then break;
              count := count +1; n := n-1;
 end;
 if count >            0 then writeln('    1 yen : ',  count);
/*  if n !=0 then break;  */
end.


text:
(*
    Test file for highlight
*)
`&+`:=proc(x::set,y::set)
    return x union y;
end proc:               # keywords
$include "utils.mpl"    # preprocessor
solve(x^2+x-1);         # maple functions
print(`&+`);            # quoted name
print("123");           # string
print('x+y');           # unvaluated expression
type(x,uneval);         # type name
fun:=proc(x::uneval)    # parameter modifier
    option inline;      # proc and module options
    print(procname);    # language variable
    return true;        # language constants
end proc:


text:
var $p <* i32>
func $Foo (var %i i32) void {
 iassign <* i32>(
    dread a64 $p,
    add i32(
     iread i32 <* i32> ( dread a64 $p ),
     dread i32 %i))
  return()
}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


text:
foo:bar;
|G G G G|G G G G|
[
	|F G G G|G|
	[
		|F G G F|
		|G G G F|
		|G G G F|
		|G D D D|
		|D E|
	]
	|F G F G|
	|F E F F|G|
	[ D ]
	|D E|
]
|F F C|F E E E|
|C G G G|G G G G|
|C C G G|G C F F|
|C D E C|D C G G|
|G C E E|E E E E|
|C E E E|E E E E|
|E C F F|G C F G|
|G C|

//+ = G
//- = E
//< = F
//> = D
//. = C
//

max_stars_repo_name: dwysocki/music-programming-language

max_issues_repo_name: dwysocki/music-programming-language

max_forks_repo_name: dwysocki/music-programming-language


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

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

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

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

dk87_betag := 7/(432*Pi*(6*Pi^2)^(1/3))/X_FACTOR_C:

dk87_f := x -> 1 + dk87_betag*x^2*(1 + params_a_a1*x^params_a_alpha)/(1 + params_a_b1*x^2):

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


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

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

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

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

$ifdef mgga_x_gx_params
params_a_c0       :=  0.827411:
params_a_c1       := -0.643560:
params_a_alphainf :=  0.852:
$endif

gx_alpha := (x, t) -> (t - x^2/8)/K_FACTOR_C:

gx_cx0 := 4/3*(2/Pi)^(1/3):
gx_cx1 := X_FACTOR_C:

gx_gx0 := a ->
  + gx_cx0/gx_cx1
  + a*(params_a_c0 + params_a_c1*a)/(1.0 + (params_a_c0 + params_a_c1 - 1)*a) * (1 - gx_cx0/gx_cx1):

gx_gx1 := a ->
  1 + (1 - params_a_alphainf)*(1 - a)/(1 + a):

gx_f_a := a->
  + gx_gx0(a)*Heaviside(1 - a)
  + gx_gx1(a)*Heaviside(a - 1):

gx_f := (x, u, t) ->
  gx_f_a(gx_alpha(x, t)):

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


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

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

(* type: work_mgga_x *)

a := 2.413:
b := 0.348:

params_a_b      := 0.40:
params_a_c      := 1.59096:
params_a_e      := 1.537:
params_a_mu     := 0.21951:

alpha           := (x, t) -> (t - x^2/8)/K_FACTOR_C:

(* Equation (8) *)

mkappa := (x, t) -> 2*Pi/(3*sqrt(5)) * \
          sqrt(alpha(x, t) + 1)/sqrt(a + log(alpha(x, t) + b)):
ff     := 2:

$include "tpss_x.mpl"

a1  := (x, t) -> mkappa(x, t)/(mkappa(x, t) + fx(x, t)):
f   := (rs, x, t, u) -> 1 + mkappa(x, t)*(1 - a1(x, t)):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
with(combinat, cartprod): ss := [`$`(0 .. ORDER)]:
cp := cartprod([ss, ss, ss, ss, ss, ss, ss]):
ret := []:
while not cp[finished] do
    ret := [op(ret), cp[nextvalue]()]
end do:
vals := select(proc (x) options operator, arrow; convert(x, `+`) <= ORDER; end proc, ret):
valsm1 := vals:
ret := 'ret';
# read "../vals.mpl":

DLIST := [dL, dS, dR, dLS, dLR, dRS, dLRS]:

xition := solve(
    [pLp=zLp + zLSp + zLRp + zLRSp,
     pSp=zSp + zLSp + zRSp + zLRSp,
     pRp=zRp + zLRp + zRSp + zLRSp,
     cLSp=(zLSp + zLRSp) - pLp*pSp,
     cLRp=(zLRp + zLRSp) - pLp*pRp,
     cRSp=(zRSp + zLRSp) - pSp*pRp,
     cLRSp=zLRSp - pLp*cRSp - pSp*cLRp - pRp*cLSp - pLp*pSp*pRp],
     [zLp, zSp, zRp, zLSp, zLRp, zRSp, zLRSp]
)[1]:

pL := zL + zLS + zLR + zLRS;
pS := zS + zLS + zRS + zLRS;
pR := zR + zLR + zRS + zLRS;
cLS := (zLS + zLRS) - pL*pS:
cLR := (zLR + zLRS) - pL*pR:
cRS := (zRS + zLRS) - pS*pR:
cLRS := zLRS - pL*cRS - pS*cLR - pR*cLS - pL*pS*pR:

a_S_0 := s * (h + (1 - 2 * h) * pS) / (1 + s * pS * (2 * h + (1 - 2 * h) * pS)):
a_S_S := (1 - 2 * h) * s / (1 + s * pS * (2 * h + (1 - 2 * h) * pS)):
qS := 1 - pS:
qL := 1 - pL:
qR := 1 - pR:

rLRS := rLS_R + rRS_L + rLR_S:

# Deltas
delta_pS := a_S_0 * pS * qS:
delta_pL := a_S_0 * cLS:
delta_pR := a_S_0 * cRS:

g := (r) -> -r + a_S_0 * (1 - r) * (1 - 2 * pS) + a_S_S * r * pS * qS;

deltat_cLS := g(rLS) * cLS:
deltat_cRS := g(rRS) * cRS:
deltat_cLR := -rLR * cLR + a_S_0 * (1 - rLR) * cLRS + a_S_S * rLR * cLS * cRS:
deltat_cLRS := (-rLRS + a_S_0 * (1 - rLRS) * (1 - 2 * pS) + a_S_S * rLR_S * pS * qS) * cLRS -
    a_S_0 * (rLS_R + rRS_L) * pS * qS * cLR -
    (a_S_0 * (2 - rLS_R - rRS_L) - a_S_S * (1 - 2 * pS) * (rLS_R + rRS_L)) * cLS * cRS:

delta_cLS := deltat_cLS - a_S_0^2 * pS * qS * cLS:
delta_cRS := deltat_cRS - a_S_0^2 * pS * qS * cRS:
delta_cLR := deltat_cLR - a_S_0^2 * cLS * cRS:
delta_cLRS := deltat_cLRS -
    a_S_0 * (pS * qS * deltat_cLR + cLS * deltat_cRS + cRS * deltat_cLS) +
    2 * a_S_0^3 * pS * qS * cLS * cRS:

B := eval(xition, [pLp=pL + delta_pL, pSp=pS + delta_pS, pRp=pR + delta_pR,
     cLSp=cLS + delta_cLS, cLRp=cLR + delta_cLR, cRSp=cRS + delta_cRS,
     cLRSp=cLRS + delta_cLRS]
):

if ARG0 + ARG1 + ARG2 + ARG3 + ARG4 + ARG5 + ARG6 = 0 then
CodeGeneration:-C(B, optimize=true, coercetypes=false, deducetypes=false, defaulttype=float, output="_transition.cii"):
end if:

mgff := (p1*exp(t1)+p2*exp(t2)+p3*exp(t3)+p4*exp(t4)+p5*exp(t5)+p6*exp(t6)+p7*exp(t7)+(1-p1-p2-p3-p4-p5-p6-p7))^n:
multimgf := proc (i1, i2, i3, i4, i5, i6, i7) options operator, arrow;
    s := i1 + i2 + i3 + i4 + i5 + i6 + i7:
    if s = 0 then 1
    else
        n^(-s)*simplify(subs(t1=0, t2=0, t3=0, t4=0, t5=0, t6=0, t7=0,
                             diff(mgff, `$`(t1, i1), `$`(t2, i2), `$`(t3, i3),
                                        `$`(t4, i4), `$`(t5, i5), `$`(t6, i6),
                                        `$`(t7, i7))));
    end if;
end proc:

mgfsubs := proc (iL, iS, iR, iLS, iLR, iRS, iLRS) options operator, arrow;
    eval(multimgf(iL, iS, iR, iLS, iLR, iRS, iLRS),
         [p1=zLp, p2=zSp, p3=zRp, p4=zLSp, p5=zLRp, p6=zRSp, p7=zLRSp]);
    eval(%, B):
end proc:

deltasubs := proc (iL, iS, iR, iLS, iLR, iRS, iLRS)
    ms := mgfsubs(iL, iS, iR, iLS, iLR, iRS, iLRS):
    eval(%, [zL=zL + dL, zS=zS + dS, zR=zR + dR,
             zLS=zLS + dLS, zLR=zLR + dLR, zRS=zRS + dRS,
             zLRS=zLRS + dLRS]):
end proc:

llhs := proc (aL, aS, aR, aLS, aLR, aRS, aLRS)
rary := []:
ms := deltasubs(aL, aS, aR, aLS, aLR, aRS, aLRS):
ind := 1:
for vv in valsm1 do
   print(vv):
   coeftayl(ms, DLIST = [0, 0, 0, 0, 0, 0, 0], vv):
   cg := eval(%, [zL=dtL, zS=dtS, zR=dtR, zLS=dtLS, zLR=dtLR, zRS=dtRS, zLRS=dtLRS]):
   rary := [op(rary), ret[ind] = cg * Edelta_lhs(t - 1, p, vv[1], vv[2], vv[3], vv[4], vv[5], vv[6], vv[7])]:
   ind := ind + 1:
end do:
return rary:
end proc:

rrhs := proc(aL, aS, aR, aLS, aLR, aRS, aLRS)
    local e, c, xx;
    xx := expand(
        (dtL + dL)^aL *
        (dtS + dS)^aS *
        (dtR + dR)^aR *
        (dtLS + dLS)^aLS *
        (dtLR + dLR)^aLR *
        (dtRS + dRS)^aRS *
        (dtLRS + dLRS)^aLRS
        - dL^aL * dS^aS * dR^aR * dLS^aLS * dLR^aLR * dRS^aRS * dLRS^aLRS);
    c := coeffs(xx, DLIST, 'e');
    convert(zip(proc (ee, cc)
        options operator, arrow;
        degs := map((term) -> degree(ee, term), DLIST):
        Edelta(t, p, degs[1], degs[2], degs[3], degs[4], degs[5], degs[6], degs[7])*cc;
    end proc, [e], [c]), `+`);
end proc:

aL := ARG0;
aS := ARG1;
aR := ARG2;
aLS := ARG3;
aLR := ARG4;
aRS := ARG5;
aLRS := ARG6;
print(aL, aS, aR, aLS, aLR, aRS, aLRS);
ll := llhs(aL, aS, aR, aLS, aLR, aRS, aLRS):
outl := cat("_lhs_", aL, aS, aR, aLS, aLR, aRS, aLRS, ".cii"):
CodeGeneration:-C(ll, optimize=true, coercetypes=false, deducetypes=false, defaulttype=float, output=outl):

rr := rrhs(aL, aS, aR, aLS, aLR, aRS, aLRS):
outl := cat("_rhs_", aL, aS, aR, aLS, aLR, aRS, aLRS, ".cii"):
CodeGeneration:-C([ret = rr], optimize=true, coercetypes=false, deducetypes=false, defaulttype=float, output=outl):

mgfsubs(aL, aS, aR, aLS, aLR, aRS, aLRS):
mm := eval(%, [zL=dtL, zS=dtS, zR=dtR, zLS=dtLS, zLR=dtLR, zRS=dtRS, zLRS=dtLRS]):
outl := cat("_mgf_", aL, aS, aR, aLS, aLR, aRS, aLRS, ".cii"):
CodeGeneration:-C([ret = mm], optimize=true, coercetypes=false, deducetypes=false, defaulttype=float, output=outl):

max_stars_repo_name: terhorst/EandR-timeseries

max_issues_repo_name: terhorst/EandR-timeseries

max_forks_repo_name: terhorst/EandR-timeseries


text:
with(Hakaru):
with(NewSLO):
with(Partition):
with(KB):

TestHakaru(Partition(a<b, VAL, b>a, VAL), VAL, label="Partition with identical bodies");

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
export Shape := module ()
    export toConstraints := proc(sh_ :: DomShape, $)
           ::specfunc({boolean,relation,specfunc(`Or`)}, `And`);
        local sh := sh_;
        if sh :: specfunc(`DConstrain`) then
            And( op(sh) );
        elif sh :: specfunc(`DSum`) then
            sh := Or(op(map(toConstraints, sh)));
            Domain:-simpl_relation(sh, norty='CNF');
        elif sh :: specfunc(`DInto`) then
            toConstraints(op(3, sh));
        else
            error "don't know how to convert to constraints %1", 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 *)

$include "vwn.mpl"

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

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


text:
type $RR <class {@e1 i32, @e2 f32, @e3 f64}>
type $SS <class <$RR> {@f1 i32, @f2 i8, @f3 i8}>
type $TT <class {@e1 i32, @e2 f32, @e3 f64}>
type $UU <class <$RR> {@f1 i32, @f2 i8, @f3 i8}>
func &foo ( var %x <$SS>) void {
  var %y <$UU>
  dassign %y 0 ( dread agg %x 0 )
  return () }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
# Machine Learning

##
## Copyright (C) 2018 Maplesoft
## Authors:    Behzad Samadi <bsamadi@maplesoft.com>
## Created:    January 2018
## Version:    0.1
## Keywords:   Machine Learning, Maple, Code Generation, Reinforcement Learning, Gaussian Process Regression, Bayesian Optimizaiton
##
## Procedures:
##    Simulate
##

MachineLearning := module()
description "Machine learning";
option package;

export CostFunction,
       CovarianceMatrix,
       DotProduct,
       KernelFunction,
       Mean,
       Mean1D,
       Norm2,
       PDMatrixInverse,
       StandardDeviation,
       StandardDeviation1D,
       SquaredDistance,
       Simulate,
       SolveLT,
       SolveUT;

$include "Simulation/Simulate.mm"
$include "LinearAlgebra/DotProduct.mm"
$include "LinearAlgebra/Norm2.mm"
$include "LinearAlgebra/PDMatrixInverse.mm"
$include "LinearAlgebra/SolveLT.mm"
$include "LinearAlgebra/SolveUT.mm"
$include "GaussianProcess/CostFunction.mm"
$include "GaussianProcess/CovarianceMatrix.mm"
$include "GaussianProcess/KernelFunction.mm"
$include "GaussianProcess/Mean.mm"
$include "GaussianProcess/SquaredDistance.mm"
$include "GaussianProcess/StandardDeviation.mm"

end module; # MachineLearning

max_stars_repo_name: SymbolicNMPC/Machine-Learning

max_issues_repo_name: SymbolicNMPC/Machine-Learning

max_forks_repo_name: SymbolicNMPC/Machine-Learning


text:
make_alternating_five := proc()
 global alternating_five;
 local G,x2,y2,z2,x3,x5,y5,id,o,oo,inv,idm,om,oom,E,ccl,
       t0,t1,i,j,k,l,g,h,w,m,x,u0,u1,u2,u3,u,v,M,theta;

 G := table():

 x2 := [2,1,4,3,5];
 y2 := [3,4,1,2,5];
 z2 := [1,3,2,5,4];
 x3 := [2,3,1,4,5];
 x5 := [3,5,4,2,1];
 y5 := [2,3,4,5,1];

 id := [1,2,3,4,5];
 o  := (u,v) -> [seq(u[v[i]],i=1..5)];
 oo := proc()
  if   nargs = 0 then return(id);
  elif nargs = 1 then return(args[1]);
  elif nargs = 2 then return(o(args));
  else return o(args[1],oo(args[2..-1]));
  fi;
 end:
 inv := proc(g)
  local gi;
  gi := table():
  for i from 1 to 5 do gi[g[i]] := i; od;
  return [seq(gi[i],i=1..5)];
 end:

 G["id"] := id;
 G["o"]  := eval(o);
 G["oo"]  := eval(oo);
 G["inv"] := eval(inv);

 G["x2"] := x2;  G["x2c"] := convert(x2,disjcyc);
 G["y2"] := y2;  G["y2c"] := convert(y2,disjcyc);
 G["z2"] := z2;  G["z2c"] := convert(z2,disjcyc);
 G["x3"] := x3;  G["x3c"] := convert(x3,disjcyc);
 G["x5"] := x5;  G["x5c"] := convert(x5,disjcyc);
 G["y5"] := y5;  G["y5c"] := convert(y5,disjcyc);

 G["H"] := table():

 G["H"][ 1] := [id];
 G["H"][ 2] := [id,x2];
 G["H"][ 3] := [id,x3,o(x3,x3)];
 G["H"][ 4] := [id,x2,y2,o(x2,y2)];
 G["H"][ 5] := [id,x5,o(x5,x5),oo(x5 $ 3),oo(x5 $ 4)];
 G["H"][ 6] := [seq(seq(o(u,v),u in [id,z2]),v in G["H"][3])];
 G["H"][10] := [seq(seq(o(u,v),u in G["H"][ 2]),v in G["H"][5])];
 G["H"][12] := [seq(seq(o(u,v),u in G["H"][ 4]),v in G["H"][3])];
 G["H"][60] := [seq(seq(o(u,v),u in G["H"][12]),v in G["H"][5])];

 E := G["H"][60];
 G["elements"] := E;

 G["element_index"] := table();
 for i from 1 to 60 do
  G["element_index"][E[i]] := i;
 od:

 G["gens_rule"] := {X2 = x2, Y2 = y2, Z2 = z2, X3 = x3, X5 = x5, Y5 = y5};

 # words in x2, y2, z2 for all elements
 G["words_a"] := [
 [],[X2],[Y2],[X2,Y2],[X2,Z2,Y2,Z2,X2,Z2],[Z2,Y2,Z2,X2,Z2],[Z2,Y2,Z2,X2,Z2,X2],
 [Y2,Z2,Y2,Z2,X2,Z2],[Y2,Z2,X2,Z2,Y2,Z2],[Z2,X2,Z2,Y2,Z2,Y2],[Z2,X2,Z2,Y2,Z2],
 [X2,Z2,X2,Z2,Y2,Z2],[X2,Y2,Z2,Y2],[Y2,Z2,Y2],[X2,Z2,Y2],[Z2,Y2],[X2,Y2,Z2,Y2,Z2,X2,Y2],
 [Y2,Z2,Y2,Z2,X2,Y2],[X2,Z2,Y2,Z2,X2,Y2],[Z2,Y2,Z2,X2,Y2],[X2,Y2,Z2,X2,Z2],[Y2,Z2,X2,Z2],
 [X2,Z2,X2,Z2],[Z2,X2,Z2],[X2,Y2,Z2,X2,Z2,Y2],[Y2,Z2,X2,Z2,Y2],[X2,Z2,X2,Z2,Y2],
 [Z2,X2,Z2,Y2],[X2,Y2,Z2],[Y2,Z2],[X2,Z2],[Z2],[X2,Y2,Z2,Y2,Z2,X2],[Y2,Z2,Y2,Z2,X2],
 [X2,Z2,Y2,Z2,X2],[Z2,Y2,Z2,X2],[Z2,Y2,Z2,Y2,Z2,X2,Z2],[X2,Z2,Y2,Z2,Y2,Z2,X2,Z2],
 [Z2,X2,Z2,X2,Y2],[X2,Z2,X2,Z2,X2,Y2],[Z2,X2],[X2,Z2,X2],[Y2,Z2,X2],[X2,Y2,Z2,X2],
 [X2,Z2,Y2,Z2],[Z2,Y2,Z2],[X2,Y2,Z2,Y2,Z2],[Y2,Z2,Y2,Z2],[Y2,Z2,X2,Y2],[Z2,X2,Y2,Z2],
 [Z2,X2,Y2],[X2,Z2,X2,Y2],[Z2,Y2,Z2,Y2],[X2,Z2,Y2,Z2,Y2],[Y2,Z2,Y2,Z2,Y2],
 [X2,Y2,Z2,Y2,Z2,Y2],[X2,Z2,X2,Z2,X2],[Z2,X2,Z2,X2],[X2,Y2,Z2,X2,Z2,X2],[Y2,Z2,X2,Z2,X2]];

 # words in x2, y2, x3, x5 for all elements
 G["words_b"] := [seq(seq(seq(seq([X2$i,Y2$j,X3$k,X5$l],i=0..1),j=0..1),k=0..2),l=0..4)];

 G["rels"] := [
  [X2$2],[Y2$2],[Z2$2],[X3$3],[X5$5],[Y5$5],
  map(op,[[X2,Y2]$2]),map(op,[[X2,Y2,Z2]$3]),
  map(op,[[Z2,X2,Z2,Y2]$3]),map(op,[[X2,Z2]$5]),map(op,[[Y2,Z2]$5]),
  [X3,Y2,Z2,X2,Z2,Y2,Z2],[Y5,Y2,Z2,X2,Z2,X2],[X5,Y2,X5,X5,X5,X3,X3]];

 ccl := (x) -> sort([op({seq(oo(g,x,inv(g)), g in E)})]);

 G["conjugacy_classes"] := [
  ccl(id),
  ccl(x2),
  ccl(x3),
  ccl(o(x5,x5)),
  ccl(x5)
 ];

 G["conjugacy_class_index"] := table();
 for i from 1 to 5 do
  for g in G["conjugacy_classes"][i] do
   G["conjugacy_class_index"][g] := i;
  od;
 od;

 G["element_order"] := table();

 for g in E do
  h := g;
  i := 1;
  while h <> id do i := i + 1; h := o(g,h); od:
  G["element_order"][g] := i;
 od:

 t0 := (sqrt(5)+1)/2;
 t1 := (sqrt(5)-1)/2;

 # Dodecahedron representation

 idm := [[1,0,0],[0,1,0],[0,0,1]];
 om := (u,v) -> expand(convert(Matrix(u) . Matrix(v),listlist));
 oom := proc()
  if   nargs = 0 then return(idm);
  elif nargs = 1 then return(args[1]);
  elif nargs = 2 then return(om(args));
  else return om(args[1],oom(args[2..-1]));
  fi;
 end:

 G["idm"] := idm;
 G["om"]  := eval(om);
 G["oom"] := eval(oom);

 G["matrix_gens_rule"] := {
   X1 = [[ 1,0,0],[0, 1,0],[0,0, 1]],
   X2 = [[ 1,0,0],[0,-1,0],[0,0,-1]],
   Y2 = [[-1,0,0],[0,-1,0],[0,0, 1]],
   Z2 = [[t1,-t0,-1],[-t0,-1,t1],[-1,t1,-t0]] /~ 2,
   X3 = [[0,0,1],[1,0,0],[0,1,0]],
   X5 = [[t1,-t0,1],[t0,1,t1],[-1,t1,t0]] /~ 2,
   Y5 = [[t0,1,-t1],[-1,t1,-t0],[-t1,t0,1]] /~ 2
 };

 G["matrix"] := table();
 G["axis"]   := table();
 G["angle"]  := table();

 for i from 1 to 60 do
  g := E[i];
  G["matrix"][g] := oom(op(subs(G["matrix_gens_rule"],G["words_b"][i])));
  if i = 1 then
   G["axis"][g] := [0,0,0];
   G["angle"][g] := 0;
  else
   M := Matrix(G["matrix"][g]);
   u0 := convert(NullSpace(M - IdentityMatrix(3))[1],list);
   u0 := simplify(expand(rationalize(u0 /~ sqrt(add(u0[j]^2,j=1..3)))));
   u1 := [1,2,3];
   u1 := u1 -~ add(u1[j] * u0[j],j=1..3) *~ u0;
   u1 := simplify(expand(rationalize(u1)));
   u2 := [u0[2] * u1[3] - u0[3] * u1[2],
          u0[3] * u1[1] - u0[1] * u1[3],
          u0[1] * u1[2] - u0[2] * u1[1]];
   u2 := simplify(expand(rationalize(u2)));
   u3 := simplify(expand(rationalize(convert(M . Vector(u1), list))));

   theta := evalf(arctan(add(u3[j] * u2[j],j=1..3),add(u3[j] * u1[j],j=1..3)));
   theta := round(evalf(60 * theta / Pi)) * Pi / 60;

   if evalf(theta) < 0 then u0 := -~ u0; theta := -theta; fi;
   G["axis"][g] := u0;
   G["angle"][g] := theta;
  fi;
 od;

 alternating_five := eval(G);
 return eval(G);
end:

make_alternating_five():

check_alternating_five := proc()
 local A5,f,fm,is_rot,ok,i,j,ij,g,h,gh,gm,hm,ghm0,ghm1;

 A5 := eval(make_alternating_five()):

 f := (u) -> A5["oo"](op(subs(A5["gens_rule"],u)));
 fm := (u) -> A5["oom"](op(subs(A5["matrix_gens_rule"],u)));

 _ASSERT(
  evalb({op(map(f,A5["rels"]))} = {A5["id"]}),
  "permutation relations");

 _ASSERT(
  evalb({op(map(fm,A5["rels"]))} = {A5["idm"]}),
  "matrix relations");

 _ASSERT(
  `and`(seq(evalb(f(A5["words_a"][i]) = A5["elements"][i]),i=1..60)),
  "words_a");

 _ASSERT(
  `and`(seq(evalb(f(A5["words_b"][i]) = A5["elements"][i]),i=1..60)),
  "words_b");

 is_rot := proc(m)
  local M,N;
  M := Matrix(m);
  if (simplify(expand(Determinant(M) - 1)) <> 0) then
   return false;
  fi;

  N := map(simplify,map(expand,Transpose(M) . M - IdentityMatrix(3)));
  if not(Equal(N,Matrix(3,3))) then return false; fi;
  return true;
 end:

 _ASSERT(
  `and`(seq(is_rot(A5["matrix"][g]),g in A5["elements"])),
  "matrices are rotations"
 );

 ok := true;
 for i from 1 to 60 do
  for j from 1 to 60 do
   g := A5["elements"][i];
   h := A5["elements"][j];
   gh := A5["o"](g,h);
   ij := A5["element_index"][gh];
   gm := A5["matrix"][g];
   hm := A5["matrix"][h];
   ghm0 := A5["matrix"][gh];
   ghm1 := A5["om"](gm,hm);
   if ghm0 <> ghm1 then ok := false; break; fi;
  od:
  if not(ok) then break; fi;
 od:

 _ASSERT(ok,"matrix representation is a homomorphism");
end:


text:
// A COMMEnt
/* Another
   comment */
var X : int := 4 + (6 * 2);
print X;
print "\n";

max_stars_repo_name: nicohi/mini-pl-interpreter

max_issues_repo_name: nicohi/mini-pl-interpreter

max_forks_repo_name: nicohi/mini-pl-interpreter