text:
#
# Erstelle Liste aller Kinematikparameter
# Init
# Dieses Arbeitsblatt erstellt einen Vektor mit allen Kinematikparameter
# Es werden nur die Kinematikparameter eingetragen, die ungleich Null sind, also z.B. keine MDH-Parameter die Länge oder Winkel Null sind.
#
# TODO: Liste mit Kinematikparametern muss immer die selbe Reihenfolge haben. Am besten in der Form a,alpha,qoffset,d,b,beta,kc. Momentan scheint die Reihenfolge zufällig, aber reproduzierbar.
#
# Moritz Schappler, schappler@irt.uni-hannover.de, 2017-02
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
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=100): # Damit der ganze Parametervektor ausgegeben werden kann
;
with(LinearAlgebra):
read "../helper/proc_MatlabExport":
read "../helper/proc_convert_t_s":
codegen_act := true:
# Lese Umgebungsvariable für Codegenerierung.
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_twist_definitions", robot_name):
robot_name_OL := robot_name: # Zusätzliche Variable zur Abgrenzung für implizite Zwangsbedingungen
;
#  Schalte um zwischen impliziten ZB und normalen Systemen
if FileTools[Exists]("../workdir/tbmode") then
  read "../workdir/tbmode": # Datei tbmode wird von Bash-Skripten erstellt und zeigt den aktuellen Modus an
else
  tbmode := "serial":
end if:
if tbmode = "implicit" then
  printf("Modus für IC ist aktiv. Bestimme Kinematikparameter für System mit impliziten ZB.\n"):
  read "../robot_codegen_definitions/robot_env_IC":
end if:
printf("Generiere Kinematik-Parametervektor für %s\n",robot_name):
# Parameter der Zwangsbedingungen lesen
kin_constraints_exist := false:
kc_symbols2 := []:
# Für explizite Zwangsbedingungen
constrfile := sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name_OL):
if FileTools[Exists](constrfile) then
  read constrfile:
  printf("Symbole der expliziten Zwangsbedingungen aus %s gelesen.\n", constrfile):
  kc_symbols2 := kc_symbols:
  kin_constraints_exist := true:
end if:
# Parameter für implizite Zwangsbedingungen (dadurch können auch neue Konstanten hinzugefügt werden)
constrfile := sprintf("../codeexport/%s/tmp/kinematic_implicit_constraints_symbols_list_maple", robot_name):
# Für implizite Zwangsbedingungen ist "robot_name" das System mit IC und nicht mit OL o.ä.
if FileTools[Exists](constrfile) then
  read constrfile:
  printf("Symbole der impliziten Zwangsbedingungen aus %s gelesen.\n", constrfile):
  kc_symbols2 := <kc_symbols2| kc_symbols>:
  kin_constraints_exist := true:
end if:
# Liste der Kinematikparameter erstellen
# Alle einzelnen MDH-Parametervektoren stapeln
pkin_tmp1 := <a;alpha;d;theta; qoffset;b;beta>:
# Kinematikparameter für kinematische Zwangsbedingungen hinzufügen (falls vorhanden)
if kin_constraints_exist then
  pkin_tmp1 := <pkin_tmp1;Transpose(kc_symbols2)>:
  np_kc := ColumnDimension(kc_symbols2):
else
  np_kc := 0:# Anzahl der Kinematikparameter für die Zwangsbedingungen
end if:
# Alle Symbole herausfinden
# Durch den Befehl fallen die vektoriell angesprochenen Gelenkkoordinaten (qJ) weg.
# TODO: Liste mit besserer Reihenfolge (siehe Dateikopf)
nms:=convert(indets(pkin_tmp1,name),list):
# Laufende Nummer für Parameterliste
kk := 0:
# Matrix mit Platz für alle Kinematikparameter (falls alle ungleich Null)
pkin_tmp2 := Matrix(NJ*6+np_kc,1):
# Terme aus Liste entfernen
for i from 1 to ColumnDimension(nms) do
  if nms[i] = Pi then # (wird im Code sowieso automatisch eingesetzt)
    next:
  end if:
  if nms[i] = t then # sollte eigentlich gar nicht drin sein können
    next:
  end if:
  # Parameter zur Liste hinzufügen
  kk := kk + 1:
  pkin_tmp2[kk,1] := nms[i]:
end do:
# Konstante Winkel, die aus MDH-Definition und nicht aus den ZB kommen umwandeln.
q_s := Matrix(1,1,0): q_t :=q_s:
pkin_tmp1 := convert_t_s(pkin_tmp1):
# Konvertierung nochmals durchführen (damit nicht delta8 und delta8s als verschiedene Symbole auftauchen)
pkin_tmp2 :=convert_t_s(pkin_tmp2):
# Ausgabevariable belegen
pkin := Matrix(pkin_tmp2[1..kk,1],kk,1):
# Symbole nochmals gruppieren
pkin := Transpose(Matrix(convert(indets(pkin,name),list))):
# Ergebnis speichern
# Für Generierung des Kinematikparametervektors in Matlab
MatlabExport(pkin, sprintf("../codeexport/%s/tmp/parameter_kin_matlab.m", robot_name), 2):
# Für schnelle Erkennung der Dimension zum Auslesen durch Bash-Skripte
save pkin, sprintf("../codeexport/%s/tmp/parameter_kin", robot_name):
# Zuordnung zwischen Kinematikparametern und MDH-Parametern
# Mit der Zuordnung kann der Vektor der Kinematikparameter (der an jede Funktion übergeben wird) aus den MDH-Parametern bestimmt werden
# TODO: Dieser Ansatz funktioniert nur, wenn keine Additionen von Variablen in den MDH-Parametern enthalten sind.
# Schritt 1: pkin_tmp1 nachverarbeiten: Gelenkvariablen und Pi entfernen (damit der solve-Befehl weniger zu verarbeiten hat)
pkin_tmp3 := pkin_tmp1: # Variable, die gestapelt eine nachbearbeitete Version aller MDH-Parameter enthält
;
for i from 1 to RowDimension(pkin_tmp3) do
  # q entfernen
  for j from 1 to RowDimension(qJ_t) do
    pkin_tmp3(i,1) := subs({qJ_t(j,1)=0}, pkin_tmp3(i,1)):
  end do:
  # Term Null setzen, wenn keine Variablen drin stehen (z.B. nur Pi)
  nms:=convert(indets(pkin_tmp3(i),name),list):
  var_da := 0: # Zähler, ob Variablen enthalten sind
  for k from 1 to ColumnDimension(nms) do
    if nms[k] = Pi then # Pi ist Konstante und funktioniert später nicht im Solve-Befehl
      next:
    end if:
    var_da := 1:
  end do:
  if var_da = 0 then
    pkin_tmp3(i,1) := 0:
  end if:
end do:

pkin_subs_mdh:=copy(pkin): # ohne "copy" wird anscheinend verlinkt und der Rest geht nicht
;
pkt3i := Matrix(RowDimension(pkin_tmp3), 1): # Index-Vektor für pkin_tmp3
;
for i from 1 to RowDimension(pkin_tmp3) do
  pkt3i(i,1):=parse(sprintf("pp%03d", i)): # Hilfsgröße zum späteren Ersetzen
end do:
# Ersetze

for i from 1 to RowDimension(pkin_tmp3) do
  if pkin_tmp3(i) = 0 then
    next:
  end if:
  for j from 1 to RowDimension(pkin) do
    # printf("i=%d, j=%d\n", i, j):
    erg := solve({pkin_tmp3(i)=pkt3i(i)}, pkin(j)):
    if ArrayTools:-NumElems(erg) = 0 then
      next:
    end if:
    # printf("i=%d\n", i):
    pkin_subs_mdh(j) := subs({lhs(erg[1])=rhs(erg[1])}, pkin_subs_mdh(j)):
    break: # Der erste Fund reicht
  end do:
end do:

# Ergebnis nachverarbeiten
# Platzhalter ("ph") generieren (für die Eingabeargumente der Matlab-Funktion zur Berechnung der Parameter)
beta_ph := Matrix(NJ,1):
b_ph := Matrix(NJ,1):
alpha_ph := Matrix(NJ,1):
a_ph := Matrix(NJ,1):
theta_ph := Matrix(NJ,1):
d_ph := Matrix(NJ,1):
qoffset_ph := Matrix(NJ,1):
for i from 1 to NJ do
  beta_ph(i) := parse(sprintf("beta_mdh(%d)", i)):
  b_ph(i) := parse(sprintf("b_mdh(%d)", i)):
  alpha_ph(i) := parse(sprintf("alpha_mdh(%d)", i)):
  a_ph(i) := parse(sprintf("a_mdh(%d)", i)):
  theta_ph(i) := parse(sprintf("theta_mdh(%d)", i)):
  d_ph(i) := parse(sprintf("d_mdh(%d)", i)):
  qoffset_ph(i) := parse(sprintf("qoffset_mdh(%d)", i)):
end do:
# Platzhalter-Vektor zur Ersetzung der gestapelten MDH-Parameter
pkt1_ph := <a_ph;alpha_ph;d_ph;theta_ph; qoffset_ph;b_ph;beta_ph>:
# Kinematikparameter für kinematische Zwangsbedingungen hinzufügen (falls vorhanden)
#TODO (aktuell wird das Problem mit der NaN-Ersetzung unten umgangen).
;
# Parameter-Indizes mit den Platzhaltern für die Matlab-Funktion ersetzen
for i from 1 to RowDimension(pkin_subs_mdh) do
  for j from 1 to RowDimension(pkt1_ph) do
    pkin_subs_mdh(i) := subs( {pkt3i(j)=pkt1_ph(j)}, pkin_subs_mdh(i) ):
  end do:
end do:
# Ersetze alle noch nicht ersetzten Kinematikparameter mit NaN. Ansonsten ist der Matlab-Code nicht lauffähig und es treten nachgelagerte Probleme auf.
for i from 1 to RowDimension(pkin_subs_mdh) do
  for j from 1 to RowDimension(pkt3i) do
    if has(pkin_subs_mdh(i), pkt3i(j)) then
      printf("Die Variable %s kommt noch in der Umrechnung mdh->pkin vor. Vermutlich Fehler beim substituieren. Setze NaN.\n", String(pkt3i(j))):
      pkin_subs_mdh(i) := subs( {pkt3i(j)=NaN}, pkin_subs_mdh(i) ):
    end if:
  end do:
  for j from 1 to RowDimension(kintmp_s) do
    if has(pkin_subs_mdh(i), kintmp_s(j)) then
      printf("Die Variable %s kommt noch in der Umrechnung mdh->pkin vor. Vermutlich Fehler beim substituieren. Setze NaN.\n", String(kintmp_s(j))):
      pkin_subs_mdh(i) := subs( {kintmp_s(j)=NaN}, pkin_subs_mdh(i) ):
    end if:

  end do:
end do:

# Folgender Term muss immer erzeugt werden (wird für Bash-Skripte benötigt. Sonst dort Fehler und Abbruch).
interface(warnlevel=0): # Unterdrücke die folgende Warnung (weil MDH-Parameter als Funktionsnamen interpretiert werden. Code funktioniert trotzdem.)
MatlabExport(pkin_subs_mdh, sprintf("../codeexport/%s/tmp/parameter_kin_from_mdh_matlab.m", robot_name), 2):
interface(warnlevel=3):

# Ausgabe
# MDH-Tabelle ausgeben
interface(rtablesize=100):
Test := <<seq(i, i=1..NJ)> | sigma| mu| beta| b | alpha | a | theta | d | v>:
Test:=<<"i" | "sigma" | "mu"|"beta"|"b"|"alpha"|"a"|"theta"|"d"|"v">,Test>:
printf("MDH-Tabelle für %s:\n", robot_name):
Test;
printf("Kinematik-Parameter für %s: %dx%d\n", robot_name, RowDimension(pkin), ColumnDimension(pkin)):
Transpose(pkin);

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


text:
bc :=
[
    [
        [2-7*x1+x1^2*x2-1/2*(x3-x1),
         6*x1-x1^2*x2-5*(x4-x2),
         2-7*x3+x3^2*x4-1/2*(x1-x3),
         6*x3-x3^2*x4+1+1/2*(x2-x4)
        ],
        [x4,x3,x2,x1]
    ],
    [
        [1-c*x-x*y^2-x*z^2,
         1-c*y-y*x^2-y*z^2,
         1-c*z-z*x^2-z*y^2,
         8*c^6+378*c^3-27
        ],
        [z, y, x, c]
    ],
    [
        [a+b+c+d+h1,
         a*b+b*c+c*d+d*h1+h1*a,
         a*b*c+b*c*d+c*d*h1+d*h1*a+h1*a*b,
         a*b*c*d+b*c*d*h1+c*d*h1*a+d*h1*a*b+h1*a*b*c,
         a*b*c*d*h1-1
        ],
        [h1, d, c, b, a]
    ],
    [
        [2*x1*(2-x1-y1)+x2-x1,
         2*x2*(2-x2-y2)+x1-x2,
         2*y1*(5-x1-2*y1)+y2-y1,
         y2*(3-2*x2-4*y2)+y1-y2
        ],
        [x2, x1, y2, y1]
    ],
    [
        [1/100-4*s*(s-1)*(s-b)*(s-c),
         1/5-b*c,
         2*s-1-b-c
        ],
        [s, c, b]
    ],
    [
        [x^2+y^2-x*y-1,
         y^2+z^2-y*z-a^2,
         z^2+x^2-z*x-b^2,
         a^2-1+b-b^2,
         3*b^6+56*b^4-122*b^3+56*b^2+3
        ],
        [z, y, x, b, a]
    ]
]:

ext: mpl

max_stars_repo_name: lihaokun/StrongSfTriDec

max_issues_repo_name: lihaokun/StrongSfTriDec

max_forks_repo_name: lihaokun/StrongSfTriDec


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

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

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

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

(* beta is taken from the params *)
params_a_gamma := (1 - log(2))/Pi^2:
params_a_BB    := 1:
$include "gga_c_pbe.mpl"

$include "tpss_c.mpl"

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


text:
#
# These procedures come from the source code of GFun.
#
getname:=proc(yofz::function(name), y, z)
  y:=op(0,yofz);
  if type(y,'procedure') then error `not an unassigned name`,y fi;
  z:=op(yofz)
end proc:

