text:
`is_element/cactus_trees` := (A::set) -> proc(J)
 global reason;
 local tag,B,C,E,a,j;

 tag := "is_element/cactus_trees";

 if A = {} then
  reason := [tag,"A is empty"];
  return false;
 fi;

 if not(type(J,set(set)) and map(op,J) minus A = {}) then
  reason := [tag,"J is not a set of subsets of A",J,A];
  return false;
 fi;

 if min(op(map(nops,J))) < 2 then
  reason := [tag,"J contains a set of size < 2",J];
  return false;
 fi;

 B := {};
 C := {A[1]};
 while nops(B) < nops(C) do
  B := C;
  C := map(op,{B,map(op,select(X -> X intersect B <> {},J))});
 od;
 if C <> A then
  reason := [tag,"The graph is not connected",A,J,C];
  return false;
 fi;

 E := {seq(seq([a,j],a in j),j in J)};
 if nops(E) <> (nops(A) + nops(J) - 1) then
  reason := [tag,"The graph is not a tree",A,E];
  return false;
 fi;

return true;
end:

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

`is_leq/cactus_trees` := (A::set) -> proc(J0,J1)
 local f,G,j0,j1,jj1,P,Q,p,q,U;

 # f will be the map J0 -> J1 with j0 a subset of f(j0)
 f := table():

 # G(j1) = f^{-1}{j1}
 G := table():

 for j1 in J1 do G[j1] := {}; od:

 for j0 in J0 do
  jj1 := select(u -> (j0 minus u) = {},J1);
  # Assuming that J1 really is a cactus tree, the set jj1 will
  # have size 0 or 1.
  if nops(jj1) <> 1 then
   return false;
  fi;
  j1 := op(1,jj1);
  f[j0] := j1;
  G[j1] := {op(G[j1]),j0};
 od:

 for j1 in J1 do
  P := G[j1];
  p := nops(P);
  if p = 0 then return false; fi; # f should be surjective
  Q := {P[1]};
  q := 0;
  while q < nops(Q) do
   q := nops(Q);
   U := map(op,Q);
   Q := select(j0 -> ((j0 intersect U) <> {}),P);
  od:
  if Q <> P then
   return false;
  fi;
 od;

 return true;
end:

`top/cactus_trees` := proc(A::set)
 if nops(A) > 1 then
  return {A};
 elif nops(A) = 1 then
  return {};
 else
  return FAIL;
 fi;
end:

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

`random_element/cactus_trees` := (A::set) -> proc()
 local a,b,B,J,K,k;

 if nops(A) = 0 then
  return FAIL;
 elif nops(A) = 1 then
  return {};
 elif nops(A) = 2 then
  return {A};
 fi;

 a := random_element_of(A);
 B := A minus {a};
 K := `random_element/cactus_trees`(B)();
 if rand(0..1)() = 0 then
  k := random_element_of(K);
  J := K minus {k} union {k union {a}};
 else
  b := random_element_of(B);
  J := K union {{a,b}};
 fi;

 return J;
end:

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

`list_elements/cactus_trees` := proc(A::set)
 option remember;
 local n,a,b,u,v,B,PP,pp,K,L,m,M,P,Q,i,k,f;

 n := nops(A);
 if n = 0 then
  return [];
 elif n = 1 then
  return [{}];
 elif n = 2 then
  return [{A}];
 fi;

 a := A[1];
 B := A minus {a};
 PP := map(pp -> [op(pp)],`list_elements/partitions`(B));
 L := NULL;
 for pp in PP do
  m := nops(pp);
  M := map(`list_elements/cactus_trees`,pp);
  P := [seq([seq(seq([K,{b}],b in pp[i]),K in M[i]),
             seq(seq([K,k],k in K),K in M[i])],i=1..m)];
  Q := [[]];
  for i from 1 to m do
   Q := [seq(seq([op(u),v],u in Q),v in P[i])];
  od:
  f := (q) -> {seq(op(q[i][1] minus {q[i][2]}),i=1..m),seq({a,op(q[i][2])},i=1..m)};
  L := L,op(map(f,Q));
 od:
 return [L];
end:

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

# This is A030019 in OEIS

`count_elements/cactus_trees` := proc(A::set)
 local n,i;
 n := nops(A);
 add(Stirling2(n-1, i)*n^(i-1),i=0..n-1);
end:


text:
`is_element/morse_matchings` := (A::set) -> (R) -> proc(M)
 global reason;
 local e,H,F;

 if not(type(M,set)) then
  reason := ["is_element/morse_matchings","M is not a set",M];
  return false;
 fi;

 for e in M do
  if not(type(e,list)) and nops(e) = 2 then
   reason := ["is_element/morse_matchings","M is not a set of lists of length 2",M,e];
   return false;
  fi;

  if {op(e)} minus A <> {} then
   reason := ["is_element/morse_matchings","M is not a set of lists of length 2 in A",M,e];
   return false;
  fi;
 od:

 H := `hasse_diagram/partord`(A)(R);

 if M minus H <> {} then
  reason := ["is_element/morse_matchings","M is not a subset of the Hasse diagram",M,H];
  return false;
 fi;

 for e in M do
  if select(f -> {op(f)} intersect {op(e)} <> {}, M) <> {e} then
   reason := ["is_element/morse_matchings","edges in M are not disjoint",M,e];
   return false;
  fi;
 od:

 F := M union map(e -> [e[2],e[1]],H minus M);
 F := `transitive_closure/autorel`(A)(F);
 F := F intersect `op/autorel`(A)(F);
 if F minus `id/autorel`(A) <> {} then
  reason := ["is_element/morse_matchings","M has cycles",M,F];
  return false;
 fi;

 return true;
end:

ext: mpl


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

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

(* type: work_mgga_x *)
(* prefix:
  mgga_x_m08_params *params;

  assert(pt->params != NULL);
  params = (mgga_x_m08_params * ) (pt->params);
*)

params_a_rpbe_kappa := 0.552:
params_a_rpbe_mu    := MU_GE:
$include "gga_x_rpbe.mpl"

params_a_kappa := KAPPA_PBE:
params_a_mu    := 0.21951:
$include "gga_x_pbe.mpl"

m08_f := (a, b, x, t) ->
  + f_pbe(x) *mgga_series_w(a, 12, t)
  + f_rpbe(x)*mgga_series_w(b, 12, t):

f_m08 := (x, t) ->
  m08_f(params_a_a, params_a_b, x, t):

f := (rs, x, t, u) ->
  f_m08(x, t):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
with(Groebner):
with(PolynomialIdeals):

J := PolynomialIdeal({x0^2 - x0 + 2*x1^2 + 2*x2^2 + 2*x3^2 + 2*x4^2 + 2*x5^2 + 2*x6^2 + 2*x7^2 + 2*x8^2,
2*x0*x1 + 2*x1*x2 - x1 + 2*x2*x3 + 2*x3*x4 + 2*x4*x5 + 2*x5*x6 + 2*x6*x7 + 2*x7*x8, 2*x0*x2 + x1^2 + 2*x1*x3 + 2*x2*x4 - x2 + 2*x3*x5 + 2*x4*x6 + 2*x5*x7 + 2*x6*x8,
2*x0*x3 + 2*x1*x2 + 2*x1*x4 + 2*x2*x5 + 2*x3*x6 - x3 + 2*x4*x7 + 2*x5*x8,
2*x0*x4 + 2*x1*x3 + 2*x1*x5 + x2^2 + 2*x2*x6 + 2*x3*x7 + 2*x4*x8 - x4,
2*x0*x5 + 2*x1*x4 + 2*x1*x6 + 2*x2*x3 + 2*x2*x7 + 2*x3*x8 - x5,
2*x0*x6 + 2*x1*x5 + 2*x1*x7 + 2*x2*x4 + 2*x2*x8 + x3^2 - x6,
2*x0*x7 + 2*x1*x6 + 2*x1*x8 + 2*x2*x5 + 2*x3*x4 - x7,
x0 + 2*x1 + 2*x2 + 2*x3 + 2*x4 + 2*x5 + 2*x6 + 2*x7 + 2*x8 - 1}):
st := time[real]():
Basis(J, tdeg(x0,x1,x2,x3,x4,x5,x6)):
time[real]() - st;


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

sys := [
diff(S(t), t) = -bi*S(t)*II(t) - S(t)*mu - S(t)*bw*W(t) + R(t)*a + mu,
diff(II(t), t) = bi*S(t)*II(t) - gam*II(t) + S(t)*bw*W(t) - mu*II(t),
diff(W(t), t) = xi*II(t) - xi*W(t),
diff(R(t), t) = gam*II(t) - R(t)*mu - R(t)*a,
y(t) = k*II(t),
y2(t) = S(t) + R(t) + II(t)
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


text:
# Roboter aus Viergelenkkette und Schubgelenk
# Voraussetzung:
# Viergelenkkette muss mit der Toolbox berechnet worden sein (Arbeitsblatt "vgk_kinematic_constraints.mw")
# Autor
# Abderahman Bejaoui
# Studienarbeit bei: Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2018-08
# (C) Institut fuer mechatronische Systeme, Leibniz Universitaet Hannover
# Quellen
# TODO
#
# Initialisierung
# Import
restart:
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
#with(ListTools):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions",robot_name):
# Variable mit Winkeln der Nebenstruktur nur in Abhängigkeit der verallgemeinerten Koordinaten
kintmp_qs := Matrix(RowDimension(kintmp_s),1):
qJ_t:= <qJ1(t),qJ2(t),qJ3(t),qJ4(t),qJ5(t)>:
qJ_s:= <qJ1s,qJ2s,qJ3s,qJ4s,qJ5s>:
# Konstante Winkel bereits hineinschreiben
for i from 1 to RowDimension(kintmp_s) do
  if diff(kintmp_s(i,1), t) = 0 then
    kintmp_qs(i,1) := kintmp_s(i,1):
  end if:
end do:
# Variablen definieren für die Hilfswinkel
#
# Ersetzungsausdrücke definieren.
# Speichere Sinus und Cosinus der Winkel direkt ab, da diese in den Rotationsmatrizen direkt auftreten.
# Spalte 1: Zu suchender Ausdruck (sin oder cos eines Winkels)
# Spalte 2: Einzusetzender Ausdruck.
# Dadurch werden arctan-Ausdrücke in der direkten Kinematik reduziert.
# Ähnliches Vorgehen wie in [1].
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
for i from 1 to RowDimension(kintmp_s) do
  kintmp_subsexp(2*i-1, 1) := sin(kintmp_s(i,1)):
  kintmp_subsexp(2*i,   1) := cos(kintmp_s(i,1)):
  # Initialisierung der rechten Spalte mit gleichen Werten. Später nur Ersetzung, wenn Vorteilhaft.
  kintmp_subsexp(2*i-1, 2) := kintmp_subsexp(2*i-1, 1):
  kintmp_subsexp(2*i,   2) := kintmp_subsexp(2*i,   1):
end do:
read sprintf("../codeexport/fourbar1TE/tmp/fourbar1TE_angles_trig_elim"):
# winkel
;
winkel(..,2):
# Viergelenkkette : A-B-C-D
winkel1:=<<sin_rho34|sin(rho34_s)>;<cos_rho34|cos(rho34_s)>;<sin_rho12|sin(rho12_s)>;<cos_rho12|cos(rho12_s)>;<sin_rho1|sin(rho1_s)>;<cos_rho1|cos(rho1_s)>>:
for i from 1 to 6 do
	winkel1(i,2):=winkel(i,2):
     cos_phi_s :=sin(qJ5s):
	sin_phi_s := cos(qJ5s):
	winkel1(i,2):=subs({sin(qJ1s)=sin_phi_s},winkel1(i,2)):
	winkel1(i,2):=subs({cos(qJ1s)=cos_phi_s},winkel1(i,2)):
     winkel1(i,2) := subs({l1=AB}, winkel1(i,2)):
     winkel1(i,2) := subs({l2=AD}, winkel1(i,2)):
     winkel1(i,2) := subs({l3=qJ2s+HC}, winkel1(i,2)):
     winkel1(i,2) := subs({l4=CB}, winkel1(i,2)):
end do:
sin_rho34_s:=winkel1(1,2):
cos_rho34_s:=winkel1(2,2):
sin_rho12_s:=winkel1(3,2):
cos_rho12_s:=winkel1(4,2):
sin_rho1_s:=winkel1(5,2):
cos_rho1_s:=winkel1(6,2):
sin_rho48_s:=sin_rho1_s:
cos_rho48_s:=cos_rho1_s:
# Kitmp_subsexp
kintmp_subsexp(1,2):=sin_rho12_s:
kintmp_subsexp(2,2):=cos_rho12_s:
kintmp_subsexp(3,2):=sin_rho34_s:
kintmp_subsexp(4,2):=cos_rho34_s:
kintmp_subsexp(5,2):=sin_rho48_s:
kintmp_subsexp(6,2):=cos_rho48_s:
# Kintmp_qs
kintmp_qs(1,1):=arctan(sin_rho12_s,cos_rho12_s):
kintmp_qs(2,1):=arctan(sin_rho34_s,cos_rho34_s):
kintmp_qs(3,1):=arctan(sin_rho48_s,cos_rho48_s):
# Export
kintmp_qt := convert_s_t(kintmp_qs):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name):
#printf("Ausdrücke für kintmp_subsexp gespeichert (Maple). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
for i from 1 to RowDimension(kintmp_s) do
  tmp := kintmp_qs(i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d",robot_name, i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d.m", robot_name, i):
end do:
save kin_constraints_exist, kintmp_qs, kintmp_qt,kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert" ,robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
save kintmp_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_qs_maple_inert", robot_name):
#printf("Ausdrücke mit Inert-Arctan exportiert (Matlab). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp ));
#kc_symbols :=Transpose(kc_symbols);
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m",robot_name),2);

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