#
# returns the smallest i such that u(n+i) appears in a recurrence
#
minindex := proc(rec,u,n)
  min(op(map(op,indets(rec,'specfunc'('linear'(n),u)))))-n
end proc:

#
# returns the largest i such that u(n+i) appears in a recurrence
#
maxindex := proc(rec,u,n)
  max(op(map(op,indets(rec,'specfunc'('linear'(n),u)))))-n
end proc:

#
# A recurrence of the form a(n+d) = p(n)/q(n) a(n) is represented through a record:
# OneTermRecurrence : record(order, numerator, denominator)
#
`type/OneTermRecurrence` := 'record(order, numerator, denominator)':

#
#checkOneTermRecurrence
# Input: a recurrence rec (either with or without initial conditions).
#          If it has initial conditions, they are ignored.
#        a(n): the name of the sequence and the name of the variable.
#
# Output:
# This procedure checks that rec is a recurrence of the form a(n+d) = p(n)/q(n) a(n)
# If the check succeeds, it returns the corresponding record. If it fails, an error is
# returned.
#
checkOneTermRecurrence := proc(rec, aofn)::OneTermRecurrence;
  local r, d, a, n, term1, term2, res;

  getname(aofn, a, n):
  if type(rec, 'set') then
    r:=select(has, rec, n);
    if nops(r)>1
      then error `invalid recurrence`, rec
    fi:
    if nops(r)=0
      then error "%1 does not appear in the recurrence", n
    fi:
    r := op(r):
  else r:=rec:
  fi:
  if type(r,'`=`')
    then r:=op(1,r)-op(2,r)
  fi:
  if indets(r,'specfunc'('anything',a)) <> indets(r,'specfunc'('linear'(n),a))
    then error "the recurrence contains elements that are not linear in %1", n
  fi:
  if nops(r) <> 2
    then error "the recurrence contains %1 terms (expected 2)", nops(r)
  fi:
  r := subs(n=n-minindex(r, a, n), r):
  d := maxindex(r, a, n):

  term1 := select(has, r, a(n)):
  term2 := select(has, r, a(n+d)):

  res := factor( -(term1/a(n)) / (term2/a(n+d)) ):

  Record( 'order'=d, 'numerator' = numer(res), 'denominator' = denom(res) )
end proc:

#
# my_factors factorizes p the same way as factors(p) would do except that the constant part is computed
# differently. We assume here that p has integer coefficients, and we want to factorize it over polynomials
# with integer coefficients. my_factors ensures that the factors have integer coefficients.
#
my_factors := proc(p)
  local L, c, fact, i, my_c, my_fact, q:
  L := factors(p):
  c := L[1]: fact := L[2]:
  my_c := c: my_fact := []:
  for i from 1 to nops(fact) do
    q := denom(fact[i][1]):
    my_fact := [ op(my_fact), [ fact[i][1]*q, fact[i][2] ] ]:
    my_c := my_c / (q^fact[i][2]):
  od:
  [ my_c, my_fact]:
end proc:

#
# This procedure decomposes a one-term recurrence with the following form:
# a(n+d) = c * s1(n)/s1(n+d) * s2(n+d)/s2(n) * p(n)/q(n) * a(n)
#
# Known issue: this procedure assumes that the only variables involved are n and x with their usual meaning.
#
decomposeOneTermRecurrence := proc(formalRec::OneTermRecurrence, res_cste, res_s1, res_s2, res_p, res_q)
  local p, q, cste, s1, s2, d, L, i, tmp, exponent, r, polyring;
  p := formalRec:-numerator:
  q := formalRec:-denominator:
  d := formalRec:-order:
  s1 := 1:
  L := op(2,my_factors(p)): # L contains the non trivial factors of p
  for i from 1 to nops(L) do
    tmp := L[i][1]: exponent := L[i][2]:
    r := gcd(tmp^exponent, subs(n=n-d, q)):
    p := quo(p,r,n): q := quo(q, subs(n=n+d, r),n): s1 := s1 * r:
  od:

  s2 := 1:
  L := op(2,my_factors(p)): # L contains the *remaining* non trivial factors of p
  for i from 1 to nops(L) do
    tmp := L[i][1]: exponent := L[i][2]:
    r := gcd(tmp^exponent, subs(n=n+d, q)):
    p := quo(p, r, n): q := quo(q, subs(n=n-d, r), n): s2 := s2 * r:
  od:

  # Finally we look for the constant part (with respect to n) of p/q
  cste := op(1, my_factors(p))/op(1, my_factors(q)):
  p := p/op(1, my_factors(p)): q := q/op(1, my_factors(q)):
  polyring := RegularChains[PolynomialRing]([n,x]):
  L := op(2, my_factors(p)):
  for i from 1 to nops(L) do
    if RegularChains[MainVariable](L[i][1], polyring) = x
      then cste := cste * L[i][1]^L[i][2]: p := quo(p,L[i][1]^L[i][2],x):
    fi:
  od:
  L := op(2, my_factors(q)):
  for i from 1 to nops(L) do
    if RegularChains[MainVariable](L[i][1], polyring) = x
      then cste := cste / L[i][1]^L[i][2]: q := quo(q,L[i][1]^L[i][2],x):
    fi:
  od:

  res_cste := cste;
  res_s1 := s1;
  res_s2 := s2;
  res_p := simplify(p);
  res_q := simplify(q);
end proc:

#
#coeffrecToTermsrec
# Input: a linear recurrence rec (either with or without initial conditions).
#           a(n): the name of the sequence and the name of the variable.
#           x: a value or symbolic name
#
# Output:
# The recurrence satisfied by a(n)*x^n. Note that this recurrence is also denoted by a(n).
# If initial conditions were provided, corresponding initial conditions are computed.
#
coeffrecToTermsrec := proc(rec, aofn, x)
  local a,n,L,r,cond,d,i,tmp,c,res;
  getname(aofn, a, n):
  if type(rec, 'set') then
    L := selectremove(has, rec, n):
    r := L[1]:
    if nops(r)>1
      then error `invalid recurrence`, rec
    fi:
    if nops(r)=0
      then error "%1 does not appear in the recurrence", n
    fi:
    r := op(r):
    cond := L[2]:
  else r := rec:
  fi:
  d := maxindex(r, a, n):
  L := indets(r,'specfunc'('linear'(n),a)):
  if indets(r,'specfunc'('anything',a)) <> L
    then error "the recurrence contains elements that are not linear in %1", n
  fi:
  L := map(op, L):
  for i from 1 to nops(L) do
    r := subs(a(op(i,L))=a(op(i,L))*x^(d-op(i,L)+n), r):
  od:
  if cond<>'cond' then
    c := {}:
    for i from 1 to nops(cond) do
      tmp := op(i, cond): # tmp should have the form 'a(k) = cste'
      if not type(tmp,'`=`') then error "Invalid initial condition: %1", tmp: fi:
      L := selectremove(has, {op(tmp)}, a):
      if (nops(L[1]) <> 1) or (nops(L[2])<>1)
        then error "Invalid initial condition: %1", tmp:
      fi:
      tmp := op(1, L[1]): # tmp has the form 'a(k)'
      c := {op(c), tmp = op(1, L[2])*x^op(tmp)}:
    od:
    res := {r, op(c)}:
  else res := r:
  fi:
  res:
end proc:

#
# This procedure removes the conditions of the form a(k)=0 from the initial conditions of rec
# It returns a list L = [L1, L2, ...] where Li = [k, expr] representing the condition a(k)=expr.
# Moreover, it asserts that the Li are ordered by increasing k.
#
removeTrivialConditions := proc(rec, aofn)
  local a,n,i,L,tmp,c,cond,k:
  getname(aofn, a, n):
  if not type(rec, 'set') then
    error "%1 is not a recurrence with initial conditions", rec
  else
    L := selectremove(has, rec, n):
    cond := L[2]:
    if nops(cond)=0
      then error "%1 does not contain initial conditions", rec
    fi:
  fi:
  c := []:
  for i from 1 to nops(cond) do
    tmp := op(i, cond): # tmp should have the form 'a(k) = cste'
    if not type(tmp,'`=`') then error "Invalid initial condition: %1", tmp: fi:
    L := selectremove(has, {op(tmp)}, a):
    if (nops(L[1]) <> 1) or (nops(L[2])<>1)
      then error "Invalid initial condition: %1", tmp:
    fi:
    if op(1, L[2])<>0 then c := [op(c), [op(op(1, L[1])), op(1, L[2])]]: fi:
  od:
  # We check that the conditions are ordered by increasing k.
  if (nops(c)=0) then return c: fi:
  k := c[1][1]:
  for i from 2 to nops(c) do
    if (c[i][1]<=k)
    then error "Unexpected error in removeTrivialConditions: the conditions are not correctly ordered (%1)\n", c
    else k := c[i][1]
    fi:
  od:
  c:
end proc:

#
# findFixpointOfDifferences: takes a set L of integer and returns the smallest set S
# containing L and such that for each i, S[i]-S[i-1] \in S
findFixpointOfDifferences := proc(L)
  local res, i:
  res := L:
  for i from 2 to nops(L) do
    res := { op(res), L[i]-L[i-1] }:
  od:
  if (res=L) then return res else return findFixpointOfDifferences(res) fi:
end proc:

#
# error_counter functions allows one to follow the accumulation of errors in each variable.
#   an error_counter is a list of the form [[var1, c1], [var2, c2], ... ]
#   where the vari are variable names and the ci indicate how many approximation errors
#   are accumulated in vari.
#

#
# This procedure initializes the counter associated with variable var to 1 (and creates it if needed.)
# It returns an up-to-date error_counter.
init_error_counter := proc (var, error_counter)
  local i, res:
  res  := error_counter:
  for i from 1 to nops(res) do
    if (res[i][1]=var)
    then res[i][2] := 1:
         return res:
    fi
  od:
  res := [op(res), [var, 1]]:
end:

#
# This procedure adds a given number to the counter associated with variable var.
# It returns an up-to-date error_counter.
add_to_error_counter := proc (var, n, error_counter)
  local i, res:
  res  := error_counter:
  for i from 1 to nops(res) do
    if (res[i][1]=var)
    then res[i][2] := res[i][2]+n:
         return res:
    fi
  od:
  res := [op(res), [var, n]]:
end proc:

#
# This procedure sets the value of the counter associated with variable var.
# It returns an up-to-date error_counter.
set_error_counter := proc(var, n, error_counter)
  local i,err:
  err  := error_counter:
  for i from 1 to nops(err) do
    if (err[i][1]=var)
    then err[i][2] := n:
         return err:
    fi
  od:
  err := [op(err), [var, n]]:
end proc:

#
# This procedure initializes the counter associated to the multiplication of var2 and var3,
# putting the result in variable var1.
# It returns an up-to-date error_counter.
error_counter_of_a_multiplication := proc (var1, var2, var3, error_counter)
  local i, res, c2, c3:
  c2 := 0: c3 := 0:
  for i from 1 to nops(error_counter) do
    if (error_counter[i][1]=var2) then c2 := error_counter[i][2] fi:
    if (error_counter[i][1]=var3) then c3 := error_counter[i][2] fi:
    if (error_counter[i][1]=var1)
    then
      res := [ op(error_counter[1..i-1]), op(error_counter[i+1..nops(error_counter)]) ]
    fi:
  od:
  if (res = 'res') then res := error_counter fi:
  res := [op(res), [var1, c2+c3+1]]:
end:

#
# Copies the error counter of var2 into var1
error_counter_on_copy := proc(var1, var2, error_counter)
  local i, err, c2:
  c2 := 0:
  for i from 1 to nops(error_counter) do
    if (error_counter[i][1] = var2) then c2 := error_counter[i][2] fi:
    if (error_counter[i][1] = var1)
    then
      err := [ op(error_counter[1..i-1]), op(error_counter[i+1..nops(error_counter)]) ]
    fi:
  od:
  if (err = 'err') then err := error_counter fi:
  if (c2 <> 0) then err := [op(res), [var1, c2]] fi:
end proc:

#
# Returns the value of the error counter associated to a variable
find_in_error_counter := proc(var, error_counter)
  local i:
  for i from 1 to nops(error_counter) do
    if (error_counter[i][1] = var) then return error_counter[i][2] fi:
  od:
  return 0:
end proc:

#
# generate_multiply_rational(fd, var1, var2, r, error_counter, indent) generates code for performing
# var1 = var2*r in MPFR
#   fd is the file descriptor in which the code shall be produced.
#   var1 and var2 are strings representing variable names. r is a Maple rational number.
#   error_counter is an error_counter (as described above).
#   indent is an optional argument. It is a string used to correctly indent the code. It is prefixed to any
#   generated line. Hence, if indent="  ", the generated code will be indented by 2 spaces.
# An up-to-date error_counter is returned.
generate_multiply_rational := proc(fd, var1, var2, r, error_counter, indent:="")
  local p,q,err:
  err := error_counter:
  if (whattype(r)<>'fraction') and (whattype(r)<>'integer')
  then error "generate_multiply_rational used with non rational number %1", r: fi:
  if (abs(r)=1)
  then
    if (var1=var2)
    then
      if (r<>1) then fprintf(fd, "%sMPFR_CHANGE_SIGN (%s);\n", indent, var1) fi:
      return err:
    else
      if (r=1)
        then fprintf(fd, "%smpfr_set (%s, %s, MPFR_RNDN);\n", indent, var1, var2):
        else fprintf(fd, "%smpfr_neg (%s, %s, MPFR_RNDN);\n", indent, var1, var2):
      fi:
      return error_counter_on_copy(var1, var2, err):
    fi
  fi:
  # Now, r is a rational number different from 1.
  p := numer(r): q := denom(r):
  if (abs(p)<>1)
  then
    fprintf(fd, "%smpfr_mul_si (%s, %s, %d, MPFR_RNDN);\n", indent, var1, var2, p):
    err := error_counter_of_a_multiplication(var1, var2, "", err):
    if(q<>1)
    then
      fprintf(fd, "%smpfr_div_si (%s, %s, %d, MPFR_RNDN);\n", indent, var1, var1, q):
      err := error_counter_of_a_multiplication(var1, var1, "", err):
    fi:
  else
    fprintf(fd, "%smpfr_div_si (%s, %s, %d, MPFR_RNDN);\n", indent, var1, var2, p*q):
    err := error_counter_of_a_multiplication(var1, var2, "", err):
  fi:
  return err:
end proc:

#
# generate_multiply_poly is the same as generate_multiply_rational but when r is a rational fraction.
# The fraction r must have the form p/q where p and q are polynomials with integer coefficients.
# Moreover, the gcd of the coefficients of p must be 1. Idem for q.
# The procedure returned a list [m, d, err] where m is the set of indices k such that
# a mpfr_mul_sik function is needed and idem for d with mpfr_div_sik.
# err is an up-to-date error counter.
generate_multiply_poly := proc(fd, var1, var2, r, error_counter, indent:="")
  local p,q,Lp,Lq,n,i,j,var, required_mulsi, required_divsi, err:
  err := error_counter:
  required_mulsi := {}:
  required_divsi := {}:
  p := numer(r): q := denom(r):
  Lp := my_factors(p): Lq := my_factors(q):
  if (Lp[1] <> 1)
    then error "generate_multiply_poly: an integer can be factored out of %1", p:
  fi:
  if (Lq[1] <> 1)
    then error "generate_multiply_poly: an integer can be factored out of %1", q:
  fi:
  Lp := Lp[2]: Lq := Lq[2]:
  var := var2:
  if (nops(Lp) <> 0)
  then
    n := 0:
    for i from 1 to nops(Lp) do n := n + Lp[i][2] od:
    if (n=1)
    then
      fprintf(fd, "%smpfr_mul_si (%s, %s", indent, var1, var):
    else
      required_mulsi := { op(required_mulsi), n }:
      fprintf(fd, "%smpfr_mul_si%d (%s, %s", indent, n, var1, var):
    fi:
    for i from 1 to nops(Lp) do
      for j from 1 to Lp[i][2] do
        fprintf(fd, ", %a", Lp[i][1]):
      od:
    od:
    fprintf(fd, ", MPFR_RNDN);\n"):
    err := set_error_counter(var1, n+find_in_error_counter(var, err) , err):
    var := var1:
  fi:
  if (nops(Lq) <> 0)
  then
    n := 0:
    for i from 1 to nops(Lq) do n := n + Lq[i][2] od:
    if (n=1)
    then
      fprintf(fd, "%smpfr_div_si (%s, %s", indent, var1, var):
    else
      required_divsi := { op(required_divsi), n }:
      fprintf(fd, "%smpfr_div_si%d (%s, %s", indent, n, var1, var)
    fi:
    for i from 1 to nops(Lq) do
      for j from 1 to Lq[i][2] do
        fprintf(fd, ", %a", Lq[i][1])
      od:
    od:
    fprintf(fd, ", MPFR_RNDN);\n"):
    err := set_error_counter(var1, n+find_in_error_counter(var, err) , err):
    var := var1:
  fi:
  if (var1 <> var) then
    fprintf(fd, "%smpfr_set (%s, %s, MPFR_RNDN);\n", indent, var1, var):
    err := set_error_counter(var1, find_in_error_counter(var, err) , err):
  fi:
  return [required_mulsi, required_divsi, err]:
end proc:

#
# This function generates the code of a procedure mpfr_mul_uin or mpfr_div_uin
#
generate_muldivsin := proc(op, n)
  local i, var:
  if ((op <> "mul") and (op <> "div"))
    then error "Invalid argument to generate_muldivuin (%1). Must be 'mul' or 'div'", op
  fi:
  if (whattype(n) <> 'integer')
  then error "Invalid argument to generate_muldivuin (%1). Must be an integer.", n
  fi:

  if (op="mul") then var := "MUL" else var := "DIV" fi:

  printf("__MPFR_DECLSPEC void mpfr_div_si%d _MPFR_PROTO((mpfr_ptr, mpfr_srcptr,\n", n):
  for i from n to 2 by -2 do
    printf("                                               long int, long int,\n"):
  od:
  if (i=1)
  then
    printf("                                               long int, mpfr_rnd_t));\n"):
  else
    printf("                                               mpfr_rnd_t));\n")
  fi:

  printf("\n\n\n"):
  printf("void\n"):
  printf("mpfr_%s_si%d (mpfr_ptr y, mpfr_srcptr x,\n", op, n):
  for i from n to 2 by -2 do
    printf("              long int v%d, long int v%d,\n", n-i+1, n-i+2):
  od:
  if (i=1)
  then
    printf("              long int v%d, mpfr_rnd_t mode)\n", n):
  else
    printf("              mpfr_rnd_t mode)\n")
  fi:
  printf("{\n"):
  printf("  long int acc = v1;\n"):
  printf("  mpfr_set (y, x, mode);\n"):
   for i from 2 to n do
    printf("  MPFR_ACC_OR_%s (v%d);\n", var, i):
  od:
  printf("  mpfr_%s_si (y, y, acc, mode);\n", op):
  printf("}\n"):
  return:
end proc:


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

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

n12_omega_x := 2.5:
n12_gamma_x := 0.004:

n12_rss := (rs, z) -> rs * 2^(1/3) * opz_pow_n(z,-1/3):

n12_vx := rs -> 1/(1 + (1/(RS_FACTOR*n12_omega_x))*rs):
n12_ux := x -> n12_gamma_x*x^2/(1 + n12_gamma_x*x^2):

n12_FN12 := (rs, z, x) ->
  + add(params_a_CC_0_[i+1]*n12_ux(x)^i, i=0..3)
  + add(params_a_CC_1_[i+1]*n12_ux(x)^i, i=0..3) * n12_vx(n12_rss(rs, z))
  + add(params_a_CC_2_[i+1]*n12_ux(x)^i, i=0..3) * n12_vx(n12_rss(rs, z))^2
  + add(params_a_CC_3_[i+1]*n12_ux(x)^i, i=0..3) * n12_vx(n12_rss(rs, z))^3:

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


text:
use DocumentTools in

if parse(GetProperty("ComputeId", value)) then
	SetProperty("bypass", enabled, true):
	SetProperty("SkipSingle", enabled, true):
	SetProperty("SimplifiedGen", enabled, true):
	SetProperty("NoBound", enabled, true):
	SetProperty("Refine", enabled, true):

else
	SetProperty("bypass", enabled, false):
	SetProperty("SkipSingle", enabled, false):
	SetProperty("SimplifiedGen", enabled, false):
	SetProperty("NoBound", enabled, false):
	SetProperty("Refine", enabled, false):
	SetProperty("UsingUpTo", enabled, false):
	SetProperty("MaxPermutations", enabled, false):
	SetProperty("Permutations", enabled, false):

fi:
end use;

max_stars_repo_name: iliailmer/sian-web-app

max_issues_repo_name: iliailmer/sian-web-app

max_forks_repo_name: iliailmer/sian-web-app


text:
N := 3;
A := {1,2,3,4,5,6};
TT := `random_element/trees`(A)();

check_map_generic("phi/nonempty_subsets/prime_simplex","nonempty_subsets","prime_simplex",[[A]]);
check_map_generic("centres/cubes","cubes","F",[[N],[A]]);
check_map_generic("fatten/F","F","cubes",[[N],[A]]);
check_map_generic("inc/F/Fbar","F","Fbar",[[N],[A]]);
check_map_generic("extend/tree_height_functions","tree_height_functions","height_functions",
                  [[A],[TT]],[[A],[TT]],[[A]]);
check_map_generic("critical_tree/Fbar","Fbar","full_trees",[[N],[A]],[[N],[A]],[[A]]);

N := 2;
A := {1,2,3,4,5};

check_map_generic("phi/SEM/SCQ","SEM","SCQ",[[N],[A]]);
check_map_generic("psi/SCQ/SEM","SCQ","SEM",[[N],[A]]);
check_map_generic("sigma/SEM/F","SEM","F",[[N],[A]]);
check_map_generic("mu/F/SEM","F","SEM",[[N],[A]]);
check_map_generic("phi2/Fbar/P2","Fbar","P2",[[N],[A]]);
check_map_generic("psi/Fbar/Q","Fbar","Q",[[N],[A]],[[N],[A]],[[A]]);


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

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

(* type: lda_exc *)

$define xc_dimensions_2d

amgb_aa := [ -0.1925,     0.117331,    0.0234188 ]:
amgb_bb := [  0.0863136, -3.394e-2,   -0.037093  ]:
amgb_cc := [  0.0572384, -7.66765e-3,  0.0163618 ]:
amgb_ee := [  1.0022,     0.4133,      1.424301  ]:
amgb_ff := [ -0.02069,    0,         0       ]:
amgb_gg := [  0.33997,    6.68467e-2,  0       ]:
amgb_hh := [  1.747e-2,   7.799e-4,    1.163099  ]:

amgb_ax = -4/(3*Pi*sqrt(2)):
amgb_beta := 1.3386:

amgb_dd := [seq(-amgb_aa[i]*amgb_hh[i], i=1..3)]:

amgb_alpha := (i, rs) -> amgb_aa[i]
  + (amgb_bb[i]*rs + amgb_cc[i]*rs^2 + amgb_dd[i]*rs^3)
  * log(1 + 1/(amgb_ee[i]*rs + amgb_ff[i]*rs^1.5 + amgb_gg[i]*rs^2 + amgb_hh[i]*rs^3)):

amgb_ex6 := (rs, z) -> -4*sqrt(2)/(3*Pi*rs)
  * (f_zeta_2d(z) - 1 - 3/8*z^2 - 3/128*z^4):

f_amgb := (rs, z) ->
  amgb_alpha(1, rs) + amgb_alpha(2, rs)*z^2 + amgb_alpha(3, rs)*z^4 + (exp(-amgb_beta*rs) - 1)*amgb_ex6(rs, z):

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


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

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

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

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

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

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

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

`mu_aux0/BA` := proc(u,v)
 local r,s,i,ii,j;
 r := nops(u);
 s := nops(v);
 if r = 0 then return v; fi;
 if s = 0 then return u; fi;
 ii := [seq([i],i=1..r+s)];
 for i from 2 to r do
  ii := map(`ext/BA`,ii,r+s);
 od;
 return add(`mu_aux1/BA`(u,v,ii[j]),j=1..nops(ii));
end:

`ext/BA` := proc (i, r)
 local j;
 seq([op(i), j], j = i[nops(i)]+1 .. r);
end:

`mu_aux1/BA` := proc (u, v, i)
 local r,s,q,j,k,m;
 r := nops(u);
 s := nops(v);
 q := [EMPTY$r+s];
 for j to r do
  q[i[j]] := op(j, u)
 od;
 k := 1;
 for j to r+s do
  if q[j] = EMPTY then
   q[j] := op(k, v);
   k := k+1;
  fi;
 od;
 m := `+`(op(i))-(1/2)*r*(r+1);
 return (-1)^m*`bar/BA`(op(q))
end:

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

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

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

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

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

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

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

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

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

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

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

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

`e/EC` := proc(n::posint,k::nonnegint,i::nonnegint)
 '`e/EC`'(n,k,modp(i,n));
end:

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

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

`mu/EC` := proc(a,b)
 local xx,ca,cb,pa,pb,na,nb,ka,kb,la,lb,ia,ib;
 if nargs > 2 then
  return `mu/EC`(`mu/EC`(a,b),args[3..-1]);
 fi;
 if type(a,numeric) or type(b,numeric) then
  return expand(a*b);
 elif type(a,`+`) then
  return map(`mu/EC`,a,b);
 elif type(b,`+`) then
  return map2(`mu/EC`,a,b);
 fi;
 if type(a,`*`) then
  xx := selectremove(type,a,numeric);
  ca := xx[1];
  pa := xx[2];
 else
  ca := 1;
  pa := a;
 fi;
 if type(b,`*`) then
  xx := selectremove(type,b,numeric);
  cb := xx[1];
  pb := xx[2];
 else
  cb := 1;
  pb := b;
 fi;
 if type(pa,specfunc(anything,`e/EC`)) and
    type(pb,specfunc(anything,`e/EC`)) then
  na,ka,ia := op(pa);
  nb,kb,ib := op(pb);
  if nb <> na then return FAIL; fi;
  if modp(ka,2) = 1 and modp(kb,2) = 1 then
   return 0;
  fi;
  la := floor(ka/2);
  lb := floor(kb/2);
  return ca * cb * binomial(la+lb,la) * `e/EC`(na,ka+kb,modp(ia+ib,na));
 else
  return('`mu/EC`'(args))
 fi;
end:

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

`e/BC` := proc(n::posint,k::nonnegint)
 '`e/BC`'(n,k);
end:

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

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

`mu/BC` := proc(a,b)
 local xx,ca,cb,pa,pb,na,nb,ka,kb,la,lb,ia,ib;
 if nargs > 2 then
  return `mu/BC`(`mu/BA`(a,b),args[3..-1]);
 fi;
 if type(a,numeric) or type(b,numeric) then
  return expand(a*b);
 elif type(a,`+`) then
  return map(`mu/BC`,a,b);
 elif type(b,`+`) then
  return map2(`mu/BC`,a,b);
 fi;
 if type(a,`*`) then
  xx := selectremove(type,a,numeric);
  ca := xx[1];
  pa := xx[2];
 else
  ca := 1;
  pa := a;
 fi;
 if type(b,`*`) then
  xx := selectremove(type,b,numeric);
  cb := xx[1];
  pb := xx[2];
 else
  cb := 1;
  pb := b;
 fi;
 if type(pa,specfunc(anything,`e/BC`)) and
    type(pb,specfunc(anything,`e/BC`)) then
  na,ka := op(pa);
  nb,kb := op(pb);
  if nb <> na then return FAIL; fi;
  if modp(ka,2) = 1 and modp(kb,2) = 1 then
   return 0;
  fi;
  la := floor(ka/2);
  lb := floor(kb/2);
  return ca * cb * binomial(la+lb,la) * `e/BC`(na,ka+kb);
 else
  return('`mu/BC`'(args))
 fi;
end:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

 n := 5;

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

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

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

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

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

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


text:
######################################################################
# pair_subsets(A) is the set of all subsets B of A with |B| = 2

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

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

`is_leq/pair_subsets` := NULL;

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

 if nops(A) < 2 then return FAIL; fi;

 r := rand(1..nops(A));
 B := {};

 while nops(B) < 2 do
  B := [A[r()],A[r()]];
  B := {op(B)};
 od;

 return B;
end;

`list_elements/pair_subsets` := proc(A::set)
 local i,j;
 return [seq(seq({A[i],A[j]},j=i+1..nops(A)),i=1..nops(A)-1)];