text:
# Inverse Dynamics for Robot with IC
# Einleitung
# Berechnung der inversen Dynamik für Roboter mit impliziten Zwangsbedingungen
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# implicit_constraints -> Berechnung für Roboter mit impliziten Zwangsbedingungen
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# dynamics -> Berechnung der Dynamik
# worldframe -> Berechnung basierend auf Energien aus Welt-KS (KS W)
# par1 -> Parametersatz 1 (Schwerpunkt als Parameter: SX,SY,SZ)
# reg2 -> Regressorform
# Autor
# Tim Job (HiWi-Job bei Moritz Schappler), 2019-05
# Moritz Schappler, moritz.schappler@imes.uni-hannover.de
# (C) Institut für Mechatronische Systeme, Universität Hannover
# Sources
# [ParkChoPlo1999] Park, FC and Choi, Jihyeon and Ploen, SR: Symbolic formulation of closed chain dynamics in independent coordinates
# [Docquier2013] Docquier, Nicolas and Poncelet, Antoine and Fisette, Paul: ROBOTRAN: a powerful symbolic gnerator of multibody models (2013)
# [DoThanhKotHeiOrt2009b] Do Thanh et al.: On the inverse dynamics problem of general parallel robots (2009)
# Initialization
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
# Einstellungen für Code-Export: Optimierungsgrad (2=höchster) und Aktivierung jedes Terms.
#codegen_act := true: # noch nicht implementiert
codegen_debug := false:
codegen_opt := 2:
read "../helper/proc_MatlabExport":
read "../robot_codegen_definitions/robot_env_IC":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name_OL):
read "../robot_codegen_definitions/robot_env_IC":
if base_method_name="twist" then # Basis-Methode "twist" wird (hier) nur für fixed Base benutzt
  expstring:="fixb":
elif base_method_name="eulxyz" then
  expstring:="floatb_eulxyz":
else
  printf("Nicht behandelte Basis-Methode: %s\n", base_method_name):