end:

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


text:
Q := [b*c+a, a*c-b, b*c-2*c^2+a-2, b*c-c^2+a-1, b*c+c^2+a+1, b*c+2*c^2+a+2, 2*
b*c-c^2+2*a-1, 2*b*c+c^2+2*a+1, a*c-2*c^2-b-2, a*c-c^2-b-1, a*c+c^2-b+1, a*c+2
*c^2-b+2, 2*a*c-c^2-2*b-1, 2*a*c+c^2-2*b+1, a*c-2*b*c-2*a-b, a*c-b*c-a-b, a*c+
b*c+a-b, a*c+2*b*c+2*a-b, 2*a*c-b*c-a-2*b, 2*a*c+b*c+a-2*b, a^2+b^2-3*c^2-3, a
^2+b^2-c^2-1, 3*a^2+3*b^2-c^2-1, a^2+b^2-2*b*c-2*a, a^2+b^2-b*c-a, a^2+b^2+b*c
+a, a^2+b^2+2*b*c+2*a, 2*a^2+2*b^2-b*c-a, 2*a^2+2*b^2+b*c+a, a^2-2*a*c+b^2+2*b
, a^2-a*c+b^2+b, a^2+a*c+b^2-b, a^2+2*a*c+b^2-2*b, 2*a^2-a*c+2*b^2+b, 2*a^2+a*
c+2*b^2-b, a*c-2*b*c-c^2-2*a-b-1, a*c-2*b*c+c^2-2*a-b+1, a*c-b*c-2*c^2-a-b-2,
a*c-b*c-c^2-a-b-1, a*c-b*c+c^2-a-b+1, a*c-b*c+2*c^2-a-b+2, a*c+b*c-2*c^2+a-b-2
, a*c+b*c-c^2+a-b-1, a*c+b*c+c^2+a-b+1, a*c+b*c+2*c^2+a-b+2, a*c+2*b*c-c^2+2*a
-b-1, a*c+2*b*c+c^2+2*a-b+1, 2*a*c-b*c-c^2-a-2*b-1, 2*a*c-b*c+c^2-a-2*b+1, 2*a
*c+b*c-c^2+a-2*b-1, 2*a*c+b*c+c^2+a-2*b+1, a^2+b^2-4*b*c-3*c^2-4*a-3, a^2+b^2-\
4*b*c-c^2-4*a-1, a^2+b^2-4*b*c+c^2-4*a+1, a^2+b^2-4*b*c+3*c^2-4*a+3, a^2+b^2-2
*b*c-3*c^2-2*a-3, a^2+b^2-2*b*c-c^2-2*a-1, a^2+b^2-b*c-c^2-a-1, a^2+b^2+b*c-c^
2+a-1, a^2+b^2+2*b*c-3*c^2+2*a-3, a^2+b^2+2*b*c-c^2+2*a-1, a^2+b^2+4*b*c-3*c^2
+4*a-3, a^2+b^2+4*b*c-c^2+4*a-1, a^2+b^2+4*b*c+c^2+4*a+1, a^2+b^2+4*b*c+3*c^2+
4*a+3, 3*a^2+3*b^2-4*b*c-c^2-4*a-1, 3*a^2+3*b^2-4*b*c+c^2-4*a+1, 3*a^2+3*b^2-2
*b*c-c^2-2*a-1, 3*a^2+3*b^2+2*b*c-c^2+2*a-1, 3*a^2+3*b^2+4*b*c-c^2+4*a-1, 3*a^
2+3*b^2+4*b*c+c^2+4*a+1, a^2-4*a*c+b^2-3*c^2+4*b-3, a^2-4*a*c+b^2-c^2+4*b-1, a
^2-4*a*c+b^2+c^2+4*b+1, a^2-4*a*c+b^2+3*c^2+4*b+3, a^2-2*a*c+b^2-3*c^2+2*b-3,
a^2-2*a*c+b^2-c^2+2*b-1, a^2-a*c+b^2-c^2+b-1, a^2+a*c+b^2-c^2-b-1, a^2+2*a*c+b
^2-3*c^2-2*b-3, a^2+2*a*c+b^2-c^2-2*b-1, a^2+4*a*c+b^2-3*c^2-4*b-3, a^2+4*a*c+
b^2-c^2-4*b-1, a^2+4*a*c+b^2+c^2-4*b+1, a^2+4*a*c+b^2+3*c^2-4*b+3, 3*a^2-4*a*c
+3*b^2-c^2+4*b-1, 3*a^2-4*a*c+3*b^2+c^2+4*b+1, 3*a^2-2*a*c+3*b^2-c^2+2*b-1, 3*
a^2+2*a*c+3*b^2-c^2-2*b-1, 3*a^2+4*a*c+3*b^2-c^2-4*b-1, 3*a^2+4*a*c+3*b^2+c^2-\
4*b+1, a^2-2*a*c+b^2-b*c-a+2*b, a^2-2*a*c+b^2+b*c+a+2*b, a^2-a*c+b^2-2*b*c-2*a
+b, a^2-a*c+b^2-b*c-a+b, a^2-a*c+b^2+b*c+a+b, a^2-a*c+b^2+2*b*c+2*a+b, a^2+a*c
+b^2-2*b*c-2*a-b, a^2+a*c+b^2-b*c-a-b, a^2+a*c+b^2+b*c+a-b, a^2+a*c+b^2+2*b*c+
2*a-b, a^2+2*a*c+b^2-b*c-a-2*b, a^2+2*a*c+b^2+b*c+a-2*b, 2*a^2-a*c+2*b^2-b*c-a
+b, 2*a^2-a*c+2*b^2+b*c+a+b, 2*a^2+a*c+2*b^2-b*c-a-b, 2*a^2+a*c+2*b^2+b*c+a-b,
a^2-4*a*c+b^2-4*b*c+c^2-4*a+4*b+1, a^2-4*a*c+b^2-2*b*c-c^2-2*a+4*b-1, a^2-4*a*
c+b^2-2*b*c+c^2-2*a+4*b+1, a^2-4*a*c+b^2-2*b*c+3*c^2-2*a+4*b+3, a^2-4*a*c+b^2+
2*b*c-c^2+2*a+4*b-1, a^2-4*a*c+b^2+2*b*c+c^2+2*a+4*b+1, a^2-4*a*c+b^2+2*b*c+3*
c^2+2*a+4*b+3, a^2-4*a*c+b^2+4*b*c+c^2+4*a+4*b+1, a^2-2*a*c+b^2-4*b*c-c^2-4*a+
2*b-1, a^2-2*a*c+b^2-4*b*c+c^2-4*a+2*b+1, a^2-2*a*c+b^2-4*b*c+3*c^2-4*a+2*b+3,
a^2-2*a*c+b^2-2*b*c-3*c^2-2*a+2*b-3, a^2-2*a*c+b^2-2*b*c-c^2-2*a+2*b-1, a^2-2*
a*c+b^2-2*b*c+c^2-2*a+2*b+1, a^2-2*a*c+b^2-b*c+c^2-a+2*b+1, a^2-2*a*c+b^2+b*c+
c^2+a+2*b+1, a^2-2*a*c+b^2+2*b*c-3*c^2+2*a+2*b-3, a^2-2*a*c+b^2+2*b*c-c^2+2*a+
2*b-1, a^2-2*a*c+b^2+2*b*c+c^2+2*a+2*b+1, a^2-2*a*c+b^2+4*b*c-c^2+4*a+2*b-1, a
^2-2*a*c+b^2+4*b*c+c^2+4*a+2*b+1, a^2-2*a*c+b^2+4*b*c+3*c^2+4*a+2*b+3, a^2-a*c
+b^2-2*b*c+c^2-2*a+b+1, a^2-a*c+b^2-b*c-c^2-a+b-1, a^2-a*c+b^2+b*c-c^2+a+b-1,
a^2-a*c+b^2+2*b*c+c^2+2*a+b+1, a^2+a*c+b^2-2*b*c+c^2-2*a-b+1, a^2+a*c+b^2-b*c-
c^2-a-b-1, a^2+a*c+b^2+b*c-c^2+a-b-1, a^2+a*c+b^2+2*b*c+c^2+2*a-b+1, a^2+2*a*c
+b^2-4*b*c-c^2-4*a-2*b-1, a^2+2*a*c+b^2-4*b*c+c^2-4*a-2*b+1, a^2+2*a*c+b^2-4*b
*c+3*c^2-4*a-2*b+3, a^2+2*a*c+b^2-2*b*c-3*c^2-2*a-2*b-3, a^2+2*a*c+b^2-2*b*c-c
^2-2*a-2*b-1, a^2+2*a*c+b^2-2*b*c+c^2-2*a-2*b+1, a^2+2*a*c+b^2-b*c+c^2-a-2*b+1
, a^2+2*a*c+b^2+b*c+c^2+a-2*b+1, a^2+2*a*c+b^2+2*b*c-3*c^2+2*a-2*b-3, a^2+2*a*
c+b^2+2*b*c-c^2+2*a-2*b-1, a^2+2*a*c+b^2+2*b*c+c^2+2*a-2*b+1, a^2+2*a*c+b^2+4*
b*c-c^2+4*a-2*b-1, a^2+2*a*c+b^2+4*b*c+c^2+4*a-2*b+1, a^2+2*a*c+b^2+4*b*c+3*c^
2+4*a-2*b+3, a^2+4*a*c+b^2-4*b*c+c^2-4*a-4*b+1, a^2+4*a*c+b^2-2*b*c-c^2-2*a-4*
b-1, a^2+4*a*c+b^2-2*b*c+c^2-2*a-4*b+1, a^2+4*a*c+b^2-2*b*c+3*c^2-2*a-4*b+3, a
^2+4*a*c+b^2+2*b*c-c^2+2*a-4*b-1, a^2+4*a*c+b^2+2*b*c+c^2+2*a-4*b+1, a^2+4*a*c
+b^2+2*b*c+3*c^2+2*a-4*b+3, a^2+4*a*c+b^2+4*b*c+c^2+4*a-4*b+1, 3*a^2-4*a*c+3*b
^2-2*b*c+c^2-2*a+4*b+1, 3*a^2-4*a*c+3*b^2+2*b*c+c^2+2*a+4*b+1, 3*a^2-2*a*c+3*b
^2-4*b*c+c^2-4*a+2*b+1, 3*a^2-2*a*c+3*b^2-2*b*c-c^2-2*a+2*b-1, 3*a^2-2*a*c+3*b
^2+2*b*c-c^2+2*a+2*b-1, 3*a^2-2*a*c+3*b^2+4*b*c+c^2+4*a+2*b+1, 3*a^2+2*a*c+3*b
^2-4*b*c+c^2-4*a-2*b+1, 3*a^2+2*a*c+3*b^2-2*b*c-c^2-2*a-2*b-1, 3*a^2+2*a*c+3*b
^2+2*b*c-c^2+2*a-2*b-1, 3*a^2+2*a*c+3*b^2+4*b*c+c^2+4*a-2*b+1, 3*a^2+4*a*c+3*b
^2-2*b*c+c^2-2*a-4*b+1, 3*a^2+4*a*c+3*b^2+2*b*c+c^2+2*a-4*b+1, -b*c+a, a*b+c,
-2*b^2-b*c+a-2, -b^2-2*b*c+2*a-1, -b^2-b*c+a-1, b^2-2*b*c+2*a+1, b^2-b*c+a+1,
2*b^2-b*c+a+2, a*b-2*b*c+2*a+c, a*b-b*c+a+c, a*b+b*c-a+c, a*b+2*b*c-2*a+c, 2*a
*b-b*c+a+2*c, 2*a*b+b*c-a+2*c, a*b-2*b^2+c-2, a*b-b^2+c-1, a*b+b^2+c+1, a*b+2*
b^2+c+2, 2*a*b-b^2+2*c-1, 2*a*b+b^2+2*c+1, a^2-2*b*c+c^2+2*a, a^2-b*c+c^2+a, a
^2+b*c+c^2-a, a^2+2*b*c+c^2-2*a, 2*a^2-b*c+2*c^2+a, 2*a^2+b*c+2*c^2-a, a^2-3*b
^2+c^2-3, a^2-b^2+c^2-1, 3*a^2-b^2+3*c^2-1, a^2-2*a*b+c^2-2*c, a^2-a*b+c^2-c,
a^2+a*b+c^2+c, a^2+2*a*b+c^2+2*c, 2*a^2-a*b+2*c^2-c, 2*a^2+a*b+2*c^2+c, a*b-2*
b^2-b*c+a+c-2, a*b-2*b^2+b*c-a+c-2, a*b-b^2-2*b*c+2*a+c-1, a*b-b^2-b*c+a+c-1,
a*b-b^2+b*c-a+c-1, a*b-b^2+2*b*c-2*a+c-1, a*b+b^2-2*b*c+2*a+c+1, a*b+b^2-b*c+a
+c+1, a*b+b^2+b*c-a+c+1, a*b+b^2+2*b*c-2*a+c+1, a*b+2*b^2-b*c+a+c+2, a*b+2*b^2
+b*c-a+c+2, 2*a*b-b^2-b*c+a+2*c-1, 2*a*b-b^2+b*c-a+2*c-1, 2*a*b+b^2-b*c+a+2*c+
1, 2*a*b+b^2+b*c-a+2*c+1, a^2-3*b^2-4*b*c+c^2+4*a-3, a^2-3*b^2-2*b*c+c^2+2*a-3
, a^2-3*b^2+2*b*c+c^2-2*a-3, a^2-3*b^2+4*b*c+c^2-4*a-3, a^2-b^2-4*b*c+c^2+4*a-\
1, a^2-b^2-2*b*c+c^2+2*a-1, a^2-b^2-b*c+c^2+a-1, a^2-b^2+b*c+c^2-a-1, a^2-b^2+
2*b*c+c^2-2*a-1, a^2-b^2+4*b*c+c^2-4*a-1, a^2+b^2-4*b*c+c^2+4*a+1, a^2+b^2+4*b
*c+c^2-4*a+1, a^2+3*b^2-4*b*c+c^2+4*a+3, a^2+3*b^2+4*b*c+c^2-4*a+3, 3*a^2-b^2-\
4*b*c+3*c^2+4*a-1, 3*a^2-b^2-2*b*c+3*c^2+2*a-1, 3*a^2-b^2+2*b*c+3*c^2-2*a-1, 3
*a^2-b^2+4*b*c+3*c^2-4*a-1, 3*a^2+b^2-4*b*c+3*c^2+4*a+1, 3*a^2+b^2+4*b*c+3*c^2
-4*a+1, a^2-2*a*b-b*c+c^2+a-2*c, a^2-2*a*b+b*c+c^2-a-2*c, a^2-a*b-2*b*c+c^2+2*
a-c, a^2-a*b-b*c+c^2+a-c, a^2-a*b+b*c+c^2-a-c, a^2-a*b+2*b*c+c^2-2*a-c, a^2+a*
b-2*b*c+c^2+2*a+c, a^2+a*b-b*c+c^2+a+c, a^2+a*b+b*c+c^2-a+c, a^2+a*b+2*b*c+c^2
-2*a+c, a^2+2*a*b-b*c+c^2+a+2*c, a^2+2*a*b+b*c+c^2-a+2*c, 2*a^2-a*b-b*c+2*c^2+
a-c, 2*a^2-a*b+b*c+2*c^2-a-c, 2*a^2+a*b-b*c+2*c^2+a+c, 2*a^2+a*b+b*c+2*c^2-a+c
, a^2-4*a*b-3*b^2+c^2-4*c-3, a^2-4*a*b-b^2+c^2-4*c-1, a^2-4*a*b+b^2+c^2-4*c+1,
a^2-4*a*b+3*b^2+c^2-4*c+3, a^2-2*a*b-3*b^2+c^2-2*c-3, a^2-2*a*b-b^2+c^2-2*c-1,
a^2-a*b-b^2+c^2-c-1, a^2+a*b-b^2+c^2+c-1, a^2+2*a*b-3*b^2+c^2+2*c-3, a^2+2*a*b
-b^2+c^2+2*c-1, a^2+4*a*b-3*b^2+c^2+4*c-3, a^2+4*a*b-b^2+c^2+4*c-1, a^2+4*a*b+
b^2+c^2+4*c+1, a^2+4*a*b+3*b^2+c^2+4*c+3, 3*a^2-4*a*b-b^2+3*c^2-4*c-1, 3*a^2-4
*a*b+b^2+3*c^2-4*c+1, 3*a^2-2*a*b-b^2+3*c^2-2*c-1, 3*a^2+2*a*b-b^2+3*c^2+2*c-1
, 3*a^2+4*a*b-b^2+3*c^2+4*c-1, 3*a^2+4*a*b+b^2+3*c^2+4*c+1, a^2-4*a*b-b^2-2*b*
c+c^2+2*a-4*c-1, a^2-4*a*b-b^2+2*b*c+c^2-2*a-4*c-1, a^2-4*a*b+b^2-4*b*c+c^2+4*
a-4*c+1, a^2-4*a*b+b^2-2*b*c+c^2+2*a-4*c+1, a^2-4*a*b+b^2+2*b*c+c^2-2*a-4*c+1,
a^2-4*a*b+b^2+4*b*c+c^2-4*a-4*c+1, a^2-4*a*b+3*b^2-2*b*c+c^2+2*a-4*c+3, a^2-4*
a*b+3*b^2+2*b*c+c^2-2*a-4*c+3, a^2-2*a*b-3*b^2-2*b*c+c^2+2*a-2*c-3, a^2-2*a*b-\
3*b^2+2*b*c+c^2-2*a-2*c-3, a^2-2*a*b-b^2-4*b*c+c^2+4*a-2*c-1, a^2-2*a*b-b^2-2*
b*c+c^2+2*a-2*c-1, a^2-2*a*b-b^2+2*b*c+c^2-2*a-2*c-1, a^2-2*a*b-b^2+4*b*c+c^2-\
4*a-2*c-1, a^2-2*a*b+b^2-4*b*c+c^2+4*a-2*c+1, a^2-2*a*b+b^2-2*b*c+c^2+2*a-2*c+
1, a^2-2*a*b+b^2-b*c+c^2+a-2*c+1, a^2-2*a*b+b^2+b*c+c^2-a-2*c+1, a^2-2*a*b+b^2
+2*b*c+c^2-2*a-2*c+1, a^2-2*a*b+b^2+4*b*c+c^2-4*a-2*c+1, a^2-2*a*b+3*b^2-4*b*c
+c^2+4*a-2*c+3, a^2-2*a*b+3*b^2+4*b*c+c^2-4*a-2*c+3, a^2-a*b-b^2-b*c+c^2+a-c-1
, a^2-a*b-b^2+b*c+c^2-a-c-1, a^2-a*b+b^2-2*b*c+c^2+2*a-c+1, a^2-a*b+b^2+2*b*c+
c^2-2*a-c+1, a^2+a*b-b^2-b*c+c^2+a+c-1, a^2+a*b-b^2+b*c+c^2-a+c-1, a^2+a*b+b^2
-2*b*c+c^2+2*a+c+1, a^2+a*b+b^2+2*b*c+c^2-2*a+c+1, a^2+2*a*b-3*b^2-2*b*c+c^2+2
*a+2*c-3, a^2+2*a*b-3*b^2+2*b*c+c^2-2*a+2*c-3, a^2+2*a*b-b^2-4*b*c+c^2+4*a+2*c
-1, a^2+2*a*b-b^2-2*b*c+c^2+2*a+2*c-1, a^2+2*a*b-b^2+2*b*c+c^2-2*a+2*c-1, a^2+
2*a*b-b^2+4*b*c+c^2-4*a+2*c-1, a^2+2*a*b+b^2-4*b*c+c^2+4*a+2*c+1, a^2+2*a*b+b^
2-2*b*c+c^2+2*a+2*c+1, a^2+2*a*b+b^2-b*c+c^2+a+2*c+1, a^2+2*a*b+b^2+b*c+c^2-a+
2*c+1, a^2+2*a*b+b^2+2*b*c+c^2-2*a+2*c+1, a^2+2*a*b+b^2+4*b*c+c^2-4*a+2*c+1, a
^2+2*a*b+3*b^2-4*b*c+c^2+4*a+2*c+3, a^2+2*a*b+3*b^2+4*b*c+c^2-4*a+2*c+3, a^2+4
*a*b-b^2-2*b*c+c^2+2*a+4*c-1, a^2+4*a*b-b^2+2*b*c+c^2-2*a+4*c-1, a^2+4*a*b+b^2
-4*b*c+c^2+4*a+4*c+1, a^2+4*a*b+b^2-2*b*c+c^2+2*a+4*c+1, a^2+4*a*b+b^2+2*b*c+c
^2-2*a+4*c+1, a^2+4*a*b+b^2+4*b*c+c^2-4*a+4*c+1, a^2+4*a*b+3*b^2-2*b*c+c^2+2*a
+4*c+3, a^2+4*a*b+3*b^2+2*b*c+c^2-2*a+4*c+3, 3*a^2-4*a*b+b^2-2*b*c+3*c^2+2*a-4
*c+1, 3*a^2-4*a*b+b^2+2*b*c+3*c^2-2*a-4*c+1, 3*a^2-2*a*b-b^2-2*b*c+3*c^2+2*a-2
*c-1, 3*a^2-2*a*b-b^2+2*b*c+3*c^2-2*a-2*c-1, 3*a^2-2*a*b+b^2-4*b*c+3*c^2+4*a-2
*c+1, 3*a^2-2*a*b+b^2+4*b*c+3*c^2-4*a-2*c+1, 3*a^2+2*a*b-b^2-2*b*c+3*c^2+2*a+2
*c-1, 3*a^2+2*a*b-b^2+2*b*c+3*c^2-2*a+2*c-1, 3*a^2+2*a*b+b^2-4*b*c+3*c^2+4*a+2
*c+1, 3*a^2+2*a*b+b^2+4*b*c+3*c^2-4*a+2*c+1, 3*a^2+4*a*b+b^2-2*b*c+3*c^2+2*a+4
*c+1, 3*a^2+4*a*b+b^2+2*b*c+3*c^2-2*a+4*c+1, a*c+b, a*b-c, a*c-2*b^2-2*c^2+b,
a*c-b^2-c^2+b, a*c+b^2+c^2+b, a*c+2*b^2+2*c^2+b, 2*a*c-b^2-c^2+2*b, 2*a*c+b^2+
c^2+2*b, a*b-2*b^2-2*c^2-c, a*b-b^2-c^2-c, a*b+b^2+c^2-c, a*b+2*b^2+2*c^2-c, 2
*a*b-b^2-c^2-2*c, 2*a*b+b^2+c^2-2*c, a*b-2*a*c-2*b-c, a*b-a*c-b-c, a*b+a*c+b-c
, a*b+2*a*c+2*b-c, 2*a*b-a*c-b-2*c, 2*a*b+a*c+b-2*c, a^2-3*b^2-3*c^2+1, a^2-b^
2-c^2+1, 3*a^2-b^2-c^2+3, a^2-2*a*c-2*b+1, a^2-a*c-b+1, a^2+a*c+b+1, a^2+2*a*c
+2*b+1, 2*a^2-a*c-b+2, 2*a^2+a*c+b+2, a^2-2*a*b+2*c+1, a^2-a*b+c+1, a^2+a*b-c+
1, a^2+2*a*b-2*c+1, 2*a^2-a*b+c+2, 2*a^2+a*b-c+2, a*b-2*a*c-b^2-c^2-2*b-c, a*b
-2*a*c+b^2+c^2-2*b-c, a*b-a*c-2*b^2-2*c^2-b-c, a*b-a*c-b^2-c^2-b-c, a*b-a*c+b^
2+c^2-b-c, a*b-a*c+2*b^2+2*c^2-b-c, a*b+a*c-2*b^2-2*c^2+b-c, a*b+a*c-b^2-c^2+b
-c, a*b+a*c+b^2+c^2+b-c, a*b+a*c+2*b^2+2*c^2+b-c, a*b+2*a*c-b^2-c^2+2*b-c, a*b
+2*a*c+b^2+c^2+2*b-c, 2*a*b-a*c-b^2-c^2-b-2*c, 2*a*b-a*c+b^2+c^2-b-2*c, 2*a*b+
a*c-b^2-c^2+b-2*c, 2*a*b+a*c+b^2+c^2+b-2*c, a^2-4*a*c-3*b^2-3*c^2-4*b+1, a^2-4
*a*c-b^2-c^2-4*b+1, a^2-4*a*c+b^2+c^2-4*b+1, a^2-4*a*c+3*b^2+3*c^2-4*b+1, a^2-\
2*a*c-3*b^2-3*c^2-2*b+1, a^2-2*a*c-b^2-c^2-2*b+1, a^2-a*c-b^2-c^2-b+1, a^2+a*c
-b^2-c^2+b+1, a^2+2*a*c-3*b^2-3*c^2+2*b+1, a^2+2*a*c-b^2-c^2+2*b+1, a^2+4*a*c-\
3*b^2-3*c^2+4*b+1, a^2+4*a*c-b^2-c^2+4*b+1, a^2+4*a*c+b^2+c^2+4*b+1, a^2+4*a*c
+3*b^2+3*c^2+4*b+1, 3*a^2-4*a*c-b^2-c^2-4*b+3, 3*a^2-4*a*c+b^2+c^2-4*b+3, 3*a^
2-2*a*c-b^2-c^2-2*b+3, 3*a^2+2*a*c-b^2-c^2+2*b+3, 3*a^2+4*a*c-b^2-c^2+4*b+3, 3
*a^2+4*a*c+b^2+c^2+4*b+3, a^2-4*a*b-3*b^2-3*c^2+4*c+1, a^2-4*a*b-b^2-c^2+4*c+1
, a^2-4*a*b+b^2+c^2+4*c+1, a^2-4*a*b+3*b^2+3*c^2+4*c+1, a^2-2*a*b-3*b^2-3*c^2+
2*c+1, a^2-2*a*b-b^2-c^2+2*c+1, a^2-a*b-b^2-c^2+c+1, a^2+a*b-b^2-c^2-c+1, a^2+
2*a*b-3*b^2-3*c^2-2*c+1, a^2+2*a*b-b^2-c^2-2*c+1, a^2+4*a*b-3*b^2-3*c^2-4*c+1,
a^2+4*a*b-b^2-c^2-4*c+1, a^2+4*a*b+b^2+c^2-4*c+1, a^2+4*a*b+3*b^2+3*c^2-4*c+1,
3*a^2-4*a*b-b^2-c^2+4*c+3, 3*a^2-4*a*b+b^2+c^2+4*c+3, 3*a^2-2*a*b-b^2-c^2+2*c+
3, 3*a^2+2*a*b-b^2-c^2-2*c+3, 3*a^2+4*a*b-b^2-c^2-4*c+3, 3*a^2+4*a*b+b^2+c^2-4
*c+3, a^2-2*a*b-a*c-b+2*c+1, a^2-2*a*b+a*c+b+2*c+1, a^2-a*b-2*a*c-2*b+c+1, a^2
-a*b-a*c-b+c+1, a^2-a*b+a*c+b+c+1, a^2-a*b+2*a*c+2*b+c+1, a^2+a*b-2*a*c-2*b-c+
1, a^2+a*b-a*c-b-c+1, a^2+a*b+a*c+b-c+1, a^2+a*b+2*a*c+2*b-c+1, a^2+2*a*b-a*c-
b-2*c+1, a^2+2*a*b+a*c+b-2*c+1, 2*a^2-a*b-a*c-b+c+2, 2*a^2-a*b+a*c+b+c+2, 2*a^
2+a*b-a*c-b-c+2, 2*a^2+a*b+a*c+b-c+2, a^2-4*a*b-4*a*c+b^2+c^2-4*b+4*c+1, a^2-4
*a*b-2*a*c-b^2-c^2-2*b+4*c+1, a^2-4*a*b-2*a*c+b^2+c^2-2*b+4*c+1, a^2-4*a*b-2*a
*c+3*b^2+3*c^2-2*b+4*c+1, a^2-4*a*b+2*a*c-b^2-c^2+2*b+4*c+1, a^2-4*a*b+2*a*c+b
^2+c^2+2*b+4*c+1, a^2-4*a*b+2*a*c+3*b^2+3*c^2+2*b+4*c+1, a^2-4*a*b+4*a*c+b^2+c
^2+4*b+4*c+1, a^2-2*a*b-4*a*c-b^2-c^2-4*b+2*c+1, a^2-2*a*b-4*a*c+b^2+c^2-4*b+2
*c+1, a^2-2*a*b-4*a*c+3*b^2+3*c^2-4*b+2*c+1, a^2-2*a*b-2*a*c-3*b^2-3*c^2-2*b+2
*c+1, a^2-2*a*b-2*a*c-b^2-c^2-2*b+2*c+1, a^2-2*a*b-2*a*c+b^2+c^2-2*b+2*c+1, a^
2-2*a*b-a*c+b^2+c^2-b+2*c+1, a^2-2*a*b+a*c+b^2+c^2+b+2*c+1, a^2-2*a*b+2*a*c-3*
b^2-3*c^2+2*b+2*c+1, a^2-2*a*b+2*a*c-b^2-c^2+2*b+2*c+1, a^2-2*a*b+2*a*c+b^2+c^
2+2*b+2*c+1, a^2-2*a*b+4*a*c-b^2-c^2+4*b+2*c+1, a^2-2*a*b+4*a*c+b^2+c^2+4*b+2*
c+1, a^2-2*a*b+4*a*c+3*b^2+3*c^2+4*b+2*c+1, a^2-a*b-2*a*c+b^2+c^2-2*b+c+1, a^2
-a*b-a*c-b^2-c^2-b+c+1, a^2-a*b+a*c-b^2-c^2+b+c+1, a^2-a*b+2*a*c+b^2+c^2+2*b+c
+1, a^2+a*b-2*a*c+b^2+c^2-2*b-c+1, a^2+a*b-a*c-b^2-c^2-b-c+1, a^2+a*b+a*c-b^2-
c^2+b-c+1, a^2+a*b+2*a*c+b^2+c^2+2*b-c+1, a^2+2*a*b-4*a*c-b^2-c^2-4*b-2*c+1, a
^2+2*a*b-4*a*c+b^2+c^2-4*b-2*c+1, a^2+2*a*b-4*a*c+3*b^2+3*c^2-4*b-2*c+1, a^2+2
*a*b-2*a*c-3*b^2-3*c^2-2*b-2*c+1, a^2+2*a*b-2*a*c-b^2-c^2-2*b-2*c+1, a^2+2*a*b
-2*a*c+b^2+c^2-2*b-2*c+1, a^2+2*a*b-a*c+b^2+c^2-b-2*c+1, a^2+2*a*b+a*c+b^2+c^2
+b-2*c+1, a^2+2*a*b+2*a*c-3*b^2-3*c^2+2*b-2*c+1, a^2+2*a*b+2*a*c-b^2-c^2+2*b-2
*c+1, a^2+2*a*b+2*a*c+b^2+c^2+2*b-2*c+1, a^2+2*a*b+4*a*c-b^2-c^2+4*b-2*c+1, a^
2+2*a*b+4*a*c+b^2+c^2+4*b-2*c+1, a^2+2*a*b+4*a*c+3*b^2+3*c^2+4*b-2*c+1, a^2+4*
a*b-4*a*c+b^2+c^2-4*b-4*c+1, a^2+4*a*b-2*a*c-b^2-c^2-2*b-4*c+1, a^2+4*a*b-2*a*
c+b^2+c^2-2*b-4*c+1, a^2+4*a*b-2*a*c+3*b^2+3*c^2-2*b-4*c+1, a^2+4*a*b+2*a*c-b^
2-c^2+2*b-4*c+1, a^2+4*a*b+2*a*c+b^2+c^2+2*b-4*c+1, a^2+4*a*b+2*a*c+3*b^2+3*c^
2+2*b-4*c+1, a^2+4*a*b+4*a*c+b^2+c^2+4*b-4*c+1, 3*a^2-4*a*b-2*a*c+b^2+c^2-2*b+
4*c+3, 3*a^2-4*a*b+2*a*c+b^2+c^2+2*b+4*c+3, 3*a^2-2*a*b-4*a*c+b^2+c^2-4*b+2*c+
3, 3*a^2-2*a*b-2*a*c-b^2-c^2-2*b+2*c+3, 3*a^2-2*a*b+2*a*c-b^2-c^2+2*b+2*c+3, 3
*a^2-2*a*b+4*a*c+b^2+c^2+4*b+2*c+3, 3*a^2+2*a*b-4*a*c+b^2+c^2-4*b-2*c+3, 3*a^2
+2*a*b-2*a*c-b^2-c^2-2*b-2*c+3, 3*a^2+2*a*b+2*a*c-b^2-c^2+2*b-2*c+3, 3*a^2+2*a
*b+4*a*c+b^2+c^2+4*b-2*c+3, 3*a^2+4*a*b-2*a*c+b^2+c^2-2*b-4*c+3, 3*a^2+4*a*b+2
*a*c+b^2+c^2+2*b-4*c+3];


text:
program IfstTC;
var ch : char;
    int : integer;
    boolx, booly : boolean;
begin
    boolx := true;
    booly := false;
    ch := 'a';
    int := 66;
    write(integer(ch));
    write(integer(int));
    write(integer(boolx), integer(booly));
    writeln;
    writeln(char(ch), char(int));
    writeln(boolean(ch), integer(int), integer(boolx), integer(booly));
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing

max_forks_repo_name: Hiroya-W/lang-processing


text:
# Class: Event
#
# Description:
#  Implementation of events -- a real algebraic number and list of associated quadrics.
#
# Author:
#  Kacper Pluta - kacper.pluta@esiee.fr
#  Laboratoire d'Informatique Gaspard-Monge - LIGM, A3SI, France
#
# Date:
#  29/12/2016
#
# License:
#  Simplified BSD License
#
# Copyright (c) 2016, Kacper Pluta
# All rights reserved.

# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   * Redistributions of source code must retain the above copyright
#     notice, this list of conditions and the following disclaimer.
#   * Redistributions in binary form must reproduce the above copyright
#     notice, this list of conditions and the following disclaimer in the
#     documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL Kacper Pluta BE LIABLE FOR ANY
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
EventType := module()
  option object;

  (*Real Algebraic Number*)
  local ranum::RealAlgebraicNumber;

  (*List of indices of associeted quadrics*)
  local quadrics::list;