fi:
# Kennung des Parametersatzes, für den die Dynamikfunktionen erstellt werden sollen. Muss im Repo und in der mpl-Datei auf 1 gelassen werden, da die folgende Zeile mit einem Skript verarbeitet wird.
codegen_dynpar := 1:
regressor_modus := "regressor":
herleitungsverfahren := "lagrange":
codeexport_invdyn := false:
codeexport_grav_inertia := false:
codeexport_coriolisvec := false:
codeexport_regressor := false:
if herleitungsverfahren = "lagrange" then #Dynamik aus Lagrange
# Ergebnisse der inversen Dynamik für par1/2 laden
  if codeexport_invdyn then
    read sprintf("../codeexport/%s/tmp/invdyn_fixb_par%d_maple.m", robot_name_OL, codegen_dynpar):
    tau := Matrix(taus_fixb(7..NQ,1)):
    printf("%s. Generiere inverse Dynamik (%s) für %s basierend auf IC-Jacobi und OL-Dynamik mit Parametersatz %d  und %s\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), herleitungsverfahren, robot_name, codegen_dynpar, base_method_name):
  fi:
# Ergebnisse der inversen Dynamik in Regressorform laden
  if codeexport_regressor then
    read sprintf("../codeexport/%s/tmp/invdyn_%s_%s_maple.m", robot_name_OL, expstring, regressor_modus):
    tau_regressor_s := tau_regressor_s(7..NQ,..):
    printf("%s. Generiere Dynamik in Regressorform (%s) für %s basierend auf IC-Jacobi und OL-Dynamik\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), herleitungsverfahren, robot_name):
  fi:
# Ergebnisse der Massenmatrix und des Gravitationsvektors laden
  if codeexport_grav_inertia then
    read sprintf("../codeexport/%s/tmp/gravload_par%d_maple.m", robot_name_OL, codegen_dynpar):
    taug := Matrix(taug_s(7..NQ,..)):
    read sprintf("../codeexport/%s/tmp/inertia_par%d_maple.m", robot_name_OL, codegen_dynpar):
    MM := Matrix(MM_s(7..NQ,7..NQ)):
    read sprintf("../codeexport/%s/tmp/coriolisvec_par%d_maple.m", robot_name_OL, codegen_dynpar):
    tauCC := Matrix(tauCC_s(7..NQ,..)):
    printf("%s. Generiere Gravitationsvektor und Massenmatrix (%s) für %s basierend auf IC-Jacobi und OL-Dynamik mit Parametersatz %d und %s\n",\
      FormatTime("%Y-%m-%d %H:%M:%S"), herleitungsverfahren, robot_name, codegen_dynpar, base_method_name):
  fi:

# Ergebnisse des Coriolisvektors laden. Noch nicht implementiert.
  if codeexport_coriolisvec then
    read sprintf("../codeexport/%s/tmp/coriolisvec_par%d_maple.m", robot_name_OL, codegen_dynpar):
    tauCC := Matrix(tauCC_s(7..NQ,..)):
    read sprintf("../codeexport/%s/tmp/inertia_par%d_maple.m", robot_name_OL, codegen_dynpar):
    MM := Matrix(MM_s(7..NQ,7..NQ)):
    printf("%s. Generiere Coriolisvektor (%s) für %s basierend auf IC-Jacobi und OL-Dynamik mit Parametersatz %d und %s\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), herleitungsverfahren, robot_name, codegen_dynpar, base_method_name):
  fi:
else #Dynamik aus Newton-Euler
# Ergebnisse der inversen Dynamik für par1/2 laden
  invdynfile:=sprintf("../codeexport/%s/tmp/invdyn_twist_NewtonEuler_linkframe_par%d_maple.m", robot_name_OL, codegen_dynpar):
  if FileTools[Exists](invdynfile) then
    read invdynfile:
  else
    printf("%s. Inversdynamik der offenen Kette wurde nicht berechnet. Abbruch der Newton-Euler-Berechnung (Regressorform).\n", \
      FormatTime("%Y-%m-%d %H:%M:%S")):
    quit: # Funktioniert in GUI nicht richtig...
    robot_name := "": # ...Daher auch Löschung des Roboternamens.
  end if:
  tau := Matrix(tau_J(1..NQJ,..)):
  taug := copy(tau):
  for i from 1 to NQ do
    taug := subs({qD_s[i,1]=0},taug):
    taug := subs({qDD_s[i,1]=0},taug):
  end do:
  printf("%s. Generiere Gravitationsvektor und gesamte inverse Dynamik (%s) für %s basierend auf IC-Jacobi und OL-Dynamik mit Parametersatz %d\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), herleitungsverfahren, robot_name, codegen_dynpar):
fi:
# Ergebnisse der impliziten Zwangsbedingungen laden
read sprintf("../codeexport/%s/tmp/kinconstr_impl_projection_jacobian_maple", robot_name):
B21 := B21:
read sprintf("../codeexport/%s/tmp/kinconstr_impl_projection_jacobian_derivative_maple", robot_name):
B21D := B21D:
# Ergebnisse des Positionsvektors und der aktiven und passiven Gelenke laden
read sprintf("../codeexport/%s/tmp/positionVector_NQJ_maple.m", robot_name):
posNQJ := posNQJ:
# Pertubationsmatrizen berechnen
# Entferne aktive/passive Koordinaten zur Umrechnung
# Berechne die Anzahl aller aktiven und passiven Gelenke (inkl virtueller Schnittgelenke)
NAJ := add(mu(k), k=1..NJ):
NPJ := NJ-NAJ:
# Bestimme Indizes der aktiven und passiven Gelenke
# Entspricht Partitionierung der Gelenkwinkel in q1 und q2 in [Docquier2013]. Hier kann die Reihenfolge von aktiven und passiven Winkeln auch durchmischt sein.
IndAct := Matrix(NAJ, 1):
IndPass:= Matrix(NPJ, 1):
ka := 1: kp := 1:
for i from 1 to NJ do
  if mu(i) = 1 then
    IndAct(ka) := i:
    ka := ka + 1:
  else
    IndPass(kp) := i:
    kp := kp + 1:
  end if:
end do:
# Entfernt die Einträge aus IndPass/IndAct, die virtuelle Schnittglenke sind
for i from 1 to NJ do
  if posNQJ(i) = 0 then
    posAct := ListTools[Search](i, convert(IndAct, list)):
    if posAct <> 0 then
      IndAct := DeleteRow(IndAct,posAct):
    end if:
    #posPass := ListTools[Search](i, convert(IndPass, list)):
    #if posPass <> 0 then
    #  IndPass := DeleteRow(IndPass,posPass):
    #end if:
  end if:
end do:
NQJA := RowDimension(IndAct):
NQJP := NQJ-NQJA: #RowDimension(IndPass):
IndPass:= Matrix(NQJP, 1):
kp := 1: ka := 1:
for i from 1 to NQJ do
  if ListTools[Search](i, convert(IndAct, list)) = 0 then
    IndPass(kp) := i:
    kp := kp + 1:
  end if:
end do:
# Permutationsmatrizen berechnen zur Umrechnung von aufgeteilten Koordinaten
# (aktiv/passiv nacheinander) in Koordinaten der Baumstruktur (aktiv/passiv gemischt)
P1 := Matrix(NQJA,NQJ):
P2 := Matrix(NQJP,NQJ):
for i from 1 to NQJA do
  P1(i,IndAct(i)) := 1:
end do:
for i from 1 to NQJP do
  P2(i,IndPass(i)) := 1:
end do:
# Inverse Dynamik
# explizit [Docquier2013], Gl. 12
W := Transpose(P1 + Transpose(B21).P2):
WD := Transpose(Transpose(B21D).P2):
#q1D_s := Transpose(W).qJD_s:
tauIC := Transpose(W).tau:
taugIC := Transpose(W).taug:
# [DoThanhKotHeiOrt2009b], Gl. (23); [ParkChoPlo1999] Gl. (56)
#tauCCIC := Transpose(W).tauCC + Transpose(W).MM.WD.q1D_s: # noch nicht implementiert
;
tauIC_regressor := P1.tau_regressor_s + Transpose(B21).P2.tau_regressor_s:
if herleitungsverfahren = "newton" then
  tau_MM := copy(tau):
  for i to NQ do
    tau_MM := subs(qD_s[i, 1] = 0,tau_MM);
  end do:

  MM := Matrix(NQJ, NQJ):
  for i to NQJ do
    for j to NQJ do
      MM[i, j] := diff(tau_MM[i, 1], qDD_s[6+j, 1]):
    end do:
  end do:
end if:
# [ParkChoPlo1999], Gl. 55; [DoThanhKotHeiOrt2009b], Gl. (23)
MMIC := Transpose(W).MM.W:
# Export
# Export der Belastung der Gelenke
# Floating Base
if codeexport_invdyn and not(base_method_name="twist") then
  MatlabExport(tauIC, sprintf("../codeexport/%s/tmp/invdyn_floatb_%s_par%d_ic_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
if codeexport_regressor and not(base_method_name="twist") then
  MatlabExport(tauIC_regressor, sprintf("../codeexport/%s/tmp/invdyn_floatb_%s_%s_ic_matlab.m", robot_name, base_method_name, regressor_modus), codegen_opt):
end if:
# Fixed Base
if codeexport_invdyn then
  tauIC_fixb:=tauIC:
  for i from 1 to NQB do
    tauIC_fixb := subs({X_base_s[i,1]=0},tauIC_fixb):
  end do:
  for i from 1 to 6 do
    tauIC_fixb := subs({V_base_s[i,1]=0},tauIC_fixb):
    tauIC_fixb := subs({VD_base_s[i,1]=0},tauIC_fixb):
  end do:
  MatlabExport(tauIC_fixb, sprintf("../codeexport/%s/tmp/invdyn_fixb_par%d_ic_matlab.m", robot_name, codegen_dynpar), codegen_opt):
end if:
if codeexport_grav_inertia then
  taugIC_fixb:=taugIC:
  MMIC_fixb:=MMIC:
  for i from 1 to NQB do
    MMIC_fixb := subs({X_base_s[i,1]=0},MMIC_fixb):
  end do:
  for i from 1 to 6 do
    MMIC_fixb := subs({V_base_s[i,1]=0},MMIC_fixb):
    MMIC_fixb := subs({VD_base_s[i,1]=0},MMIC_fixb):
  end do:
  MatlabExport(taugIC_fixb, sprintf("../codeexport/%s/tmp/gravload_joint_floatb_twist_par%d_ic_matlab.m", robot_name, codegen_dynpar), codegen_opt):
  MatlabExport(MMIC_fixb, sprintf("../codeexport/%s/tmp/inertia_joint_joint_floatb_%s_par%d_ic_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Noch nicht implementiert
if codeexport_coriolisvec then
  tauCCIC_fixb:=tauCCIC:
  for i from 1 to NQB do
    tauCCIC_fixb := subs({X_base_s[i,1]=0},tauCCIC_fixb):
  end do:
  for i from 1 to 6 do
    tauCCIC_fixb := subs({V_base_s[i,1]=0},tauCCIC_fixb):
    tauCCIC_fixb := subs({VD_base_s[i,1]=0},tauCCIC_fixb):
  end do:
  MatlabExport(tauCCIC_fixb, sprintf("../codeexport/%s/tmp/coriolisvec_joint_fixb_twist_par%d_ic_matlab.m", robot_name, codegen_dynpar), codegen_opt):
end if:
if codeexport_regressor then
  tauIC_regressor_fixb:=tauIC_regressor:
  for i from 1 to NQB do
    tauIC_regressor_fixb := subs({X_base_s[i,1]=0},tauIC_regressor_fixb):
  end do:
  for i from 1 to 6 do
    tauIC_regressor_fixb := subs({V_base_s[i,1]=0},tauIC_regressor_fixb):
    tauIC_regressor_fixb := subs({VD_base_s[i,1]=0},tauIC_regressor_fixb):
  end do:
  MatlabExport(tauIC_regressor_fixb, sprintf("../codeexport/%s/tmp/invdyn_fixb_%s_ic_matlab.m", robot_name, regressor_modus), codegen_opt):
end if:
if herleitungsverfahren = "newton" then
  tauIC_fixb:=tauIC:
  taugIC_fixb:=taugIC:
  MMIC_fixb:=MMIC:
  for i from 1 to NQB do
    tauIC_fixb := subs({X_base_s[i,1]=0},tauIC_fixb):
    MMIC_fixb := subs({X_base_s[i,1]=0},MMIC_fixb):
  end do:
  for i from 1 to 6 do
    tauIC_fixb := subs({V_base_s[i,1]=0},tauIC_fixb):
    tauIC_fixb := subs({VD_base_s[i,1]=0},tauIC_fixb):
    MMIC_fixb := subs({V_base_s[i,1]=0},MMIC_fixb):
    MMIC_fixb := subs({VD_base_s[i,1]=0},MMIC_fixb):
  end do:
  MatlabExport(tauIC_fixb, sprintf("../codeexport/%s/tmp/invdyn_fixb_snew_par%d_ic_matlab.m", robot_name, codegen_dynpar), codegen_opt):
  MatlabExport(taugIC_fixb, sprintf("../codeexport/%s/tmp/gravload_joint_floatb_twist_snew_par%d_ic_matlab.m", robot_name, codegen_dynpar), codegen_opt):
  MatlabExport(MMIC_fixb, sprintf("../codeexport/%s/tmp/inertia_joint_joint_floatb_%s_snew_par%d_ic_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Export Anzahl der aktiven Koordinaten NAJ
MatlabExport(NAJ, sprintf("../codeexport/%s/tmp/NAJ_ic_matlab.m", robot_name), codegen_opt):

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


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

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

(* type: work_gga_x *)

mu    := 0.2195149727645171:
c     := (146/2025)*(4/9) - (73/405)*(2/3) + (mu - 10/81):
kappa := 0.8040:

f0_aux := s -> kappa + MU_GE * s^2 + (mu - MU_GE) * s^2 * exp(-s^2) + log(1 + c * s^4):
f0_wc  := s -> 1 + kappa*(1 - kappa/f0_aux(s)):
f_wc   := x -> f0_wc(X2S*x):

f      := x -> f_wc(x):


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

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

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

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

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

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

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

       newJ
  end proc;

  #main procedure for int/sum reparamterizations
  reparam:= proc(e::LO(symbol, algebraic), {ctx::list:= []}, $)
  local
       J:= innermostIntSum(e),   #the integral or sum
       newJ, #What J will become
       #possible subs target(s)
       oldarg::{algebraic, set(algebraic)}:= map2(op, 2, indets(e, specfunc(applyintegrand)))
  ;
       if nops(J) = 0 then
            userinfo(2, procname, "No sum or integral found.");
            return e
       end if;
       if nops(J) > 1 then
            userinfo(1, reparam, "Case of multiple innermost Int or Sums not yet handled.");
            return 'procname'(e)
       end if;
       J:= J[];
       if J::specfunc(Sum) then
            userinfo(1, reparam, "Sums not yet handled.");
            return 'procname'(e)
       end if;

       if nops(oldarg) <> 1 then
            userinfo(1, thisproc, "More than 1 reparam possible:", oldarg);
            return 'procname'(e)
       end if;
       oldarg:= oldarg[];   #Extract the reparam target.

       #If the target is simply a name, return input unchanged.
       if oldarg::symbol then
            userinfo(2, procname, "No need for a reparameterization.");
            return e
       end if;

       #Make the change of vars.
       newJ:= simplify_assuming('ChangeVarInt(J, oldarg)', build_kb(ctx,"reparam"));

       if newJ = 0 then
            WARNING("Integral is 0, likely due to improper handling of an infinity issue.");
            userinfo(
                 1, procname, "u subs:",
                 print(
                      #Reformat the ChangeVarInt command for readability.
                      subs(
                           x= ':-x',
                           'ChangeVarInt'(
                                subsindets(
                                     J,
                                     specfunc(applyintegrand),
                                     f-> ':-h'(op(2,f))
                                )
                           ),
                           oldarg
                      )
                 )
            )
       end if;

       subs(J= newJ, e)
  end proc;

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
program Write;
begin
    ;
    writeln('It'

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:
# 4-Achsroboter von KUKA
# Voraussetzung:
# Viergelenkkette muss mit der Toolbox berechnet worden sein (Arbeitsblatt "vgk_kinematic_constraints.mw")
# Autor
# Weipu Shan
# Studienarbeit bei: Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2018-08
# (C) Institut fuer mechatronische Systeme, Leibniz Universitaet Hannover
# Quellen
# SA Bejaoui: Bejaoui2018_S749; "Modellierung kinematischer Zwangsbedingungen für hybride serielle Roboter mit planaren Parallelmechanismen"
# SA Shan: Shan2019_S828, "Reduktion der Modellkomplexität von seriell-hybriden Palettierrobotern"
#
#
# Initialisierung
# Import (Library)
interface(warnlevel=0): 		# Unterdrücke die folgende Warnung.
restart: 				# Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: 		# Für Speicherung
;
with(StringTools): 			# Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
#with(ListTools):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;

# Import (hybriddyn)
read "../helper/proc_MatlabExport": 	# Exportiert einen Term als Matlab-Code in optimierter Form
read "../transformation/proc_rotx":  	# rotation um X
read "../transformation/proc_roty":  	# rotation um Y
read "../transformation/proc_rotz":  	# rotation um Z
;
read "../helper/proc_convert_s_t":   	# zeitabhängige Variable für Ableitung
read "../helper/proc_convert_t_s":   	# konstante Variable
read "../robot_codegen_constraints/proc_subs_kintmp_exp": #substitute für Spalte 1 mit Spalte 2
;
with(RealDomain): 			# Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":			#aktuelle Roboter, MDH-Tabelle
;
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions",robot_name):	#von Fourbar
;

# Variable Initialisierung
# Variable mit Winkeln der Nebenstruktur nur in Abhängigkeit der verallgemeinerten Koordinaten
# Ergebnisse von Trigonometrischer Elimination lesen
read sprintf("../codeexport/palh3m1TE/tmp/kinematic_constraints_maple_inert.m"):	# von aktuellem Roboter
;
kin_constraints_exist := kin_constraints_exist:
kintmp_qs := kintmp_qs:
kintmp_qt := kintmp_qt:
kintmp_subsexp := kintmp_subsexp:
# Variable entfernen
kintmp_subsexp:= Matrix(2*RowDimension(kintmp_s),2):
# Export
kintmp_qt := convert_s_t(kintmp_qs):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name):
#printf("Ausdrücke für kintmp_subsexp gespeichert (Maple). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
for i from 1 to RowDimension(kintmp_s) do
  tmp := kintmp_qs(i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d",robot_name, i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d.m", robot_name, i):
end do:
save kin_constraints_exist, kintmp_qs, kintmp_qt,kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert" ,robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
save kintmp_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_qs_maple_inert", robot_name):
#printf("Ausdrücke mit Inert-Arctan exportiert (Matlab). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_qs ));
#kc_symbols :=Transpose(kc_symbols);
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m",robot_name),2);

#printf("Fertig. %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
kintmp_qs(1,1):
kc_symbols(1,1):
kintmp_qs(2,1):


text:
#!/bin/bash maple
# Use `maple -q expand6.mpl -D n=100` to run

e := a0 + a1:
f := 0:

for i from 2 to (n - 1)
do
    f := f + cat(a, convert(i, string)):
    e := e + cat(a, convert(i, string)):
end do:

f := -f:
e := e ^ 2:

st := time[real]():
e := expand(e):
e := subs(a0 = f, e):
e := expand(e):
1000*(time[real]()-st);

print(e);

done

max_issues_repo_name: HQSquantumsimulations/symengine

max_forks_repo_name: HQSquantumsimulations/symengine


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

cases := [
    [
        [[a * b, b], [a]],
        [1 / a]
    ],
    [
        [[a + b, a * b], [b]],
        []
    ],
    [
        [[a * b + c, b], [a, c]],
        []
    ],
    [
        [[a * b, c * b, d * b], [a, c, d]],
        [a / c, a / d]
    ],
    [
        [[a + b + c, a^2 + b^2 + c^2, a^3 + b^3 + c^3, c], [a, b]],
        [a + b, a * b]
    ],
    [
        [[a + b, a * b, x, a * x + b * c], [a, b, c]],
        [a + b, a * b]
    ]
]:

num_passed := 0:
num_failed := 0:

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

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

max_stars_repo_name: iliailmer/AllIdentifiableFunctions

max_issues_repo_name: iliailmer/AllIdentifiableFunctions

max_forks_repo_name: iliailmer/AllIdentifiableFunctions


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

sys := [
diff(x5(t), t) = k5*x6(t) + k4*x6(t) - k6*x5(t)*x3(t),
diff(x6(t), t) = -k5*x6(t) - k4*x6(t) + k6*x5(t)*x3(t),
diff(x4(t), t) = -k3*x4(t) - k2*x4(t) + k1*x1(t)*x2(t),
diff(x2(t), t) = k3*x4(t) + k2*x4(t) + k1*x1(t)*x2(t),
diff(x1(t), t) = k4*x6(t) + k2*x4(t) - k1*x1(t)*x2(t),
diff(x3(t), t) = k5*x6(t) + k3*x4(t) - k6*x5(t)*x3(t),
y1(t) = x3(t),
y2(t) = x2(t)
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


text:
# This is from the original version of Gemma Halliwell's thesis.
# It is probably not useful.

`is_element/marked_digraphs` := (V::set) -> (m,n::nonnegint) -> proc(EIO)
 local E,VI,VO,NI,NO,v;

 if not (type(EIO,list) and nops(EIO) = 3) then
  return false;
 fi;

 E,VI,VO := op(EIO);

 if not `is_element/digraphs`(V)(E) then
  return false;
 fi;

 if not(type(VI,set) and nops(VI) = m and VI minus V = {}) then
  return false;
 fi;

 if not(type(VO,set) and nops(VO) = n and VO minus V = {}) then
  return false;
 fi;

 if VI intersect VO <> {} then return false; fi;

 NI := `in_neighbour_table/digraphs`(V)(E);
 NO := `out_neighbour_table/digraphs`(V)(E);

 for v in VI do
  if nops(NI[v]) <> 0 then return false; fi;
 od:

 for v in VO do
  if nops(NO[v]) <> 0 then return false; fi;
  if nops(NI[v]) <> 1 then return false; fi;
 od:

 return true;
end:

`internal_vertices/marked_digraphs` := (V) -> (m,n) -> proc(EIO)
 local E,VI,VO;
 E,VI,VO := op(EIO);
 return V minus (VI union VO);
end:


text:
(*
    以代表元为核心的对象
    主要以con划分对象，isol和tsol仅做参考作用
    保证RepSol对象的所有属性的相关变量都是global的。
*)
$ifndef _REPSOL_
$define _REPSOL_

RepSol:=module()
    option object;
    local
            # 局部函数
            getDisplayDcon::static, # 合并显示成立条件
            apList::static,         # 拓展list
            rmlist::static;         # 去除一层嵌套list
    export
            # 导出变量
            rep,                    # 代表元
            dcon:=[],               # 不变量方程
            acon:=[],               # 附加约束
            isol:=[],               # 代表元的通解
            tsol:=[],               # 代表元通解和特解的转化
            osol:=[],               # 对应的InvSol对象
            sid:=1,                 # 选择的最简条件
            # 导出函数
            ## 设置相关
            appendSol::static,      # 扩充一个RepSol对象成立的条件
            getCon::static,         # 获取一个RepSol对象成立的条件
            sortCon::static,        # 对成立条件进行排序
            selectCon::static,      # 选择最简成立条件
            rmCon::static,          # 删除某个成立条件的一部分
            ## 输出相关
            printRep::static,       # 简要显示代表元和所有可能的成立条件
            fullPrintRep::static,   # 显示代表元和完整的成立条件以及对应的不变量方程和变换方程的解
            ModulePrint::static,    # 简要显示代表元
            uniqueAndSort::static;      # 删去重复条件

    # 用于拓展一个代表元对象所能代表的区域
    appendSol:=proc(r::RepSol,s::InvSol)
        local ieq,sieq,isol,icon,tsols,tcons,i,n;
        if not assigned(r:-rep) then
            r:-rep:=getRep(s);
        end if;
        ieq:=s:-ieq;
        sieq:={ieq[]};
        isol:=s:-isol;
        icon:=s:-icon;
        tsols:=rmlist(s:-tsol);
        tcons:=rmlist(s:-tcon);
        n:=numelems(tsols);
        for i from 1 to n do
            apList(r:-dcon,ieq);
            apList(r:-acon,classifySolve({icon[],tcons[i][]}) minus sieq);
            apList(r:-isol,isol);
            apList(r:-tsol,tsols[i]);
            apList(r:-osol,s);
        end do;
        return;
    end proc:

    # 拓展list
    apList:=proc(lst::evaln)
        lst:=[eval(lst)[],_rest];
        return;
    end proc:

    # 获取一个RepSol对象成立的条件
    getCon:=proc(r::RepSol)
        return zip((x,y)->[getDisplayDcon(x)[],y[]],r:-dcon,r:-acon);
    end proc:

    # 不变量方程的简化显示
    getDisplayDcon:=proc(dcon)
        local n;
        n:=numelems(dcon);
        return [seq(Delta[i]=rhs(dcon[i]),i=1..n)];
    end proc:

    # 去掉一层嵌套list
    rmlist:=proc(x)
        return map(y->y[],x);
    end proc:

    # 对成立条件进行排序
    sortCon:=proc(r::RepSol)
        local con,ind;
        con:=getCon(r);
        ind:=sortByComplexity(con,index);
        r:-dcon:=r:-dcon[ind];
        r:-acon:=r:-acon[ind];
        r:-isol:=r:-isol[ind];
        r:-tsol:=r:-tsol[ind];
        r:-osol:=r:-osol[ind];
    end proc:

    # 选择最简成立条件
    selectCon:=proc(r::RepSol,sid)
        r:-sid:=sid;
        return;
    end proc:

    # 删除某个成立条件的一部分
    rmCon:=proc(r::RepSol,id::posint,con::set)
        r:-acon[id]:=r:-acon[id] minus con;
    end proc:

    # 简要显示代表元和所有可能的成立条件
    printRep:=proc(r::RepSol)
        print(r:-rep);
        print~(getCon(r));
        return ;
    end proc:

    # 显示代表元和完整的成立条件以及对应的不变量方程和变换方程的解
    fullPrintRep:=proc(r::RepSol)
        local i,n,con;
        print(r:-rep);
        con:=getCon(r);
        n:=numelems(con);
        for i from 1 to n do
            printf("[%d]------------------------------------------------------",i);
            print(con[i]);
            print(r:-isol[i]);
            print(r:-tsol[i]);
        end do;
        return;
    end proc:

    # 删去重复条件
    uniqueAndSort:=proc(r::RepSol)
        local id,con,fun,ind;
        fun:=(x,y,z)->[x,y,z];
        con:=getCon(r);
        id:=fun~(con,r:-isol,r:-tsol);
        ind:=uniqueObj(id,key=(x->x),'index');
        r:-dcon:=r:-dcon[ind];
        r:-acon:=r:-acon[ind];
        r:-isol:=r:-isol[ind];
        r:-tsol:=r:-tsol[ind];
        r:-osol:=r:-osol[ind];
        sortCon(r);
        return r;
    end proc:

    # 简要显示代表元
    ModulePrint:=proc(r::RepSol)
        return r:-rep;
    end proc:
end module:

$endif


text:
read "common-procedures.mpl":

#---------------------------------------------------------------------
# hi_mi_lo takes an arbitrary precision number x and returns three doubles such that:
# x ~ x_hi + x_mi + x_lo
hi_mi_lo:= proc(x)
local x_hi, x_mi, x_lo, res:
x_hi:= nearest(evalf(x)):
res:=x-x_hi:
if (res = 0) then
  x_mi:=0:
  x_lo:=0:
else
  x_mi:=nearest(evalf(res)):
  res:=x-(x_hi+x_mi);
  if (res = 0) then
    x_lo:=0:
  else
    x_lo:=nearest(evalf(res)):
  end if:
end if:
x_hi,x_mi,x_lo:
end:

#---------------------------------------------------------------------
# same as hi_mi_lo, but returns hexadecimal strings
ieeehexa3:=proc(x)
local x_hi, x_mi, x_lo;
  x_hi, x_mi, x_lo := hi_mi_lo(x):
  ieeehexa(x_hi), ieeehexa(x_mi), ieeehexa(x_lo):
end proc:

#---------------------------------------------------------------------
#Like poly_exact, but the first n coefficients as triple doubles,
#the next m coefiicients as double doubles and the rest as doubles
#
poly_exact32:=proc(P,n,m)
local deg,i, coef, coef_hi, coef_mi, coef_lo, Q:
Q:= 0:
convert(Q, polynom):
deg:=degree(P,x):
  for i from 0 to deg do
    coef :=coeff(P,x,i):
    coef_hi, coef_mi, coef_lo:=hi_mi_lo(coef):
    Q:= Q + coef_hi*x^i:
    if(i<m+n) then
        Q := Q + coef_mi*x^i:
    fi:
    if(i<n) then
	Q := Q + coef_lo*x^i:
    fi:
  od:
  return(Q):
end:

#---------------------------------------------------------------------
# Returns nearest double-double:
nearestDD := proc(x)
  local xh,xm:
  xh := nearest(x);
  xm := nearest(x - xh);
  evalf(xh + xm);
end:

#---------------------------------------------------------------------
# Returns nearest triple-double:
nearestTD := proc(x)
  local xh,xm,xl:
  xh := nearest(x);
  xm := nearest(x - xh);
  xl := nearest(x - (xh + xm));
  evalf(xh + xm + xl);
end:

max_stars_repo_name: squarePenguin/parvsl

max_issues_repo_name: squarePenguin/parvsl

max_forks_repo_name: squarePenguin/parvsl


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

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

(* type: mgga_exc *)

(* replace: "mbrxc_x\(" -> "xc_mgga_x_mbrxc_get_x(" *)

mbrxc_a1 := 0.074746:
mbrxc_a2 := 0.147:
mbrxc_a3 := 0.0032:

(* This is the derivative of f = (1+x)^(5/3)*exp(-2/3*x)/(x - 3) = (32*Pi)^(2/3)/(6*Q) *)
mbrxc_aux_dfdx := x -> -2/3 * (1 + x)^(2/3) * exp(-2*x/3) * (x^2 - 3*x + 6) / (x - 3)^2:

`diff/mbrxc_x` := proc(Q, g)
  - (32*Pi)^(2/3)/6 * diff(Q, g)/(Q^2 * mbrxc_aux_dfdx(mbrxc_x(Q)))
end proc:

mbrxc_Q := (x, t) ->
        mbrxc_a1*(2*t) - K_FACTOR_C + mbrxc_a2*x^2 + mbrxc_a3*x^4:

mbrxc_min_Q := 5.0e-13:
mbrxc_cQ := Q -> my_piecewise3(abs(Q) < mbrxc_min_Q,
  my_piecewise3(Q > 0, mbrxc_min_Q, -mbrxc_min_Q), Q):

mbrxc_v := x ->
  - (32*Pi)^(1/3)/(8*X_FACTOR_C) * exp(x/3)*(8 - exp(-x)*(x^2 + 5*x + 8))/(x*(1 + x)^(1/3)):

mbrxc_f := (x, u, t) ->
  - mbrxc_v(mbrxc_x(mbrxc_cQ(mbrxc_Q(x, t))))/2:

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


text:
%def genstub
%self.each do |msg|
%next unless msg.id
void svr_[%msg.name%]([%msg.name%] param, ccf::session_responder response,
		ccf::shared_session from, ccf::auto_zone& z);

%end

void dispatch(ccf::shared_session from,
		ccf::method_t method, ccf::msgobj param,
		ccf::session_responder response, ccf::auto_zone& z)
{
	switch(method) {
	%self.each do |msg|
	%next unless msg.id
	case [%msg.name%]::method:
		svr_[%msg.name%](param.as<[%msg.name%]>(), response, from, z);
		break;
	%end
	default:
		throw std::runtime_error("unknown method");  // FIXME
	}
}
%end  # genstub
%# vim: syntax=mplex

max_stars_repo_name: frsyuki/ccf

max_issues_repo_name: frsyuki/ccf

max_forks_repo_name: frsyuki/ccf


text:
# 基于 CodeParser.jar 的代码读取包。
# 优点在于：
# 	1. 可以去除预处理器必须在行首的限制。
# 	2. include默认使用相对路径。
#	3. 可以指定文件编码。
#	4. 默认每个文件只会include一遍。
# 注：会自动在 libname 包含的目录下寻找jar文件。
JavaCodeReader:=module()
	option  package;
	export 	ReadCode,ParseCode;
	local  	parseCode,
			parseStr,
			ModuleLoad;

	ModuleLoad:=proc()
		local jarName,jarPath,path,isFinded:=false;
		jarName:="CodeParser.jar";
		for path in :-libname do
			jarPath:=FileTools:-JoinPath([path,jarName]);
			if FileTools:-Exists(jarPath) then
				isFinded:=true;
				break;
			end if;
		end do;
		if not isFinded then
			error "找不到 %1 ，请确保 %2 下存在该文件。",jarName,[libname];
		end if;
		parseCode:=define_external(
			'mapleCall',
			JAVA,CLASSPATH=jarPath,
			CLASS="org.yjt.maple.ParseCodeUtil",
			'fin'::string,
			'fout'::string,
			'inputEncode'::string,
			'outputEncode'::string);
	end proc:
	ModuleLoad();

	# 采用Base64编码传递文件路径以支持中文
	# 注意Maple内部字符串的编码为UTF8
	parseStr:=proc(s::string)
		return StringTools:-Encode(s,'encoding'='base64');
	end proc:

	ParseCode:=proc(fin::string,fout::string,inEncode:="UTF8",outEncode:="UTF8")
		parseCode(parseStr(fin),parseStr(fout),inEncode,outEncode);
	end proc:

	ReadCode:=proc(fname::string,inEncode:="UTF8",outEncode:="UTF8")
		description "similar to `read`";
		local tmpFile;
		tmpFile:=cat(fname,".tmp.mpl");
		ParseCode(fname,tmpFile,inEncode,outEncode);
		read(tmpFile);
		FileTools[Remove](tmpFile);
		return;
	end proc:
end module:


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

`eta/one_cubes_prime` := (A::set) -> `eta/cubes`(1)(A);

`gamma/one_cubes_prime` :=
 (A::set,B::set) -> (p) -> (U,V) -> `gamma/cubes`(1)(A,B)(p)(U,V);


text:
`cartesian_product/sets` := proc()
 local A,P,a,p;
 if nargs = 0 then
  return {[]};
 else
  A := args[1];
  if not type(A,set) then
   error "Argument is not a set";
  fi;
  P := `cartesian_product/sets`(args[2..-1]);
  return {seq(seq([a,op(p)],p in P),a in A)};
 fi;
end:

`cartesian_product/lists` := proc()
 local A,P,a,p;
 if nargs = 0 then
  return [[]];
 else
  A := args[1];
  if not type(A,list) then
   error "Argument is not a list";
  fi;
  P := `cartesian_product/lists`(args[2..-1]);
  return [seq(seq([a,op(p)],p in P),a in A)];
 fi;
end:


text:
# The full Eilenberg-Zilber operad EZ has EZ(A) equal to the set of
# natural maps N_*(X) \to N_*(X)^{\otimes A} for simplicial sets X.
# By consideration of universal examples, we see that any such map
# sends N_{\leq n}(X) into (N_{\leq n}(X))^{\otimes A} for all n.
# Thus, if we define EZ_{\leq n}(A) to be the set of natural maps
# N_{\leq n}(X) \to (N_{\leq n}(X))^{\otimes A}, we get a quotient
# operad of EZ.  The group EZ(A) is an infinite product of copies of
# the integers and so is not computationally tractable.  We
# therefore work with EZ_{\leq n} instead.

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

 if x = 0 then return true; fi;

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

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

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

 if nops(x) = 0 or nops(x) > n+1 then return false; fi;

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

 if map(op,{op(x)}) <> A then
  return false;
 fi;

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

 return true;
end:

`is_equal/eilenberg_zilber` := (n::nonnegint) -> (A::set) -> (x,y) -> evalb(x = y):

`random_basis_element/eilenberg_zilber` := (n::nonnegint) -> (A) -> proc(m_)
 local m,t,M,a,i;

 if A = {} then
  if nargs = 0 or m_ = 0 then
   return Z({});
  else
   error("A is empty");
  fi;
 fi;

 m := `if`(nargs > 0,min(n,m_),rand(0..n)());

 t := table():
 M := {seq(i,i=0..m)};

 for a in A do t[a] := random_nonempty_subset_of(M); od;

 for i from 0 to m-1 do
  a := random_element_of(A);
  t[a] := t[a] union {i,i+1};
 od:

 return Z(seq(select(a->member(i,t[a]),A),i=0..m));
end:

`random_element/eilenberg_zilber` := (n::nonnegint) -> (A) -> proc(coeff_range_,num_terms_)
 local coeff_range,num_terms;

 coeff_range := -3..3;
 num_terms := 5;
 if nargs > 0 then coeff_range := args[1]; fi;
 if nargs > 1 then num_terms := args[2]; fi;

 add(rand(coeff_range)() * `random_basis_element/eilenberg_zilber`(n)(A)(),i=1..num_terms);
end:

`transpose0/eilenberg_zilber` := (A::set) -> proc(x)
 local m,M,t,a;

 m := nops(x) - 1;
 M := [seq(i,i=0..m)];
 t := table():
 for a in A do
  t[a] := select(i -> member(a,op(i+1,x)),M);
 od;

 return Zt(eval(t));
end;

`transpose/eilenberg_zilber` := (A::set) ->
 apply_linear(`transpose0/eilenberg_zilber`(A));

`detranspose0/eilenberg_zilber` := (A::set) -> proc(x)
 local t,m;

 t := op(1,x);
 m := max(op(map(op,{{0},seq(t[a],a in A)})));
 return Z(seq(select(a -> member(i,t[a]),A),i=0..m));
end:

`detranspose/eilenberg_zilber` := (A::set) ->
 apply_linear(`detranspose0/eilenberg_zilber`(A));

`eta/eilenberg_zilber` := (n::nonnegint) -> (A::set) ->
 `if`(nops(A)=1,add(Z(A $ k),k=1..(n+1)),FAIL);

`source_degree/eilenberg_zilber` := (A::set) -> proc(x)
 return nops(x) - 1;
end:

`target_degrees/eilenberg_zilber` := (A::set) -> proc(x)
 local xt,d,a;

 xt := op(1,`transpose/eilenberg_zilber`(A)(x));
 d := table();
 for a in A do
  d[a] := nops(xt[a]) - 1;
 od:

 return eval(d);
end:

`gamma0/eilenberg_zilber` := (A::set,B::set) -> (p) -> proc(y,x)
 local F,s,d,t,u,a,b;

 F := fibres(A,B)(p);
 s := op(1,`transpose/eilenberg_zilber`(B)(y));
 d := `target_degrees/eilenberg_zilber`(B)(y);

 t := table():
 for b in B do
  if `source_degree/eilenberg_zilber`(F[b])(x[b]) <> d[b] then
   return 0;
  fi;

  t[b] := op(1,`transpose/eilenberg_zilber`(F[b])(x[b]));
 od;

 u := table();
 for a in A do
  b := p[a];
  u[a] := [seq(s[b][j+1], j in t[b][a])];
 od;

 return `detranspose/eilenberg_zilber`(A)(Zt(u));
end:

`gamma/eilenberg_zilber` := (n::nonnegint) -> (A::set,B::set) -> (p) -> proc(y,x)
 local c,y1,F,d0,BL,M,t,z,b,m,mt;

 if y = 0 then
  return 0;
 elif type(y,`+`) then
  return map(t -> `gamma/eilenberg_zilber`(n)(A,B)(p)(t,x),y);
 fi;

 if type(y,`*`) then
  c,y1 := selectremove(type,y,integer);
 else
  c := 1;
  y1 := y;
 fi;

 if not(type(y1,specfunc(Z))) then
  return FAIL;
 fi;

 F := fibres(A,B)(p);

 d0 := `target_degrees/eilenberg_zilber`(B)(y1);

 BL := sort([op(B)]);
 M := [[1]];
 for b in BL do
  t := map(coeff_split,sum_terms(x[b]));
  t := select(u -> `source_degree/eilenberg_zilber`(F[b])(u[2]) = d0[b],t);
  M := [seq(seq([m[1] * u[1],seq(m[i],i=2..nops(m)),u[2]],u in t),m in M)];
 od:

 z := 0;

 for m in M do
  mt := table([seq(BL[i]=m[i+1],i=1..nops(BL))]);
  z := z + c * m[1] * `gamma0/eilenberg_zilber`(A,B)(p)(y1,mt);
 od;

 return z;
end:


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

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

(* type: work_mgga_c *)

$define lda_c_pw_params
$include "lda_c_pw.mpl"

$define lda_x_params
$include "lda_x.mpl"

gamma_x  := 0.11:
gamma_ss := 1.6:
gamma_ab := 0.14:

params_a_c_x  := [0.8085,  0.6682, 0.1420]:
params_a_c_ss := [0.2606, -0.9608, 0.9023]:
params_a_c_ab := [1.2033, -2.2717, 0.9596]:

b98_q := (x, t, u) -> (u/4 - t + x^2/8 + K_FACTOR_C)/K_FACTOR_C:

b98_g := (gamma, cc, q) -> add(cc[i]*(gamma*q^2/sqrt(1 + gamma*q^2))^(i-1), i=1..3):

f_b98 := (lda_func, g_ss, cc_ss, g_ab, cc_ab, rs, z, xs0, xs1, ts0, ts1, us0, us1) ->
  + lda_stoll_par(lda_func, rs,  z,  1)
    * b98_g(g_ss, cc_ss, b98_q(xs0, ts0, us0)) * Fermi_D(xs0, ts0)
  + lda_stoll_par(lda_func, rs, -z, -1)
    * b98_g(g_ss, cc_ss, b98_q(xs1, ts1, us1)) * Fermi_D(xs1, ts1)
  + lda_stoll_perp(lda_func, rs, z)
    * b98_g(g_ab, cc_ab, (b98_q(xs0, ts0, us0) + b98_q(xs1, ts1, us1))/2):

f := (rs, z, xt, xs0, xs1, ts0, ts1, us0, us1) ->
  + f_lda_x(rs,  1)*b98_g(gamma_x, params_a_c_x, b98_q(xs0, ts0, us0))
  + f_lda_x(rs, -1)*b98_g(gamma_x, params_a_c_x, b98_q(xs1, ts1, us1))
  + f_b98(f_pw, gamma_ss, params_a_c_ss, gamma_ab, params_a_c_ab,
        rs, z, xs0, xs1, ts0, ts1, us0, us1):


text:
#################################
# Using direct einstein's tensor
# equation to built IREs
################################
read("eq_ire.mpl"):

res2:=grcomponent(HSR):
tg2:=grcomponent(TG2(up),[x]):

with(DEtools):
addcoords(compact,[t,x],[t,ctfm(x)]);

ire_fields := [res2, tg2, restt, restx, resxx, resthth, ire_rpsi_direct]:
c_ire_fields := [0, 0, 0 , 0 , 0 , 0 , 0]:

# Compactifying Einstein's equation:
for ii from 1 to nops(ire_fields) do
 c_ire_fields[ii] := PDEchangecoords(ire_fields[ii],[t,x],compact,[t,x]);
 c_ire_fields[ii] := subs({diff(ctfm(x),x)=ctfmp(x),diff(ctfm(x),x,x)=ctfmpp(x)},c_ire_fields[ii]);

 for kk from 1 to nops(grid_functions) do
        fn := grid_functions[kk];
        c_ire_fields[ii] := subs(fn(t,ctfm(x))=fn(t,x),c_ire_fields[ii]);
 end do:

end do:

printf("checking compactification... all residuals should be zero:\n");
for ii from 1 to nops(ire_fields) do
   res := simplify(ire_fields[ii] - eval(c_ire_fields[ii],{ctfm(x) = x, ctfmp(x) = 1, ctfmpp(x) = 0})):
   printf("res = %a\n",res);
end do:

for ii from 1 to nops(ire_fields) do
   c_ire_fields[ii] := collect(collect(eval(expand(c_ire_fields[ii]),{ctfmpp(x)=octfmp(x)*ctfmp(x)^2}),1/ctfmp(x)),1/ctfm(x));
end do:

res2 := c_ire_fields[1]:
tg2 := c_ire_fields[2]:
restt := c_ire_fields[3]:
restx := c_ire_fields[4]:
resxx := c_ire_fields[5]:
resthth := c_ire_fields[6]:
ire_rpsi_direct := c_ire_fields[7]:

Gen_Eval_Code(restx,input="c",proc_name="ire_restx");
Gen_Eval_Code(resxx,input="c",proc_name="ire_resxx");
Gen_Eval_Code(restt,input="c",proc_name="ire_restt");
Gen_Eval_Code(resthth,input="c",proc_name="ire_resthth");
Gen_Eval_Code(tg2-Lamx(t,x),input="c",proc_name="ire_val_lamx");
Gen_Eval_Code(ire_rpsi_direct,input="c",proc_name="ire_rpsi_direct");

pl:=table([ t=[-1,-1],x=[-1,0],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

Gen_Eval_Code(res2,input="c",proc_name="ire_hs");


text:
impulse process ACV_remote_deleteAllPortals {
/summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Blue"],Duration:2147483647}
/summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Orange"],Duration:2147483647}
/summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Purple"],Duration:2147483647}
/summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Red"],Duration:2147483647}
ACV_remote_deletePortals()
}


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

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

(* type: gga_exc *)

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

ityh_enhancement := xs -> b88_f(xs):
ityh_attenuation := a  -> attenuation_erf(a):

ityh_k_GGA := (rs, z, xs) -> sqrt(9*Pi/(2*X_FACTOR_C*ityh_enhancement(xs))) * n_spin(rs, z)^(1/3):
ityh_f_aa  := (rs, z, xs) -> ityh_attenuation(p_a_hyb_omega_0_/(2*ityh_k_GGA(rs, z, xs))):

ityh_f := (rs, z, xs) -> ityh_f_aa(rs, z, xs) * ityh_enhancement(xs):

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


text:
mobius_complex := proc()
 local T,V,E,F,R,r,v,a,mob,mob_arc,i,e;

 V := [seq(seq([i,a],a=0..2),i=0..1)];
 F := [
  seq([[0,a],[0,modp(a+1,3)],[1,a]],a=0..2),
  seq([[0,a],[0,modp(a+1,3)],[1,modp(a+1,3)]],a=0..2),
  seq([[0,a],[1,modp(a+1,3)],[1,modp(a+2,3)]],a=0..2)
 ];

 E := map(f -> ([f[1],f[2]],[f[1],f[3]],[f[2],f[3]]),F);

 T := table([]);

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

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

 v := table():
 for a from 0 to 2 do
  v[0,a] := [a/6,0];
  v[1,a] := [a/3,1];
 od:

 R := 1;
 r := 1/3:
 mob := mobius_embedding(R,r);

 mob_arc := mobius_arc(R,r);

 T["plot"] :=
  display(
   spacecurve(mob(t,1),t=0..1,colour=blue),
   plot3d(mob(t,u),t=0..1,u=0..1,style=wireframe,colour="LightGray"),
   seq(mob_arc(op(v[op(e[1])]),op(v[op(e[2])]),colour=red),e in E),
   axes=none,scaling=constrained
  );

 return eval(T);
end():

flat_mobius_complex := proc(n::posint)
 local f,g,T,V,E,F,r,v,a,i;

 f := (i) -> modp(i+1,4*n+2);
 g := (i) -> modp(i+2*n+3/2+(-1)^i/2,4*n+2);

 V := [seq(i,i=0..4*n+1)];
 E := [seq([i,f(i)],i=0..4*n+1),seq([i,g(i)],i=0..4*n+1)];
 F := [seq([i,f(i),g(i)],i=0..4*n+1)];

 T := table([]);

 T["vertices"] := V;
 T["edges"] := E;
 T["faces"] := F;
 T["max_simplices"] := T["faces"];
 T["all_simplices"] := [map(v -> [v],V),op(E),op(F)];

 T["embedding_dim"] := 3;
 T["embedding"] := table([seq(i = evalf(mobius_embedding(0.2)(i/(4*n+2),1)),i=0..4*n+1)]);

 T["plot"] := `plot/raw_simplicial_complex`(T["vertices"])(T["faces"],3,T["embedding"]);

 return eval(T);

end:


text:
# Proof by PDF that
# If R ~ rayleigh(2) then R^2 ~ stdChiSq(2)

PDF_rayleigh := (a, x) -> 2*x*exp(-x^2/a)/a;

# Transform on random variable is f(R)=R^2. Inverse function is g(x) = sqrt(x)
PDF_transform := (x) -> PDF_rayleigh(2,sqrt(x))*diff(sqrt(x),x);

PDF_chisq := (n, x) -> x^(1/2*n-1)*exp(-1/2*x)/(2^(1/2*n))/GAMMA(1/2*n);
simplify(PDF_chisq(2,x));

evalb(PDF_chisq(2,x) = PDF_transform(x));

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


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

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

(* type: gga_exc *)

$include "scaling.mpl"

(* parameters *)
wpbeh_A :=  1.0161144:
wpbeh_B := -0.37170836:
wpbeh_C := -0.077215461:
wpbeh_D :=  0.57786348:
wpbeh_E := -0.051955731:

(*
 Note that kF has a 6 and not a 3 as it should in principle
 be. This is because the HSE formula, if one would take the papers
 seriously, does not fulfill the spin sum-rule. This is probably
 an oversight from them. So, we have to choose, either a 6 or a 3.

 Nwchem seems to have the factor of 6, but VASP and espresso have
 a 3. This would amount to rescaling omega by a factor of
 cbrt(2). We follow the quantum chemistry community and put the 6.
*)

(* Cutoff criterion below which to use polynomial expansion *)
EGscut     := 0.08:
wcutoff    := 14:
expfcutoff := 700:

(* first let us calculate H(s) *)
wpbeh_Ha1 := 0.00979681:
wpbeh_Ha2 := 0.0410834:
wpbeh_Ha3 := 0.187440:
wpbeh_Ha4 := 0.00120824:
wpbeh_Ha5 := 0.0347188:
wpbeh_H := s ->
  + (wpbeh_Ha1*s^2 + wpbeh_Ha2*s^4)
  / (1 + wpbeh_Ha3*s^4 + wpbeh_Ha4*s^5 + wpbeh_Ha5*s^6):

(*
  Now we calculate F(s). We use the parameters that were in the original code,
  but these constants are:

  wpbeh_Fc1 := 4*wpbeh_A^2/(9*wpbeh_C) + (wpbeh_B - wpbeh_A*wpbeh_D)/wpbeh_C:
  wpbeh_Fc2 := -4/(3*36*wpbeh_C):
*)
wpbeh_Fc1 := 6.4753871:
wpbeh_Fc2 := 0.47965830:
wpbeh_F := s ->
  wpbeh_Fc1*wpbeh_H(s) + wpbeh_Fc2:

(* several auxiliary variables *)
eb1  := w -> my_piecewise3(w < wcutoff, 1.455915450052607, 2):
aux1 := s -> wpbeh_D + s^2*wpbeh_H(s):
aux2 := s -> 9*wpbeh_H(s)*s^2/(4*wpbeh_A):
aux3 := (w, s) -> aux1(s) + w^2:
aux4 := (w, s) -> s^2*wpbeh_H(s) + eb1(w)*w^2:
aux5 := (w, s) -> 9*aux4(w, s)/(4*wpbeh_A):
aux6 := (w, s) -> wpbeh_D + aux4(w, s):

(* and now G(s) *)
Ga := s ->
  + sqrt(Pi)*(
    + 15*wpbeh_E
    + 6*wpbeh_C*(1 + wpbeh_F(s)*s^2)*aux1(s)
    + 4*wpbeh_B*aux1(s)^2
    + 8*wpbeh_A*aux1(s)^3)
  / (16*aux1(s)^(7/2))
  - (3*Pi/4)*sqrt(wpbeh_A)*exp(aux2(s))*(1 - erf(sqrt(aux2(s)))):
Gb := s ->
  15*sqrt(Pi)*s^2/(16*aux1(s)^(7/2)):

wpbeh_EGa1 := -0.02628417880:
wpbeh_EGa2 := -0.07117647788:
wpbeh_EGa3 :=  0.08534541323:
wpbeh_EG := s -> my_piecewise3(s > EGscut,
  -(3*Pi/4 + Ga(s))/Gb(s),
  wpbeh_EGa1 + wpbeh_EGa2*s^2 + wpbeh_EGa3*s^4
):

term2 := s-> (
  + aux1(s)^2*wpbeh_B
  + aux1(s)*wpbeh_C
  + 2*wpbeh_E
  + aux1(s)*s^2*wpbeh_C*wpbeh_F(s)
  + 2*s^2*wpbeh_EG(s)
)/(2*aux1(s)^3):

term3 := (w, s) -> -w*(
  + 4*aux3(w, s)^2*wpbeh_B
  + 6*aux3(w, s)*wpbeh_C
  + 15*wpbeh_E
  + 6*aux3(w, s)*s^2*wpbeh_C*wpbeh_F(s)
  + 15*s^2*wpbeh_EG(s)
)/(8*aux1(s)*aux3(w, s)^(5/2)):

term4 := (w, s) -> -w^3*(
  + aux3(w, s)*wpbeh_C
  + 5*wpbeh_E
  + aux3(w, s)*s^2*wpbeh_C*wpbeh_F(s)
  + 5*s^2*wpbeh_EG(s)
)/(2*aux1(s)^2*aux3(w, s)^(5/2)):

term5 := (w, s) -> -w^5*(
  + wpbeh_E
  + s^2*wpbeh_EG(s)
)/(aux1(s)^3*aux3(w, s)^(5/2)):

t10 := (w, s) ->
  1/2*wpbeh_A*log(aux4(w, s)/aux6(w, s)):

(* Use simple gaussian approximation for large w *)
term1_largew := (w, s) ->
  -1/2*wpbeh_A*(-xc_E1_scaled(aux5(w, s)) + log(aux6(w, s)) - log(aux4(w, s))):

(* For everything else use the full blown expression *)
ea1 := -1.128223946706117:
ea2 :=  1.452736265762971:
ea3 := -1.243162299390327:
ea4 :=  0.971824836115601:
ea5 := -0.568861079687373:
ea6 :=  0.246880514820192:
ea7 := -0.065032363850763:
ea8 :=  0.008401793031216:

np1 := w ->
  - 1.5*ea1*sqrt(wpbeh_A)*w
  + 27*ea3*w^3/(8*sqrt(wpbeh_A))
  - 243*ea5*w^5/(32*(wpbeh_A)^(3/2))
  + 2187*ea7*w^7/(128*(wpbeh_A)^(5/2)):

np2 := w ->
  - wpbeh_A
  + 9*ea2*w^2/4.0
  - 81*ea4*w^4/(16*wpbeh_A)
  + 729*ea6*w^6/(64*wpbeh_A^2)
  - 6561*ea8*w^8/(256*wpbeh_A^3):

t1 := (w, s) ->
  1/2*(np1(w)*Pi*xc_erfcx(sqrt(aux5(w, s))) - np2(w)*xc_E1_scaled(aux5(w, s))):

f2 := (w, s) ->
  1/2*ea1*sqrt(Pi)*wpbeh_A/sqrt(aux6(w, s)):
f3 := (w, s) ->
  1/2*ea2*wpbeh_A/aux6(w, s):
f4 := (w, s) ->
  ea3*sqrt(Pi)*(-9/(8*sqrt(aux4(w, s))) + 0.25*wpbeh_A/aux6(w, s)^(3/2)):
f5 := (w, s) ->
  (ea4/128)*(-144/aux4(w, s) + 64*wpbeh_A/aux6(w, s)^2):
f6 := (w, s) ->
  ea5*(3*sqrt(Pi)*(3*aux6(w, s)^(5/2)*(9*aux4(w, s) - 2*wpbeh_A)
    + 4*aux4(w, s)^(3/2)*wpbeh_A^2))/(32*aux6(w, s)^(5/2)*aux4(w, s)^(3/2)*wpbeh_A):
f7 := (w, s) ->
  ea6*((32*wpbeh_A/aux6(w, s)^3 + (-36 + 81*s^2*wpbeh_H(s)/wpbeh_A)/aux4(w, s)^2))/32:
f8 := (w, s) ->
  ea7*(-3*sqrt(Pi)*(-40*aux4(w, s)^(5/2)*wpbeh_A^3
    + 9*aux6(w, s)^(7/2)*(27*aux4(w, s)^2 - 6*aux4(w, s)*wpbeh_A + 4*wpbeh_A^2)))
  /(128*aux6(w, s)^(7/2)*aux4(w, s)^(5/2)*wpbeh_A^2):
f9 := (w, s) -> (
  + 324*ea6*eb1(w)*aux6(w, s)^4*aux4(w, s)*wpbeh_A
  + ea8*(384*aux4(w, s)^3*wpbeh_A^3
    + aux6(w, s)^4*(-729*aux4(w, s)^2 + 324*aux4(w, s)*wpbeh_A - 288*wpbeh_A^2))
  )/(128*aux6(w, s)^4*aux4(w, s)^3*wpbeh_A^2):

t2t9 := (w, s) ->
  + f2(w, s)*w + f3(w, s)*w^2 + f4(w, s)*w^3 + f5(w, s)*w^4
  + f6(w, s)*w^5 + f7(w, s)*w^6 + f8(w, s)*w^7 + f9(w, s)*w^8:

term1 := (w, s) -> my_piecewise3(
  w > wcutoff, term1_largew(w, s),
  t1(m_min(w, wcutoff), s) + t2t9(m_min(w, wcutoff), s) + t10(m_min(w, wcutoff), s)
):

f_wpbeh0 := (w, s) -> - 8/9 *(
  term1(w, s) + term2(s) + term3(w, s) + term4(w, s) + term5(w, s)
):

f_wpbeh := (rs, z, x) ->
  f_wpbeh0(nu(rs, z), m_max(1e-15, s_scaling_2(X2S*x))):

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


text:
(*
 Copyright (C) 2020 Susi Lehtola

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

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

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

csk_p := x -> X2S^2*x^2:
csk_q := u -> X2S^2*u:

(* Equation (21) *)
csk_z  := (p, q) -> 20/9*q - 40/27*p:

(* Equation (22) *)
csk_f0 := (p, q, z) ->  1 + 5*p/3 + z*csk_I(z):

(*
   The function I(z) contains exp(-1/|z|^a) which is numerically
   challenging for small |z|. Because of this, we truncate the term
   close to z=0 for negative z, since I(z) -> 1 for small |z|.
   For positive z, I(z)=1 identically, because of the step function.
*)
csk_I_negz := z -> (1 - exp(-1/abs(z)^params_a_csk_a))^(1/params_a_csk_a):
csk_I_cutoff_small := (-log(DBL_EPSILON))^(-1/params_a_csk_a):
csk_I_cutoff_large := (-log(1 - DBL_EPSILON))^(-1/params_a_csk_a):

csk_I := z -> my_piecewise5(
      z < -csk_I_cutoff_large, 0,
      z > -csk_I_cutoff_small, 1,
      csk_I_negz(m_max(m_min(z, -csk_I_cutoff_small), -csk_I_cutoff_large))
  ):

csk_f := (x, u) ->
  csk_f0(csk_p(x), csk_q(u), csk_z(csk_p(x), csk_q(u))):

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


text:
#@ Not autoload

# This was originally written for Mathematica, and needs more
# translation.

#**************************************************************
# the rest of this does not appear in Deligne's formulaire.

# The cubical structure

#h[{{x0_,y0_,z0_},
#   {x1_,y1_,z1_},
#   {x2_,y2_,z2_}}] :=
#     Det[{{x0,y0,z0},{x1,y1,z1},{x2,y2,z2}}] z0 z1 z2 /
#     ( Det[{{x0,z0},{x1,z1}}] *
#       Det[{{x1,z1},{x2,z2}}] *
#       Det[{{x2,z2},{x0,z0}}] )

# z in terms of x where y = 1, valid mod x^10

#zx =
#  x^3 - a1 x^4 + a1^2 x^5 + a2 x^5 - a1^3 x^6 - 2 a1 a2 x^6 - a3 x^6 +
#   a1^4 x^7 + 3 a1^2 a2 x^7 + a2^2 x^7 + 3 a1 a3 x^7 + a4 x^7 - a1^5 x^8 -
#   4 a1^3 a2 x^8 - 3 a1 a2^2 x^8 - 6 a1^2 a3 x^8 - 3 a2 a3 x^8 -
#   3 a1 a4 x^8 + a1^6 x^9 + 5 a1^4 a2 x^9 + 6 a1^2 a2^2 x^9 + a2^3 x^9 +
#   10 a1^3 a3 x^9 + 12 a1 a2 a3 x^9 + 2 a3^2 x^9 + 6 a1^2 a4 x^9 +
#   3 a2 a4 x^9 + a6 x^9

# Invariant differential mod x^10, where y = 1

#omx =
#  dx (-1 + a1 x - a1^2 x^2 - a2 x^2 + a1^3 x^3 + 2 a1 a2 x^3 + 2 a3 x^3 -
#     a1^4 x^4 - 3 a1^2 a2 x^4 - a2^2 x^4 - 6 a1 a3 x^4 - 2 a4 x^4 +
#     a1^5 x^5 + 4 a1^3 a2 x^5 + 3 a1 a2^2 x^5 + 12 a1^2 a3 x^5 +
#     6 a2 a3 x^5 + 6 a1 a4 x^5 - a1^6 x^6 - 5 a1^4 a2 x^6 - 6 a1^2 a2^2 x^6 -
#     a2^3 x^6 - 20 a1^3 a3 x^6 - 24 a1 a2 a3 x^6 - 6 a3^2 x^6 -
#     12 a1^2 a4 x^6 - 6 a2 a4 x^6 - 3 a6 x^6)

# The following functions are invertible near [0:1:0], and
# they satisfy uu ww = vv^3.  Note that vv depends only on y
# and z.  Moreover, {uu != 0} \cap {x = 0} = {0} as schemes.

#uu = y^2 + a1 x y - a2 x^2 + a3 y z - a4 x z - a6 z^2
#vv = y^2 + a3 y z - a6 z^2
#ww = a1^2 x^2 y^2 + a2 x^2 y^2 - a1 x y^3 + y^4 + a2 a3 x^2 y z -
#   2 a1 a4 x^2 y z + a1^2 a2 x y^2 z + a2^2 x y^2 z - a1 a3 x y^2 z +
#   a4 x y^2 z - a1^3 y^3 z - 2 a1 a2 y^3 z + 2 a3 y^3 z + a4^2 x^2 z^2 -
#   a2 a6 x^2 z^2 + a2^2 a3 x y z^2 - 2 a1 a2 a4 x y z^2 + a3 a4 x y z^2 +
#   a1 a6 x y z^2 - 2 a1 a2 a3 y^2 z^2 + a3^2 y^2 z^2 + 3 a1^2 a4 y^2 z^2 +
#   2 a2 a4 y^2 z^2 - 2 a6 y^2 z^2 + a2 a4^2 x z^3 - a2^2 a6 x z^3 -
#   a4 a6 x z^3 + 2 a2 a3 a4 y z^3 - 3 a1 a4^2 y z^3 + 2 a1 a2 a6 y z^3 -
#   2 a3 a6 y z^3 + a4^3 z^4 - 2 a2 a4 a6 z^4 + a6^2 z^4

# We now give some formulae that identify the singular locus.

#Frule = {F  -> HomogeneousCubic,
#         Fx -> D[HomogeneousCubic,x],
#         Fy -> D[HomogeneousCubic,y],
#         Fz -> D[HomogeneousCubic,z]}

#A0 = 4 F - 2 x Fx - (2 y + a3 z) Fy
#A1 = - 2 Fx + a1 Fy
#A2 = 4 a2 z F +
#     - ((a1 a3 + a4) z^2 + (a1^2 + 2 a2) x z + 3 x^2 + a1 y z) Fx +
#     - (a2 a3 z^2 + (a1 a2 + 2 a3) x z + 2 a1 x^2 + (2 a2 z + 4 x) y) Fy
#
#AA0 = 2 x^3 - B4 x z^2 - B6 z^3
#AA1 = 6 x^2 + B2 x z + B4 z^2
#AA2 = x^4 - B4 x^2 z^2 - 2 B6 x z^3 - B8 z^4

#Dl = (-7 B2 B4 z^2 + 27 B6 z^2 - 36 B4 x z + 3 B2 x^2) AA0 +
#     (-8 B4^2 z^3 + 2 B2 B6 z^3 + B2 B4 x z^2 - 9 B6 x z^2 +
#      12 B4 x^2 z - B2 x^3) AA1 +
#      B2^2 z AA2

# The singular point, where Delta = 0 and c4 is invertible.

#xs = 18 b6 - b2 b4

#ys = (a1 a4 - 2 a2 a3) b2 + 12 a3 b4 - 9 a1 b6

#zs = c4

# This evaluates to {0,0,0,0}, showing that [xs:ys:zs] is
# indeed a singular point.

CheckSingular0 =
({F,Fx,Fy,Fz} /. Frule /. {x->xs,y->ys,z->zs}) -
 {-c6,36,0,3 b2} Delta

# The following evaluates to {0,0}.  Because points with z=0
# are smooth and zs = c4 is assumed invertible, we conclude
# that every singular point is projectively equivalent to
# [xs:ys:zs].

CheckSingular1 =
({-2 b2 z Fx -
 (6 a1 x + 12 y - a1 b2 z + 18 a3 z) Fy +
 24 z Fz,
 a1 b2 z Fx +
 (3 a1^2 x + 6 a1 y + (2 a2 b2 - 3 a1 a3 - 24 a4)z ) Fy -
 12 a1 z Fz} /. Frule) -
{ (zs x - xs z) z^2 ,
  (zs y - ys z) z^2}

# Grobner basis for the singular points over the locus where
# c4 = Delta = 2 = 0 (with y > x > a3 > a1 > a6 > a4 > a2)

AdditiveLocusMod2 = {
 y^2 + x^3 + a2 x^2 + a4 x + a6,
 x^2 y + a4 y + a3 x^2 + a2 a3 x + a1 a6 + a3 a4,
 a3 y + x^3 + a4 x,
 a1 y + x^2 + a4,
 x^4 + a3^2 x + a1^2 a6 + a4^2 + a1 a3 a4 + a2 a3^2,
 a3 x^3 + a1^3 a6 + a1 a4^2 + a1^3 a2 a4,
 a1^2 a4^2 + a3^2 x^2,
 a3 + a1 x,
 a3^3 + a1^2 a3 a4,
 a1 a3^2 + a1^3 a4,
 a1^3 a3,
 a1^4}

ReducedAdditiveLocusMod2 = {y^2 + a2 a4 + a6, x^2 + a4, a1, a3}

AdditiveLocusMod3 = {
 y - a1 x - a3,
 a1^4 - a1^2 a2 + a2^2,
 -a3^2 - a6 - a1 a3 x + a4 x - x^3,
 a1 a3 - a4 + a1^2 x + a2 x,
 a1^3 a3 + a1 a2 a3 - a1^2 a4 - a2 a4,
 -(a1^2 a3^2) - a1 a3 a4 - a4^2,
 a1^2 a3^2 + a2 a3^2 + a1^2 a6 + a2 a6 - a1 a3 x^2 + a4 x^2}

ReducedAdditiveLocusMod3 = {
 y - a1 x - a3,
 x^3 + a6 + a3^2,
 a4 - a1 a3,
 a2 + a1^2}

ReducedAdditiveLocusRational = {
 y + a1 x/2 + a3/2,
 x^2 - b2 x/6 - b4/6,
 a6 + a3^2/4 - a2^3/27 - a1^2 a2^2/36 - a1^4 a2/144 - a1^6/1728,
 a4 + a1 a3/2 - a2^2/3 - a1^2 a2/6 - a1^4/48}

#**************************************************************
# Various quotients of localisations of Z[a1,..,a6] that give
# etale maps to the elliptic moduli stack (I think)

# etale1 covers D(6 c6)

etale1 = {a1 -> 0, a2 -> 0, a3 -> 0, a6 -> 1}

# etale2 covers D(6 c4)

etale2 = {a1 -> 0, a3 -> 0, a2 -> 0, a4 -> 1}

# etale3 covers a neighbourhood of V(2) D(c4)

etale3 = {a1 -> 1, a3 -> 0, a2 -> 0, a4 -> 0}

# etale4 covers a neighbourhood of V(3) D(c4)

etale4 = {a1 -> 0, a3 -> 0, a2 -> 1, a4 -> 0}

# etale5 covers a neighbourhood of V(2,a1) ; etale over D(3)

etale5 = {a3 -> 1, a2 -> 0, a4 -> 0, a6 -> 0}

# etale6 covers a neighbourhood of V(3,a2) ; etale over D(2)

etale6 = {a1 -> 0, a3 -> 0, a4 -> 1, a6 -> 0}

# Det[Jac[etale6,a2]] (for example) is the Jacobian
# determinant for the map
#
#  Z[a1,a2,a3,a4,a6] -> Z[v,r,s,t,a2]
#
# (where v = 1/u) given by alpha.  This must be invertible for
# the map to be etale.

Jac[er_,aa_] :=
 Outer[D,
  FixedPoint[Expand[# /. er /. u -> 1/v]&, alpha /@ {a1,a3,a2,a4,a6}],
  {v,r,s,t,aa}
 ]

#**************************************************************
# Singular fibres

# The curve y^2 z + x y z = x^3 is a nodal cubic, with
# multiplicative formal group.

multcurve = {a1->1, a3->0, a2->0, a4->0, a6->0}

# There is a birational map f from P^1 to the curve, with
# inverse g.	This sends both 0 and infinity to the singular
# point [0:0:1], and the restriction to G_m is a homomorphism.
# The discriminant is zero and the j invariant is infinite.

fm[s_,t_] := {s t (s - t), t^2 s, (s - t)^3}

gm[x_,y_,z_] := {y + x, y}

# The curve y^2 = x^3 is a cuspidal cubic, with additive
# formal group.

addcurve = {a1->0, a3->0, a2->0, a4->0, a6->0}

# There is a birational map f from P^1 to the curve, with
# inverse g.	This sends infinity to the singular point
# [0:0:1] with multiplicity two, and the restriction to G_a is
# a homomorphism the discriminant is zero and the j invariant
# is undefined.

fa[s_,t_] := {t^2 s, t^3, s^3}

ga[x_,y_,z_] := {x, y}

#**************************************************************
# Curves with prescribed j invariant

# This curve is defined if jj - 1728 is invertible, and is
# nonsingular if jj is also invertible.  It has
# c4 = - c6 = jj/(jj - 1728) and Delta = jj^2 / (jj - 1728)^3
# and j = jj (unless jj = 0, in which case j is undefined).

# jcurve[0] corresponds to the curve (y + x/2)^2 =(x + 1/12)^3

jcurve[jj_] :=
 {a1 -> 1, a2 -> 0, a3 -> 0,
  a4 -> 36/(1728 - jj), a6 -> 1/(1728 - jj)}

#**************************************************************
# Legendre curve.  This is defined where 2 is invertible and
# lm != 0 or 1.  The points of order two are
# P = [0:0:1] and Q = [1:0:1] and P+Q = [lm:0:1].  An invariant
# differential is given by om = -dx/y.  This is the universal
# example of a curve with an invariant differential and a
# level two structure such that x(P) = 0 and x(Q) = 1, where
# x is the unique choice adapted to om with a1 = a3 = 0.

# c4    = 16(1 - lm + lm^2)
# c6    = 32(lm - 2)(lm + 1)(2 lm - 1)
# Delta = 16 lm^2 (lm - 1)^2
# j     = 256 (1 - lm + lm^2)^3 / ((lm - 1)^2 lm^2)

legendre = {a1 -> 0, a3 -> 0, a2 -> -(lm + 1), a4 -> lm, a6 -> 0}

#**************************************************************
# Jacobi quartic.

# c4    = 2^6 3^4 (dl^2 + 3 ep)
# c6    = 2^9 3^6 (dl^3 - 9 dl ep)
# Delta = 2^12 3^12 ep (ep - dl^2)^2
# j     = 64 (dl^2 + 3 ep)^3 / (ep (ep - dl^2)^2)

jacobi_rule :=
 {a1 = 0, a3 = 0, a2 = 0, a4 = -108 dl^2 - 324 ep, a6 = 3888 dl ep - 432 dl^3}

xXY := (6*(-2*dl - dl^2*X^2 + 3*ep*X^2 + 2*dl*Y))/(-1 + dl*X^2 + Y);
yXY := (108*(dl^2 - ep)*X)/(-1 + dl*X^2 + Y);

Xxy = (6*(12*dl - x))/y
Yxy = (2592*dl^3 - 7776*dl*ep + 216*dl^2*x + 648*ep*x -
       36*dl*x^2 + y^2)/y^2

# In quartic form, the zero element is (0,1).

jacobiquartic = Y^2 - (1 - 2 dl X^2 + ep X^4)

# invariant differential

omegaXY = -dX/(6 Y)

# This can also be written as \sum_n P_n(dl,ep) X^(2n) dX
# where P_n(dl,ep)=P_n(dl/Sqrt[ep]) ep^(n/2) and P_n(t) is the
# Legendre polynomial.  These are defined by
# (1 -2xt +t^2)^(-1/2) = \sum_n P_n(x) t^n

# There is a canonical point P of order 2: (x,y)=(12 dl,0).
# Apparently this is the universal example of a curve with an
# invariant differential and a point of exact order 2 (where
# 2 is invertible).  Apparently also
# div(X) = [2]^{-1}[0] - [2]^{-1}[P]

# It seems that if six is invertible and we have a given point
# P of order 2, we can put the curve in the form
# y^2 = x^3 + a4 x + a6 and then P = (12 dl,0) say and we can
# define ep = -(4 a4 + 432 dl^2)/1296 to put the curve in
# Jacobi form.  I'm not sure this is right, though.

#**************************************************************

# Formal group formulae, valid up to terms of weight 7.

F = x0 + x1 + a1*x0*x1 - a2*x0^2*x1 + 2*a3*x0^3*x1 - 2*a1*a3*x0^4*x1 -
 2*a4*x0^4*x1 + 2*a1^2*a3*x0^5*x1 + 2*a2*a3*x0^5*x1 + 2*a1*a4*x0^5*x1 -
 2*a1^3*a3*x0^6*x1 - 4*a1*a2*a3*x0^6*x1 - 2*a3^2*x0^6*x1 -
 2*a1^2*a4*x0^6*x1 - 2*a2*a4*x0^6*x1 - 3*a6*x0^6*x1 - a2*x0*x1^2 -
 a1*a2*x0^2*x1^2 + 3*a3*x0^2*x1^2 + a2^2*x0^3*x1^2 - a1*a3*x0^3*x1^2 -
 4*a4*x0^3*x1^2 + a1^2*a3*x0^4*x1^2 + a1*a4*x0^4*x1^2 - a1^3*a3*x0^5*x1^2 -
 a1*a2*a3*x0^5*x1^2 - a1^2*a4*x0^5*x1^2 - 9*a6*x0^5*x1^2 + 2*a3*x0*x1^3 +
 a2^2*x0^2*x1^3 - a1*a3*x0^2*x1^3 - 4*a4*x0^2*x1^3 + a1*a2^2*x0^3*x1^3 +
 a1^2*a3*x0^3*x1^3 - 2*a2*a3*x0^3*x1^3 - a2^3*x0^4*x1^3 - a1^3*a3*x0^4*x1^3 -
 2*a1*a2*a3*x0^4*x1^3 + 4*a3^2*x0^4*x1^3 - a1^2*a4*x0^4*x1^3 +
 4*a2*a4*x0^4*x1^3 - 15*a6*x0^4*x1^3 - 2*a1*a3*x0*x1^4 - 2*a4*x0*x1^4 +
 a1^2*a3*x0^2*x1^4 + a1*a4*x0^2*x1^4 - a2^3*x0^3*x1^4 - a1^3*a3*x0^3*x1^4 -
 2*a1*a2*a3*x0^3*x1^4 + 4*a3^2*x0^3*x1^4 - a1^2*a4*x0^3*x1^4 +
 4*a2*a4*x0^3*x1^4 - 15*a6*x0^3*x1^4 + 2*a1^2*a3*x0*x1^5 + 2*a2*a3*x0*x1^5 +
 2*a1*a4*x0*x1^5 - a1^3*a3*x0^2*x1^5 - a1*a2*a3*x0^2*x1^5 -
 a1^2*a4*x0^2*x1^5 - 9*a6*x0^2*x1^5 - 2*a1^3*a3*x0*x1^6 -
 4*a1*a2*a3*x0*x1^6 - 2*a3^2*x0*x1^6 - 2*a1^2*a4*x0*x1^6 - 2*a2*a4*x0*x1^6 -
 3*a6*x0*x1^6

logF =
 x - (a1*x^2)/2 + (a1^2*x^3)/3 + (a2*x^3)/3 - (a1^3*x^4)/4 - (a1*a2*x^4)/2 -
 (a3*x^4)/2 + (a1^4*x^5)/5 + (3*a1^2*a2*x^5)/5 + (a2^2*x^5)/5 +
 (6*a1*a3*x^5)/5 + (2*a4*x^5)/5 - (a1^5*x^6)/6 - (2*a1^3*a2*x^6)/3 -
 (a1*a2^2*x^6)/2 - 2*a1^2*a3*x^6 - a2*a3*x^6 - a1*a4*x^6 + (a1^6*x^7)/7 +
 (5*a1^4*a2*x^7)/7 + (6*a1^2*a2^2*x^7)/7 + (a2^3*x^7)/7 +
 (20*a1^3*a3*x^7)/7 + (24*a1*a2*a3*x^7)/7 + (6*a3^2*x^7)/7 +
 (12*a1^2*a4*x^7)/7 + (6*a2*a4*x^7)/7 + (3*a6*x^7)/7

# ms[k] is the [k] series.

ms[-1] =
 -x + a1*x^2 - a1^2*x^3 + a1^3*x^4 + a3*x^4 - a1^4*x^5 - 3*a1*a3*x^5 +
 a1^5*x^6 + 6*a1^2*a3*x^6 + a2*a3*x^6 - a1^6*x^7 - 10*a1^3*a3*x^7 -
 4*a1*a2*a3*x^7 - 2*a3^2*x^7 + a1^7*x^8 + 15*a1^4*a3*x^8 +
 10*a1^2*a2*a3*x^8 + a2^2*a3*x^8 + 10*a1*a3^2*x^8 + a3*a4*x^8 - a1^8*x^9 -
 21*a1^5*a3*x^9 - 20*a1^3*a2*a3*x^9 - 5*a1*a2^2*a3*x^9 - 30*a1^2*a3^2*x^9 -
 5*a2*a3^2*x^9 - 5*a1*a3*a4*x^9

ms[2] =
 2*x + a1*x^2 - 2*a2*x^3 - a1*a2*x^4 + 7*a3*x^4 + 2*a2^2*x^5 - 6*a1*a3*x^5 -
 12*a4*x^5 + a1*a2^2*x^6 + 7*a1^2*a3*x^6 + 2*a2*a3*x^6 + 6*a1*a4*x^6 -
 2*a2^3*x^7 - 8*a1^3*a3*x^7 - 14*a1*a2*a3*x^7 + 4*a3^2*x^7 - 8*a1^2*a4*x^7 +
 4*a2*a4*x^7 - 54*a6*x^7

ms[3] =
 3*x + 3*a1*x^2 + a1^2*x^3 - 8*a2*x^3 - 12*a1*a2*x^4 + 39*a3*x^4 -
 6*a1^2*a2*x^5 + 24*a2^2*x^5 - 9*a1*a3*x^5 - 96*a4*x^5 - a1^3*a2*x^6 +
 44*a1*a2^2*x^6 + 30*a1^2*a3*x^6 - 57*a2*a3*x^6 - 48*a1*a4*x^6 +
 30*a1^2*a2^2*x^7 - 72*a2^3*x^7 - 36*a1^3*a3*x^7 - 168*a1*a2*a3*x^7 +
 234*a3^2*x^7 - 72*a1^2*a4*x^7 + 288*a2*a4*x^7 - 936*a6*x^7

ms[5] =
 5*x + 10*a1*x^2 + 10*a1^2*x^3 - 40*a2*x^3 + 5*a1^3*x^4 - 140*a1*a2*x^4 +
 310*a3*x^4 + a1^4*x^5 - 222*a1^2*a2*x^5 + 376*a2^2*x^5 + 306*a1*a3*x^5 -
 1248*a4*x^5 - 205*a1^3*a2*x^6 + 1740*a1*a2^2*x^6 + 620*a1^2*a3*x^6 -
 2130*a2*a3*x^6 - 3120*a1*a4*x^6 - 120*a1^4*a2*x^7 + 3750*a1^2*a2^2*x^7 -
 3560*a2^3*x^7 - 90*a1^3*a3*x^7 - 9540*a1*a2*a3*x^7 + 10540*a3^2*x^7 -
 5800*a1^2*a4*x^7 + 14240*a2*a4*x^7 - 33480*a6*x^7

Null[];


text: print "Hello world!";

max_stars_repo_name: jgke/miniplc

max_issues_repo_name: jgke/miniplc

max_forks_repo_name: jgke/miniplc


text:
check_set_generic("RR");
check_set_generic("CC");
check_set_generic("RZ");
check_set_generic("R",[[5]]);
check_set_generic("C",[[5]]);
check_set_generic("CP",[[5]]);
check_set_generic("sphere",[[5]]);
check_set_generic("RZ_set");
check_set_generic("RP_points");
check_set_generic("RP_lines");
check_set_generic("RP_conics");
check_set_generic("RP_arcs");
check_set_generic("simplices",[[5], [3]]);
check_set_generic("real_functions",[[{10, 20, 30, 40, 50}]]);
check_set_generic("vector_functions",[[4], [{10, 20, 30, 40, 50}]]);
check_set_generic("maps",[[{10, 20, 30, 40, 50}, {1, 2, 3}]]);
check_set_generic("binops",[[{0, 1}, {1, 2, 3}, {1, 2, 3}]]);
check_set_generic("epi",[[{10, 20, 30, 40, 50}, {1, 2, 3}]]);
check_set_generic("subsets",[[{10, 20, 30, 40, 50}]]);
check_set_generic("nonempty_subsets",[[{10, 20, 30, 40, 50}]]);
check_set_generic("proper_nonempty_subsets",[[{10, 20, 30, 40, 50}]]);
check_set_generic("big_subsets",[[{10, 20, 30, 40, 50}]]);
check_set_generic("pairs",[[{10, 20, 30, 40, 50}]]);
check_set_generic("pair_subsets",[[{10, 20, 30, 40, 50}]]);
check_set_generic("triples",[[{10, 20, 30, 40, 50}]]);
check_set_generic("partitions",[[{10, 20, 30, 40, 50}]]);
check_set_generic("permutations",[[5]]);
check_set_generic("shuffles",[[3, 4]]);
check_set_generic("nat_partitions",[[8, 3]]);
check_set_generic("pos_partitions",[[8, 3]]);
check_set_generic("rel",[[{11, 22, 33, 44}, {1, 2, 3}]]);
check_set_generic("autorel",[[{1, 2, 3}]]);
check_set_generic("equiv",[[{10, 20, 30, 40, 50}]]);
check_set_generic("preord",[[{10, 20, 30, 40, 50}]]);
check_set_generic("partord",[[{10, 20, 30, 40, 50}]]);
check_set_generic("total_preord",[[{10, 20, 30, 40, 50}]]);
check_set_generic("ord",[[{10, 20, 30, 40, 50}]]);
check_set_generic("cycord",[[{10, 20, 30, 40, 50}]]);
check_set_generic("linord",[[{10, 20, 30, 40, 50}]]);
check_set_generic("posint_intervals");
check_set_generic("interval_partitions",[[5]]);
check_set_generic("graphs",[[{10, 20, 30, 40, 50}]]);
check_set_generic("paths",[[{11, 12, 13, 14, 15, 16, 17, 18, 19, 20}], [{[11, 15], [11, 17], [11, 18], [11, 19], [11, 20], [12, 13], [12, 16], [12, 19], [12, 20], [13, 12], [13, 15], [13, 16], [13, 19], [13, 20], [14, 15], [14, 18], [14, 20], [15, 11], [15, 13], [15, 14], [15, 16], [15, 18], [15, 20], [16, 12], [16, 13], [16, 15], [16, 19], [17, 11], [17, 19], [17, 20], [18, 11], [18, 14], [18, 15], [18, 19], [19, 11], [19, 12], [19, 13], [19, 16], [19, 17], [19, 18], [20, 11], [20, 12], [20, 13], [20, 14], [20, 15], [20, 17]}], [4]]);
check_set_generic("trails",[[{11, 12, 13, 14, 15, 16, 17, 18, 19, 20}], [{[11, 15], [11, 17], [11, 18], [11, 19], [11, 20], [12, 13], [12, 16], [12, 19], [12, 20], [13, 12], [13, 15], [13, 16], [13, 19], [13, 20], [14, 15], [14, 18], [14, 20], [15, 11], [15, 13], [15, 14], [15, 16], [15, 18], [15, 20], [16, 12], [16, 13], [16, 15], [16, 19], [17, 11], [17, 19], [17, 20], [18, 11], [18, 14], [18, 15], [18, 19], [19, 11], [19, 12], [19, 13], [19, 16], [19, 17], [19, 18], [20, 11], [20, 12], [20, 13], [20, 14], [20, 15], [20, 17]}], [4]]);
check_set_generic("digraphs",[[{11, 22, 33, 44}]]);
check_set_generic("dipaths",[[{11, 12, 13, 14, 15, 16, 17, 18, 19, 20}], [{[11, 12], [11, 13], [11, 15], [11, 19], [12, 14], [12, 15], [12, 16], [12, 19], [13, 16], [13, 17], [13, 19], [14, 15], [14, 19], [15, 13], [15, 17], [15, 19], [16, 11], [16, 15], [16, 18], [16, 19], [17, 11], [17, 12], [17, 14], [17, 19], [18, 12], [19, 18], [19, 20], [20, 13], [20, 15], [20, 18]}], [4]]);
check_set_generic("ditrails",[[{11, 12, 13, 14, 15, 16, 17, 18, 19, 20}], [{[11, 12], [11, 13], [11, 15], [11, 19], [12, 14], [12, 15], [12, 16], [12, 19], [13, 16], [13, 17], [13, 19], [14, 15], [14, 19], [15, 13], [15, 17], [15, 19], [16, 11], [16, 15], [16, 18], [16, 19], [17, 11], [17, 12], [17, 14], [17, 19], [18, 12], [19, 18], [19, 20], [20, 13], [20, 15], [20, 18]}], [4]]);
check_set_generic("trees",[[{10, 20, 30, 40, 50}]]);
check_set_generic("full_trees",[[{10, 20, 30, 40, 50}]]);
check_set_generic("binary_trees",[[{10, 20, 30, 40, 50}]]);
check_set_generic("standard_stasheff_trees",[[6]]);
check_set_generic("binary_standard_stasheff_trees",[[6]]);
check_set_generic("stasheff_trees",[[{11, 22, 33, 44}]]);
check_set_generic("standard_stasheff_star",[[5]]);
check_set_generic("stasheff_star",[[{11, 22, 33, 44}]]);
check_set_generic("W",[[4], [{11, 22, 33, 44}]]);
check_set_generic("prime_W",[[4], [{11, 22, 33, 44}]]);
check_set_generic("SW",[[4], [{11, 22, 33, 44}]]);
check_set_generic("W1",[[{11, 22, 33, 44}]]);
check_set_generic("BW1",[[{11, 22, 33, 44}]]);
check_set_generic("SWW",[[3], [{10, 20, 30, 40, 50}, {10, 30, 50}]]);
check_set_generic("single_cubes",[[3]]);
check_set_generic("cubes",[[3], [{10, 20, 30, 40, 50}]]);
check_set_generic("one_cubes_prime",[[{10, 20, 30, 40, 50}]]);
check_set_generic("simplicial_interval",[[6]]);
check_set_generic("simplicial_maps",[[6, 5]]);
check_set_generic("simplicial_epi",[[5, 3]]);
check_set_generic("simplicial_mono",[[3, 5]]);
check_set_generic("simplex",[[{10, 20, 30, 40, 50}]]);
check_set_generic("simplex_interior",[[{10, 20, 30, 40, 50}]]);
check_set_generic("ord_simplex_interior",[[{10, 20, 30, 40, 50}]]);
check_set_generic("prime_simplex",[[{10, 20, 30, 40, 50}]]);
check_set_generic("prime_simplex_boundary",[[{10, 20, 30, 40, 50}]]);
check_set_generic("height_functions",[[{11, 12, 13, 14, 15, 16, 17, 18, 19, 20}]]);
check_set_generic("tree_height_functions",[[{11, 12, 13, 14, 15, 16, 17, 18, 19, 20}], [{{11}, {12}, {13}, {14}, {15}, {16}, {17}, {18}, {19}, {20}, {11, 13, 14}, {16, 17, 18, 19, 20}, {11, 12, 13, 14, 15, 16, 17, 18, 19, 20}}]]);
check_set_generic("F",[[3], [{10, 20, 30, 40, 50}]]);
check_set_generic("Fbar",[[3], [{10, 20, 30, 40, 50}]]);
check_set_generic("E",[[6]]);
check_set_generic("P",[[6]]);
check_set_generic("SEM",[[3], [{10, 11, 12, 13, 14, 15, 16}]]);
check_set_generic("ACP",[[3], [{10, 11, 12, 13, 14, 15, 16}]]);
check_set_generic("SCP",[[3], [{10, 11, 12, 13, 14, 15, 16}]]);
check_set_generic("ICP",[[3], [A]]);
check_set_generic("SCP2",[[3], [{10, 11, 12, 13, 14, 15, 16}], [B6]]);
check_set_generic("FFbar",[[3], [{10, 20, 30, 40, 50}]]);
check_set_generic("X/F3",[[4]]);
check_set_generic("Y/F3",[[4]]);
check_set_generic("D_plus/F3",[[4]]);
check_set_generic("D_minus/F3",[[4]]);
check_set_generic("D_zero/F3",[[4]]);
check_set_generic("comm",[[{10, 20, 30, 40, 50}]]);
check_set_generic("itloc",[[{1, 2, 3}], [3]]);
check_set_generic("Phi0/itloc",[[{10, 20, 30, 40, 50}]]);
check_set_generic("M11/itloc",[[{10, 20, 30, 40, 50}], [{10, 30, 50}]]);
check_set_generic("M1/itloc",[[{1, 2, 3, 4, 5}], [[{4}, {3, 4}, {1, 3}, {3, 4}]]]);


text:
go := module ()
  local undoStack, sr;
  export here, where, up, child, fac, trm, silent, undo, ModuleApply;

  where := [];

  up := proc()
    here, where := op(1,where)(here),
                   [op(2..-1,where)];
  end proc;

  child := proc(i :: integer)
    local _;
    here, where := op(i,here),
                   [unapply(subsop(i=_,here), _), op(where)];
  end proc;

  sr := proc(operator, t :: type)
    local s, r, _;
    if here :: operator then
      s, r := selectremove(type, here, t);
      here, where := s, [unapply(operator(_,r), _), op(where)];
    elif type(here, t) then
      where := [(_->_), op(where)];
    else
      here, where := operator(), [unapply(operator(_,here), _), op(where)];
    end if
  end proc;

  silent := proc()
    local instr, _, j;
    for instr in _passed do
      if instr :: list then
        thisproc(op(instr))
      elif instr :: integer then
        child(instr)
      elif instr = up then
        up()
      elif instr :: all(anything) then
        here, where := op(1,instr),
                       [unapply(subsindets(here,identical(op(1,instr)),(x->_)),
                                _),
                        op(where)];
      elif instr :: fac(type) then
        sr(`*`, op(1,instr))
      elif instr :: trm(type) then
        sr(`+`, op(1,instr))
      elif instr :: type then
        j := [seq(`if`(op(j,here) :: instr, j, NULL), j=1..nops(here))];
        if nops(j) = 1 then
          child(op(1,j))
        else
          error "not unique: %%", op(map(op, j, here))
        end if
      elif instr :: procedure then
        here := instr(here)
      end if
    end do
  end proc;

  undoStack := [];

  undo := proc()
    here, where := op(op(1,undoStack));
    undoStack := [op(2..-1,undoStack)]:
    here, where;
  end proc;

  ModuleApply := proc()
    undoStack := [[here, where], op(undoStack)];
    try silent(_passed);
    finally print(here, where);
    end
  end proc;
end module:


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

testinput := "6,10
0,14
9,10
0,3
10,4
4,11
6,0
6,12
4,1
0,13
10,12
3,4
3,0
8,4
1,10
2,14
8,10
9,0

fold along y=7
fold along x=5";

split:=StringTools:-Split(input, "\n"):
member("", split, 'loc'); loc;

coords := split[1..loc-1]:
folds := split[loc+1..-1];
coords := {map(p->[parse(p)], coords)[]}:
folds := parse~(map(s->StringTools:-SubstituteAll(s, "fold along ", ""), folds));

vis := coords->plots:-display(plottools:-point(map(c->[c[1],-c[2]], [coords[]]),
    symbol=solidbox, symbolsize=20), axes=boxed, scaling=constrained):

newcoords := coords:
for f in folds do
    if lhs(f) = 'y' then
        (c1,c2) := selectremove(p->p[2]<rhs(f), newcoords);
        c2 := map(p->[p[1], 2*rhs(f)-p[2]], c2);
        newcoords := c1 union c2;
    else
        (c1,c2) := selectremove(p->p[1]<rhs(f), newcoords);
        c2 := map(p->[2*rhs(f)-p[1], p[2]], c2);
        newcoords := c1 union c2;
    end if;
    if f = folds[1] then
        answer1 := nops(newcoords);
    end if;
end do:
answer1;
vis(newcoords);


text:
SL2_11 := table():

SL2_11["elements"] := {
 seq(seq(seq(seq(<<e*a|b>,<c|mods(e*(1+b*c)/a,11)>>,c=-5..5),b=-5..5),a=1..5),e in {-1,1}),
 seq( seq(seq(<<0|e*b>,<mods(-e/b,11)|d>>,d=-5..5),b=1..5),e in {-1,1})
}:

SL2_11["projective_line"] := {seq(i,i=-5..5),infinity};

SL2_11["action"] := proc(g,u)
 local a,b,c,d,nn,dd;

 a := g[1,1]; b := g[1,2];
 c := g[2,1]; d := g[2,2];
 if u = infinity then
  if c = 0 then
   return infinity;
  else
   return mods(a/c,11);
  fi;
 else
  nn := mods(a * u + b,11);
  dd := mods(c * u + d,11);
  if dd = 0 then
   return infinity;
  else
   return mods(nn/dd,11);
  fi;
 fi;
end:

SL2_11["first_hexad"] := {seq(mods(i^2,11),i=0..5)};

SL2_11["hexads"] :=
 {seq(map2(SL2_11["action"],g,SL2_11["first_hexad"]),g in SL2_11["elements"])};

# `is_element/steiner_systems`(SL2_11["projective_line"])(5,6)(SL2_11["hexads"]);


text:
program whilest;	{sample26}
var
 n,
 i,
 sum
 : integer;
begin
readln(n);i:=n;sum:=0;
while i>0 do begin
sum :=             sum
 + i;i:=i-1
end;writeln('Summention of 1 - ', n, '  is ', sum)
end.

max_issues_repo_name: shouth/LanguageProcessing


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

interface(quiet=true):

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

polyDegree := 4:

ERecSqrtmMin := 0.707:
ERecSqrtmMax := 1.4143:

xErrMin := evalf((1/(ERecSqrtmMax^2))):
xErrMax := evalf((1/(ERecSqrtmMin^2))):

printf("Using a %d degree polynomial for computing the iteration seed\n",polyDegree);

polyExact := numapprox[minimax](1/sqrt(x),x=1/2..2,[polyDegree,0],1,'err'):
poly := poly_exact(polyExact):
eps := numapprox[infnorm]((poly*sqrt(x))-1,x=xErrMin..xErrMax):

printf("Relative approximation error eps = 2^(%f)\n",log[2](abs(eps))):

printf("Writing tables...\n");

filename:="TEMPSQRT/sqrt.h":
fd:=fopen(filename, WRITE, TEXT):

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

for i from 0 to polyDegree do
	fprintf(fd, "\#define SQRTPOLYC%d %1.50e\n",i,coeff(poly,x,i)):
od:
fprintf(fd, "\#define TWO52 %1.50e\n",evalf(2^(52))):

fclose(fd):

filename:="TEMPSQRT/sqrt.sed":
fd:=fopen(filename, WRITE, TEXT):

fprintf(fd, "s/_epsilonApprox/%1.50e/g\n",eps):
fprintf(fd, "s/_ERecSqrtmMin/%1.50e/g\n",ERecSqrtmMin):
fprintf(fd, "s/_ERecSqrtmMax/%1.50e/g\n",ERecSqrtmMax):

for i from 0 to polyDegree do
   fprintf(fd, "s/_SQRTPOLYC%d/%1.50e/g\n",i,coeff(poly,x,i)):
od:

fclose(fd):

printf("... done\n");


text:
program assignboolean;
var x, y : boolean;
begin
    x := true;
    y := false;
end.


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


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

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

lspbe_f0 := s -> 1 + params_a_kappa*(1 - params_a_kappa/(params_a_kappa + params_a_mu*s^2))
            - (params_a_kappa+1)*(1-exp(-params_a_alpha*s^2)):
lspbe_f  := x -> lspbe_f0(X2S*x):

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


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

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

(* type: gga_exc *)

$define gga_x_rpw86_params
$include "gga_x_pw86.mpl"

lv_alpha := 0.02178:
lv_beta  := 1.15:
lv_muLV   := 0.8491/9:

lv_f0 := s ->
   + (1 + lv_muLV*s^2)/(1 + lv_alpha*s^6)
   + lv_alpha*s^6*pw86_f0(s)/(lv_beta + lv_alpha*s^6):

lv_f  := x -> lv_f0(X2S*x):

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