# Method: ModuleCopy
#   Standard constructor / copy constructor
#
# Parameters:
#   self::EventType                   - a new object to be constructed
#   proto::EventType                  - a prototype object from which self is derived
#   ranum::RealAlgebraicNumber    - a real algebraic number
#   quadrics::rational            - a list of indices which corresponds to associated quadrics
#
# Output:
#   An object of type RealAlgebraicNumber.
#
# Exceptions:
#  "Invalid range. A range is valid when: a <= b."
#  "Degree of %1 is invalid."
#
  export ModuleCopy::static := proc( self::EventType,
                                     proto::EventType,
                                     ranum::RealAlgebraicNumber,
                                     quadrics::list, $ )
    if _passed = 2 then
      self:-ranum := proto:-ranum;
      self:-quadrics := proto:-quadrics;
    else
      self:-ranum := ranum;
      self:-quadrics := quadrics;
    fi;
  end proc;

# Method: ModulePrint
#   Standard printout of an object of type EventType.
#
# Parameters:
#   self::EventType                  - an event
#
  export ModulePrint::static := proc( self::EventType )
   nprintf("(%a, %a)", self:-ranum, self:-quadrics);
  end proc;

# Method: ModuleApply
#   Define standard constructor.
#
  export ModuleApply::static := proc()
   Object(EventType, args)
  end proc;

# Method: ModuleDeconstruct
#   Provides information how to recreate an object after being serialized.
#
# Parameters:
#   self::EventType                  - an event
#
  export ModuleDeconstruct::static := proc( self::EventType )
    ('EventType')(('RealAlgebraicNumber')(GetPolynomial(self:-ranum), GetInterval(self:-ranum)[1],
    GetInterval(self:-ranum)[2]), self:-quadrics)
  end proc;

# Method: GetRealAlgebraicNumber
#   Getter of the real algebraic number.
#
# Parameters:
#   self::EventType                  - an event
# Output:
#   Associated real algebraic number.
  export GetRealAlgebraicNumber::static := proc(self::EventType)
    return self:-ranum;
  end proc;

# Method: GetQuadrics
#   Getter of the list of indices.
#
# Parameters:
#   self::EventType                  - an event
# Output:
#   Associated quadrics' indices.
  export GetQuadrics::static := proc(self::EventType)
    return self:-quadrics;
  end proc;

# Method: Compare
#   A method used to compare two events.
#
# Parameters:
#   l::EventType      - an event
#   r::EventType      - an event
#
# Output:
#   -1 when l is smaller than r, 0 when they are equal and 1 when l is bigger than r.
#
  export Compare::static := proc( l::EventType, r::EventType, $ )
    return RealAlgebraicNumber:-Compare(l:-ranum, r:-ranum);
  end proc;

end module;

max_stars_repo_name: copyme/MapleTools

max_issues_repo_name: copyme/MapleTools

max_forks_repo_name: copyme/MapleTools


text:
with(BinomSums):
with(CodeTools):

### rser
Test(
  coeff(rser(1/(1-2*t), [t], 10),t,8),
  2^8);

Test(
  coeff(rser(1/(1-2*t)/u[1], [t,u[1]], 10),t,8),
  2^8);

Test(
  coeff(rser(1/(1-2*t)/u[1]/(u[2]-u[1]), [t,u[1],u[2]], 10),t,8),
  2^8);

### packvars
Test(
  BinomSums[packvars]([ 1/(x[4]+x[6]), x[7] ], x),
  [{x[4] = x[1], x[6] = x[2]}, {x[7] = x[3]}]
);

### geomsum
Test(
  BinomSums[geomsum](1, k=0..n),
  n+1
);

Test(
  BinomSums[geomsum](k, k=0..n),
  n*(n+1)/2,
  `expand`
);

# when the sum is infinite, geomsum add a variable
Test(
  BinomSums[geomsum](k, k=0..infinity),
  _W[_k]/(1-_W[_k])^2,
  `normal`
);

# when the sum is infinite on the left, geomsum yields an error
Test(
  BinomSums[geomsum](k, k=-infinity..0),
  "",
  `testerror`
);

# it should recognize when the input does not really depend on the summation variable
Test(
  BinomSums[geomsum]((-1+u[1])^(-n-1+j)*(-1+u[2])^(-n+j-1)*(1-u[1])^(-j)*(1-u[2])^(-j)*u[1]^(-n)*u[2]^(-n), j=0..n),
  u[1]^(-n)*u[2]^(-n)*(-1+u[1])^(-n-1)*(-1+u[2])^(-n-1)*(n+1),
  `normal`
);

max_stars_repo_name: lairez/binomsum

max_issues_repo_name: lairez/binomsum

max_forks_repo_name: lairez/binomsum


text:
(*
 Copyright (C) 2017 M.A.L. Marques
 Copyright (C) 2018 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_rtpss_params *params;

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

(* These are used within the tpss_x routine *)
tpss_ff     := z -> 2:
tpss_kappa := (x, t) -> params_a_kappa:

$include "tpss_x.mpl"

(* Equation (6) *)

rtpss_f := (x, u, t) -> 1 + tpss_kappa(x, t)*(1 - exp(-tpss_fx(x, t)/tpss_kappa(x,t))):

f := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) -> mgga_exchange(rtpss_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/.
*)

gvt4_gamm := (alpha, x, z) -> 1 + alpha*(x^2 + z):

gtv4 := (alpha, dd, x, z) ->
  dd[1]/gvt4_gamm(alpha, x, z) +
  (dd[2]*x^2 + dd[3]*z)/gvt4_gamm(alpha, x, z)^2 +
  (dd[4]*x^4 + dd[5]*x^2*z + dd[6]*z^2)/gvt4_gamm(alpha, x, z)^3:


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

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

(* type: lda_exc *)

(* This was checked against code kindly provided by Emil Proynov *)

(* Equations (24) *)
a_24_i := [-113.693369789727190, 24.00502151278711440, 49.34131295839670750,
            -23.8242372168379302, 0.944080741695104794, 0.000293039144178338]:
b_24_i := [-109.74263493216910, 16.2663129444242415, 54.4034331373908366,
            -25.154009904187990, 1.0]:

f_r := rs -> add(a_24_i[i]*rs^(i-1), i=1..6)/add(b_24_i[i]*rs^(i-1), i=1..5):

(* Equations (25) *)
c_25_i := [-0.32481568604919886, 1.180131465463191050, -1.42693041498421640,
            0.580344063812247980, -0.01099122367291440]:
d_25_i := [-0.57786103193239430, 2.09708505883490736, -2.52188183586948180, 1.0]:

f_s := z -> add(c_25_i[i]*z^(i-1), i=1..5)/add(d_25_i[i]*z^(i-1), i=1..4):

(* Equation (23) *)
(* The factor 1.28 is absent from the paper, but it is in the original code. See erratum *)
ss := (rs, z) -> f_r(rs)*f_s(z)*1.28:

(* Equation (22) *)
alpha_z := (rs, z) -> 2/(opz_pow_n(z,ss(rs, z)) + opz_pow_n(-z,ss(rs, z))):

(* Equation (21) *)
eta6  := 0.41081146652128:
eta7  := 0.599343256903515:
eta8  := 1.70939476802168:
eta9  := 0.077123208419481:
eta10 := 0.46958449007619:

alpha_n := rs ->
  + eta6
  + eta7*exp( -eta8*rs^(1/3))*rs^(2/3)
  + eta9*exp(-eta10*rs^(1/3))*rs^(1/3):

(* Equation (20) *)
alpha_eff := (rs, z) -> alpha_n(rs)*alpha_z(rs, z):

(* Equation (19) *)
eta1 :=  0.538074483500437:
eta2 := -2.226094990985190:
eta3 :=  0.837303782322808:
eta4 :=  2.619709858963178:
eta5 :=  1.036657594643520:

beta_eff := rs ->
  + eta1
  + eta2*exp(-eta3*rs^(1/3))*rs^(1/4)
  + eta4*exp(-eta5*rs^(1/3))*rs^(1/3):

(* Equation (15), see erratum *)
ax   := (3*Pi^2)^(1/3):
k_fs := (rs, z) -> ax*RS_FACTOR/rs * opz_pow_n(z,1/3):

(* Equation (17) *)
k_uu := (rs, z) -> alpha_eff(rs,  z)*k_fs(rs,  z):
k_dd := (rs, z) -> alpha_eff(rs, -z)*k_fs(rs, -z):

(* Equation (18) *)
k_ud := (rs, z) -> beta_eff(rs)
  * 2*k_fs(rs, z)*k_fs(rs, -z)/(k_fs(rs, z) + k_fs(rs, -z)):

(* Table III *)
a1  := 0.1846304394851914:
a2  := 5.93965654951900799:
a3  := 2.36958012866641818:
a4  := .51188865525958770e-1:
a5  := .9576892532004281e-1:
a6  := .283592616144882565e-1:
a7  := .226274169979695208e-1:
a8  := .531736155271654809e-2:
a9  := .1915378506400854:
a10 := .1473137771194929:
a11 := .1528250938350897:
a12 := 1.01508307543839117:
a13 := .7641254691754473e-1:
a14 := .898537460263473410:
a15 := .1795667349750801e-1:
a16 := .3461820740347690e-1:
a17 := .3591334699501599e-1:
a18 := .222017353476155799:

c1  := 132.479090287794355:
c2  := 32.4014708516771368:
c3  := 22.5664453162503806:
c4  := 11.2832226581251903:
c5  := .401060523940960082:
c6  := evalf(0.32):
c7  := .751988482389300153e-1:
c8  := 116.935042647480910:
c9  := 29.6240023046901289:
c10 := .482257181994472723:
c11 := .246903981179097557:
c12 := evalf(1/2):
c13 := .410709696778185459:
c14 := .105323524476768857:
c15 := 14.5650971711659670:
c16 := .781250000000000000:
c17 := .623347313127238558:
c18 := .146484374999999999:
c19 := 111.811548105797788:
c20 := .160041105570901272:
c21 := evalf(.78125):
c22 := .32086695060795739:
c23 := 13.2844495072998436:
c24 := .268418671319107341:
c25 := .471060597934991862:
c26 := evalf(1/4):
c27 := .252882919616989509:
c28 := .720485831127149779e-1:
c29 := 42.6490544891031073:

(* Definitions in the beginning of the appendix *)
D_1 := k -> a6*k^2 + a7*k + a8:
D_2 := k -> a1*k^2 + a10*k + a16:
D_3 := k -> a5*k^2 + a13*k + a15:
D_4 := k -> a9*k^2 + a11*k + a17:
D_5 := k -> c5*k^2 + c6*k + c7:
D_6 := k -> c12*k^2 + c13*k + c14:
D_7 := k -> c16*k^2 + c17*k + c18:
D_8 := k -> sqrt(c26*k^2 + c27*k + c28):

(* Equation (10) *)
Q_1ud := k ->  1/D_1(k) * (
  - arctan(a2*k + a3)*D_2(k)/k - log(D_1(k))*D_3(k)/k
  + log(k)*D_4(k)/k - a4*k + a12 + a14/k + a18/k^2
):

(* Equation (11) *)
Q_2ud := k ->
  - c1/k - c2/k^2 - c3*log(k)/k + c4*log(D_5(k))/k
  + c8*arctan(a2*k + a3)/k + c9*log(k + c10)/k - c11/k*log(D_6(k)):

(* Equation (12) *)
Q_3ud := k ->
  + c19*arctan(c20/(c21*k + c22))/k - c23*arctanh((c24 + c25*k)/D_8(k))/k
  - c15*log(D_7(k))/k - c29*D_8(k)/k^2:

(* Equation (9) *)
ec_opp := (rs, z) ->
  (1 - z^2)/4*(Q_1ud(k_ud(rs, z)) + Q_2ud(k_ud(rs, z)) + Q_3ud(k_ud(rs, z))):

(* Equation (13) *)
ec_par := (rs, z) ->
  + opz_pow_n( z,2)/8*(Q_1ud(k_uu(rs, z)) + Q_2ud(k_uu(rs, z)) + Q_3ud(k_uu(rs, z)))
  + opz_pow_n(-z,2)/8*(Q_1ud(k_dd(rs, z)) + Q_2ud(k_dd(rs, z)) + Q_3ud(k_dd(rs, z))):

f := (rs, z) -> n_total(rs)*(ec_opp(rs, z) + ec_par(rs, z)):


text:
######################################################################
# F(N)(A) is Inj(A,RR^N) modulo translation and scaling.
#
# Elements are represented as tables indexed by A, whose values
# are lists of length N of real numbers.  The equivalence relation
# is built into the definition of the function `is_equal/F` rather
# than the representation of the elements.

`is_element/F` := (N::posint) -> (A::set) -> proc(x)
 local a;
 global reason;

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

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

 for a in A do
  if not `is_element/R`(N)(x[a]) then
   reason := [convert(procname,string),"x[a] is not in R^N",a,x[a],N];
   return false;
  fi;
 od;

 if nops(map(op,{entries(x)})) <> nops(A) then
  reason := [convert(procname,string),"x is not injective",x];
  return false;
 fi;

 return true;
end;

######################################################################
# This function normalises the infinity norm.  This is convenient
# because it preserves rationality, but it is geometrically unpleasant.

`normalise/F` := (N::posint) -> (A::set) -> proc(x)
 local x0,u,r,a;

 x0 := table();
 u := [0$N];

 for a in A do
  x0[a] := x[a];
  u := u +~ x[a];
 od;
 u := u /~ nops(A);
 r := 0;
 for a in A do
  x0[a] := x0[a] -~ u;
  r := max(r,op(map(abs,x0[a])));
 od;
 for a in A do
  x0[a] := x0[a] /~ r;
 od;

 return(eval(x0));
end;

######################################################################
# This function normalises the euclidean norm
`normalise_2/F` := (N::posint) -> (A::set) -> proc(x)
 local x0,u,r2,r,a;

 x0 := table();
 u := [0$N];

 for a in A do
  x0[a] := x[a];
  u := u +~ x[a];
 od;
 u := u /~ nops(A);
 r2 := 0;
 for a in A do
  x0[a] := x0[a] -~ u;
  r2 := r2 + `norm_2/R`(N)(x0[a])^2;
 od;
 r := sqrt(r2);

 for a in A do
  x0[a] := x0[a] /~ r;
 od;

 return(eval(x0));
end;

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

`is_equal/F` := (N::posint) -> (A::set) -> proc(x,y)
 local x0,y0,a;
 global reason;

 x0 := `normalise/F`(N)(A)(x);
 y0 := `normalise/F`(N)(A)(y);

 for a in A do
  if not(`is_equal/R`(N)(x0[a],y0[a])) then
   reason := [convert(procname,string),"x0[a] <> y0[a]",a,x0[a],y0[a]];
   return false;
  fi;
 od;

 return true;
end;

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

`is_leq/F` := NULL;

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

`random_element/F` := (N::posint) -> (A::set) -> proc()
 local ok,x,a;

 ok := false;
 while not(ok) do
  x := table();
  for a in A do
   x[a] := `random_element/R`(N)();
  od:
  ok := `is_element/F`(N)(A)(x);
 od;

 return eval(x);
end;

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

`list_elements/F` := NULL;

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

`count_elements/F` := NULL;

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

`separation_infinity/F` := (N::posint) -> (A::set) -> proc(x)
 local x0,d,n,i,j;

 x0 := `normalise/F`(N)(A)(x);
 d := 2;
 n := nops(A);

 for i from 1 to n-1 do
  for j from i+1 to n do
   d := min(d,`d_infinity/R`(N)(x0[A[i]],x0[A[j]]));
  od;
 od;

 return d;
end;

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

`fatten/F` := (N::posint) -> (A::set) -> proc(x)
 local x0,d,m,a,i,e,u;

 x0 := `normalise/F`(N)(A)(x);
 d := `separation_infinity/F`(N)(A)(x0);
 m := 0;
 for a in A do
  for i from 1 to N do
   m := max(m,d/2+abs(x0[a][i]));
  od;
 od;

 e := [1 $ N];

 for a in A do
  u[a] := [((x0[a] -~ ((d/2) *~ e)) /~ m +~ e) /~ 2,
           ((x0[a] +~ ((d/2) *~ e)) /~ m +~ e) /~ 2];
 od;

 return eval(u);
end;

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

`inc/F/Fbar` := (N::posint) -> (A::set) -> proc(x)
 local y,P,T,t;

 y := table();

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

 for T in P do
  y[T] := table([seq(t = x[t],t in T)]);
 od:

 return eval(y);
end;

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

`gaps/F` := (N::posint) -> (A::set) -> proc(x)
 local G,V,W,m,i,j,a;

 G := NULL;
 m := nops(A);

 for i from 1 to N do
  V := sort([seq(x[a][i],a in A)]);
  W := [seq(V[j+1] - V[j],j = 1 .. m-1)];
  G := G,max(W);
 od;

 G := [G];
 return G;
end;

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

`gap/F` := (N::posint) -> (A::set) -> proc(x)
 return max(`gaps/F`(N)(A)(x));
end;

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

`normalise_gap/F` := (N::posint) -> (A::set) -> proc(x)
 local m,g,x0,a,i;

 m := [seq(min(seq(x[a][i],a in A)),i=1..N)];
 g := `gap/F`(N)(A)(x);

 x0 := table();
 for a in A do
  x0[a] := (x[a] -~ m) /~ g;
 od;

 return eval(x0);
end;

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

`delta/F` := (N::posint) -> (A::set) -> (a,b) -> proc(x)
 local u,n,i;
 u := x[b] -~ x[a];
 n := sqrt(add(u[i]^2,i=1..N));
 return u /~ n;
end:

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

`draw_unnormalised/F` := (A::set) -> proc(x)
 local a;
 display(seq(point(x[a]),a in A),scaling=constrained,axes=none);
end;

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

`draw/F` := (A::set) -> proc(x)
 display(`draw_unnormalised/F`(A)(`normalise/F`(2)(A)(x)),view=[-1..1,-1..1]);
end;


text:
######################################################################
# prime_simplex(A) is the set of maps x : A -> [0,1] with max = 1.

`is_element/prime_simplex` := (A::set) -> proc(x)
 local a,u,v;
 global reason;

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

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

 v := 0;

 for a in A do
  u := x[a];

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

  v := max(u,v);
 od;

 if simplify(v - 1) <> 0 then
  reason := [convert(procname,string),"max(x[a] : a in A) <> 1",a,v];
  return false;
 fi;

 return true;
end;

`is_equal/prime_simplex` := (A::set) -> proc(x,y)
 local a;
 global reason;

 for a in A do
  if simplify(x[a] - y[a]) <> 0 then
   reason := [convert(procname,string),"x[a] <> y[a]",a,x[a],y[a]];
   return false;
  fi;
 od;

 return true;
end;

`is_leq/prime_simplex` := NULL;

`random_element/prime_simplex` := (A::set) -> proc(d::posint := 5)
 local x,a,r1,r2,u;
 if A = {} then return FAIL; fi;

 r1 := rand(0..d);
 r2 := rand(1..d);
 x := table();
 u := 0;
 while u = 0 do
  for a in A do
   x[a] := r1()/r2();
   u := max(u,x[a]);
  od;
 od;

 for a in A do x[a] := x[a]/u; od;

 return eval(x);
end;

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

`phi/nonempty_subsets/prime_simplex` := (A::set) -> proc(U)
 local x,a;

 x := table();
 for a in A do x[a] := 0; od;
 for a in U do x[a] := 1; od;

 return eval(x);
end;

`supp/prime_simplex` := (A::set) -> proc(x)
 select(a -> x[a] > 0,A);
end:

`functor/prime_simplex` := (A::set,B::set) -> (f) -> proc(x)
 local a,b,y;

 y := table():
 for b in B do y[b] := 0; od:
 for a in A do y[f[a]] := max(y[f[a]],x[a]); od:
 return eval(y);
end:


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

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

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

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

$include "mgga_c_vsxc.mpl"
$include "mgga_c_m05.mpl"

m06l_f := (rs, z, xs0, xs1, ts0, ts1) ->
  + m05_f(rs, z, xs0, xs1, ts0, ts1)
  + vsxc_f(rs, z, xs0, xs1, ts0, ts1):

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


text:
    {   Diet.mpl   }

{
    This model determines a least cost diet which meets the daily
    allowances of nutrients for a man weighing 154 lbs.

    Reference:  G. B. Dantzig, Linear Programming and Extensions,
                Princeton University Press, Princeton, New Jersey, 1963.
}

TITLE
    Stiglers_nutrition_model

INDEX
    nutrients  =  (Calorie,Protein,Calcium,Iron,Avitamin
                   Thiamine,Riboflavin,Niacin,Cvitamin) : 10

    foods      =  (WheatFlour,CornMeal,EvapMilk,Margarine,
                   Cheese,Lard,Liver,PorkRoast,Salmon,
                   GreenBeans,Cabbage,Onions,Potatoes,Spinach,
                   SweetPot,Peaches,Prunes,LimaBeans,NavyBeans) : 10

DATA
    Required[nutrients]  =  (  3       ! Calories        [thousands]
                              70       ! Protein         [grams]
                               0.8     ! Calcium         [grams]
                              12       ! Iron            [milligrams]
                               5       ! Avitamin        [thousand ius]
                               1.8     ! Thiamine (B1)   [milligrams]
                               2.7     ! Riboflavin (B2) [milligrams]
                              18       ! Niacin          [milligrams]
                              75 ) ;   ! Cvitamin        [milligrams]

    A[foods,nutrients] = DATAFILE(nutri.dat)  !  Nutritive values of foods
                                              !  per dollar expenditure.

DECISION
    x[foods]   ->  ""    !  dollars of food to be purchased daily

MODEL

    MIN  Cost = SUM(foods: x) ;

SUBJECT TO

    NutrBal[nutrients] -> "" :  SUM(foods: A*x)  >  Required[nutrients];

END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


text:
# File: RigidMotionsRecoverNMM.mpl
#
# Description:
#  This file contains functions used to obtain an arrangement 6 dimensional parameter space of 3D
#  digitized rigid motions.
#  This code has been written for research propose and its aim is to calculate a particular
#  arrangement of quadrics. Therefore, it can or it cannot be useful in study of generic
#  arrangements. The final output are sample points of full dimensional open cells.
#
#  The code was written in relation with the paper: Kacper Pluta, Guillaume Moroz, Yukiko
#  Kenmochi, Pascal Romon, Quadric arrangement in classifying rigid motions of a 3D digital image,
#  2016, https://hal.archives-ouvertes.fr/hal-01334257 referred late on as [Quadrics:2016].
#
# Author:
#  Kacper Pluta - kacper.pluta@esiee.fr
#  Laboratoire d'Informatique Gaspard-Monge - LIGM, A3SI, France
#
# Date:
#  11/12/2015
#
# License:
#  Simplified BSD License
#
# Copyright (c) 2015, Kacper Pluta
# All rights reserved.

# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   * Redistributions of source code must retain the above copyright
#     notice, this list of conditions and the following disclaimer.
#   * Redistributions in binary form must reproduce the above copyright
#     notice, this list of conditions and the following disclaimer in the
#     documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL Kacper Pluta BE LIABLE FOR ANY
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#

RigidMotionsRecoverNMM := module()
  option package;
  uses   RigidMotionsParameterSpaceCommon, RigidMotionsMaplePrimesCode;

  (* String which represents a type of the neighborhood i.e. N1, N2 or N3. *)
  global nTypeGlobal;
  (* String which represents the path to the database. *)
  global dbPathGlobal;
  (* List of integers which are the indices of half-grid planes. *)
  global kRangeGlobal;
  (* List of the variables in which the problem is expressed. *)
  global varsGlobal;

  (*Controls how many sample points should be fetch from the database in one quary.*)
  export BUFFER_SIZE := 1000;

  export ParallelCalculateNMM, Get3DNMM, RecoverTranslationSamplePoints, GetOrderedCriticalPlanes,
  CriticalPlanes, CalculateNMM, LaunchFindDistinctSamplePoints,  LaunchComputeNMM,
  FetchSamplePointsFromDB, FetchTopologicallyDistinctSamplePointsFromDB,
  ParallelFindTopologicallyDistinctSamplePoints;

# Procedure: CriticalPlanes
#   Compute critical planes in the remainder range
#
# Parameters:
#   R                - the rotation matrix obtained from CayleyTransform
#   neighborhood     - a neighborhood for which one wants to compute NMM
#   kRange           - a range of planes to consider
#
# Output :
#   Returns a list of lists each containing critical planes for one direction
CriticalPlanes := proc(R::Matrix, neighborhood::list, kRange::list)
  # we remove the element [0, 0, 0]
  local n := subsop(ListTools:-Search([0,0,0], neighborhood)=NULL, neighborhood);
  local T := combinat:-cartprod([n, kRange]);
  local planes := [[],[],[]], params;
  while not T[finished] do
    params := T[nextvalue]();
    planes[1] := [op(planes[1]), (R . Vector(3, params[1]) - Vector(3, params[2]-1/2))[1]];
    planes[2] := [op(planes[2]), (R . Vector(3, params[1]) - Vector(3, params[2]-1/2))[2]];
    planes[3] := [op(planes[3]), (R . Vector(3, params[1]) - Vector(3, params[2]-1/2))[3]];
  end do;
 return planes;
end proc;

# Procedure: Get3DNMM
#   Compute neighbourhood motion maps
#
# Parameters:
#   neighborhood     - a neighborhood for which one wants to compute NMM
#   sampleTrans     - a midpoint of a frame in the remainder range
#   R                - the rotation matrix obtained from CayleyTransform
#
# Output:
#   Returns 3D neighborhood motion map for given vars and corresponding translations.
Get3DNMM := proc(neighborhood::list, sampleTrans::Vector, R::Matrix)
  local n, NMM := [];
  for n in neighborhood do
    NMM := [op(NMM), convert(map[inplace](round, R.Vector(3, n) + sampleTrans), list)];
  od;
  return NMM;
end proc;

# Procedure: RecoverTranslationSamplePoints
#   Compute midpoints of each frame in the remainder range
#
# Parameters:
#   planes            - ordered list of critical planes in the remainder range
#
# Output:
#   Returns centers of frames in the remainder range
RecoverTranslationSamplePoints := proc(planes::list)
  local s:
    s := proc(planes::list, i::integer, j::integer, k::integer)
      return [(1/2) * add(planes[1][i .. i+1]), (1/2)*add(planes[2][j.. j+1]),
                                            (1/2)*add(planes[3][k .. k+1])];
  end proc;
  return [seq(seq(seq(s(planes, i, j, k), k=1..nops(planes[3]-1)),
               j=1..nops(planes[2]-1)), i =1..nops(planes[1])-1)];
end proc;

# Procedure: GetOrderedCriticalPlanes
#   Compute critical planes in the remainder range
#
# Parameters:
#   vars             - a list of variables
#   samplePoint      - a rotational sample point
#   planes           - a precomputed list of planes in the remainder range
#
# Output:
#   Returns the ordered critical planes in the remainder range for X, Y and Z directions and order
#   signature.
GetOrderedCriticalPlanes := proc(vars::list, samplePoint::list, planes::list)
  local params, sdPlanes := [[],[],[]];
  local xSig, ySig, zSig, Signature;

  if nops(vars) <> 3 then
    error "Only 3D arrangement is supported.";
  fi;

  sdPlanes[1] := eval(planes[1], [vars[1] = samplePoint[1], vars[2] = samplePoint[2], vars[3] =
                                                                                samplePoint[3]]);
  sdPlanes[2] := eval(planes[2], [vars[1] = samplePoint[1], vars[2] = samplePoint[2], vars[3] =
                                                                                samplePoint[3]]);
  sdPlanes[3] := eval(planes[3], [vars[1] = samplePoint[1], vars[2] = samplePoint[2], vars[3] =
                                                                                samplePoint[3]]);

  xSig, sdPlanes[1] := Isort(sdPlanes[1]);
  ySig, sdPlanes[2] := Isort(sdPlanes[2]);
  zSig, sdPlanes[3] := Isort(sdPlanes[3]);

  # remove all out of the range regions
  sdPlanes[1] := remove(proc(x) return evalb(x <= -1/2) end proc,
  remove(proc(x) return evalb(x >= 1/2) end proc, sdPlanes[1]));
  sdPlanes[2] := remove(proc(x) return evalb(x <= -1/2) end proc,
  remove(proc(x) return evalb(x >= 1/2) end proc, sdPlanes[2]));
  sdPlanes[3] := remove(proc(x) return evalb(x <= -1/2) end proc,
  remove(proc(x) return evalb(x >= 1/2) end proc, sdPlanes[3]));
  # add region boarders
  sdPlanes[1] := [-1/2,op(sdPlanes[1]),1/2];
  sdPlanes[2] := [-1/2,op(sdPlanes[2]),1/2];
  sdPlanes[3] := [-1/2,op(sdPlanes[3]),1/2];

  Signature := cat(op(map(proc (x) sprintf("%d", x) end proc, [op(xSig), op(ySig), op(zSig)])));
  return Signature, sdPlanes;
end proc:

# Procedure: CalculateNMM
#   Reads data from hard drive and generates NMM
#
# Parameters:
#   vars             - list of variables in which the problem is expressed
#   planes           - a precomputed list of planes in the remainder range
#   buffer           - an Array which contains rotational sample points
#   N                - a given neighborhood
#   R                - a Matrix computed with Cayley transform
#   db               - an instance of ComputationRegister with open connection to the database
# Output:
#   A neighborhood motion map is saved in the database.
#
CalculateNMM := proc(vars::list, planes::list, buffer::Array, N::list, R::Matrix,
                                                         db::ComputationRegister)
  local sdPlanes, trans, x, RR, y, NMM;
  for x in buffer do
    RR := eval(R, [vars[1] = x[2], vars[2] = x[3], vars[3] = x[4]]);
    sdPlanes := GetOrderedCriticalPlanes(vars, x[2..()], planes)[2];
    trans := RecoverTranslationSamplePoints(sdPlanes);
    for y in trans do
      NMM := Get3DNMM(N, Vector(3, y), RR);
      InsertNMM(db, x[1], NMM, y);
    od;
  end do;
end proc:

# Procedure: FetchTopologicallyDistinctSamplePointsFromDB
#   Used to fetch topologically distinct rotational sample points from the database. The number of
#   sample points is controlled by BUFFER_SIZE.
#
# Parameters:
#   db::ComputationRegister       - an instance of ComputationRegister with open connection to the
#                                   database.
#   fromID::integer               - an id of the first sample point from a range to be fetched
#   last::integer                 - an id of the very last sample points to be fetched
#
# Output:
#   An Array of lists which of each represent a sample point. Note that the first element of each
#   list is an id of a given sample point.
FetchTopologicallyDistinctSamplePointsFromDB := proc(db::ComputationRegister, fromID::integer,
                                                                                last::integer)
  local n := fromID + RigidMotionsRecoverNMM:-BUFFER_SIZE - 1;
  if n > last then
    n := last;
  fi;
  return FetchTopologicallyDistinctSamplePoints(db, fromID, n);
end proc;

# Procedure: ParallelCalculateNMM
#   Uses Grid framework to generates unique NMM.
#
ParallelCalculateNMM := proc()
  local db:= Object(ComputationRegister, dbPathGlobal), n;
  local noTPoints;
  local first::integer, last::integer;
  local R := CayleyTransform(varsGlobal), N := GetNeighborhood(nTypeGlobal);
  local planes := RigidMotionsRecoverNMM:-CriticalPlanes(R, N, kRangeGlobal);
  local i::integer, buffer, samplePoint, sig::string;
  noTPoints := NumberOfTopologicallyDistinctSamplePoints(db);

  n := trunc(noTPoints / Grid:-NumNodes());
  first :=  Grid:-MyNode() * n + 1; last := (Grid:-MyNode() + 1) * n;
  if Grid:-MyNode() = Grid:-NumNodes() - 1 then
    last := noTPoints;
  fi;

  for i from first by RigidMotionsRecoverNMM:-BUFFER_SIZE to last do
    buffer := RigidMotionsRecoverNMM:-FetchTopologicallyDistinctSamplePointsFromDB(db, i, last);
    RigidMotionsRecoverNMM:-CalculateNMM(varsGlobal, planes, buffer, N, R, db);
    SynchronizeNMM(db);
  od;
  Close(db);
  Grid:-Barrier();
end proc:

# Procedure: FetchSamplePointsFromDB
#   Used to fetch rotational sample points from the database. The number of sample points is
#   controlled by BUFFER_SIZE.
#
# Parameters:
#   db::ComputationRegister       - an instance of ComputationRegister with open connection to the
#                                   database.
#   fromID::integer               - an id of the first sample point from a range to be fetched
#   last::integer                 - an id of the very last sample points to be fetched
#
# Output:
#   An Array of lists which of each represent a sample point. Note that the first element of each
#   list is an id of a given sample point.
FetchSamplePointsFromDB := proc(db::ComputationRegister, fromID::integer, last::integer)
  local n := fromID + RigidMotionsRecoverNMM:-BUFFER_SIZE - 1;
  if n > last then
    n := last;
  fi;
  return FetchSamplePointsWithoutSignature(db, fromID, n );
end proc;

# Procedure: ParallelFindTopologicallyDistinctSamplePoints
#   Finds rotational sample points which lead to unique arrangement of the critical planes in the
#   remainder range.
#
 ParallelFindTopologicallyDistinctSamplePoints := proc()
  local first::integer, last::integer;
  local R := CayleyTransform(varsGlobal), N := GetNeighborhood(nTypeGlobal);
  local planes := RigidMotionsRecoverNMM:-CriticalPlanes(R, N, kRangeGlobal);
  local db, n, i::integer, buffer, samplePoint, sig, noTPoints;
  db:= Object(ComputationRegister, dbPathGlobal);
  noTPoints := NumberOfSamplePoints(db);
  n := trunc(noTPoints / Grid:-NumNodes());
  first :=  Grid:-MyNode() * n + 1; last := (Grid:-MyNode() + 1) * n;
  if Grid:-MyNode() = Grid:-NumNodes() - 1 then
    last := noTPoints;
  fi;

  for i from first by RigidMotionsRecoverNMM:-BUFFER_SIZE to last do
    buffer := RigidMotionsRecoverNMM:-FetchSamplePointsFromDB(db, i, last);
    for samplePoint in buffer do
      sig := RigidMotionsRecoverNMM:-GetOrderedCriticalPlanes(varsGlobal,samplePoint[2..()],planes)[1];
      InsertSignature(db, samplePoint[1], sig);
    od;
  od;
  SynchronizeSamplePointsSignatures(db);
  Close(db);
  Grid:-Barrier();
end proc;

LaunchFindDistinctSamplePoints := proc(vars::list, nType::string, kRange::list, dbPath::string,
                                                        nodes:=kernelopts(numcpus))
  local db:=Object(ComputationRegister, dbPath);
  PrepareSamplePoints(db);
  Close(db);
  nTypeGlobal := nType; kRangeGlobal := kRange; dbPathGlobal := dbPath; varsGlobal := vars;

  Grid:-Setup("local");
  Grid:-Launch(RigidMotionsRecoverNMM:-ParallelFindTopologicallyDistinctSamplePoints,
  imports=['varsGlobal', 'nTypeGlobal', 'kRangeGlobal', 'dbPathGlobal'], numnodes=nodes,
  allexternal=false);
  db:=Object(ComputationRegister, dbPath);
  CloseSignaturesAddition(db);
  Close(db);
end proc:

# Procedure: LaunchOnGridGetNMM
#   Setup and run computation on a local grid
#
# Parameters:
#   vars             - list of variables in which the problem is expressed
#   nType            - size of neighborhood i.e N1, N2 and N3.
#   kRange           - a range of planes of the half-grid
#   dbPath           - a path to a database file.
#   nodes            - number of nodes used in the parallel computations
#
# Output:
#   List of unique neighborhood motion maps
LaunchComputeNMM := proc(vars::list, nType::string, kRange::list, dbPath::string,
                                                        nodes:=kernelopts(numcpus))
  local db:=Object(ComputationRegister, dbPath);
  DropRedundantSamplePoints(db);
  Close(db);
  nTypeGlobal := nType; kRangeGlobal := kRange; dbPathGlobal := dbPath; varsGlobal := vars;

  Grid:-Setup("local");
  Grid:-Launch(RigidMotionsRecoverNMM:-ParallelCalculateNMM, imports=['varsGlobal', 'nTypeGlobal',
  'kRangeGlobal', 'dbPathGlobal'], numnodes=nodes, allexternal=false);
end proc:

end module: