text:
# Base Parameter Regressor for Robot based on MDH frames
# Einleitung
# Erstellung einer parameterlinearen Regressorform in Newton Euler Bewegungsgleichung
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# chain -> Berechnung für eine serielle Struktur (nicht: Baumstruktur)
# fixb -> fixed base. Kein Floating base Modell. Dort ist diese Form der Minimalparameterform nicht möglich.
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# NewtonEuler -> Berechnung der Newton Euler Bewegungsgleichung
# regressor -> Regressorform (parameterlinear)
#
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
interface(rtablesize=100): # Zur Anzeige von größeren Vektoren
;
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
codegen_act := true:
codegen_opt := 2:
codeexport_grav := true:
codeexport_corvec := true:
codeexport_cormat := true:
codeexport_inertia := true:
codeexport_inertiaD := true:
codeexport_invdyn := true:
codeexport_act:= true:
codegen_dynpar := 2:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../robot_codegen_definitions/robot_env":
printf("Generiere Regressorform für %s\n", robot_name, codegen_dynpar):
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name, base_method_name):
read sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
kin_constraints_exist := kin_constraints_exist: # nur zum Abschätzen der Komplexität
;
# Term-Vereinfachungen einstellen
if not assigned(simplify_options) or simplify_options(10)=-1 then # Standard-Einstellungen:
  if not kin_constraints_exist then # normale serielle Ketten und Baumstrukturen
    use_simplify := 0: # Standardmäßig aus
  else # mit kinematischen Zwangsbedingungen
    use_simplify := 1: # standardmäßig simplify-Befehle anwenden
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(10): # zehnter Eintrag ist für Dynamik-Regressor

end if:
# Ergebnisse der Newton Euler- Bewegungsgleichung laden
read sprintf("../codeexport/%s/tmp/invdyn_%s_NewtonEuler_linkframe_par%d_maple.m", robot_name, base_method_name, codegen_dynpar):
f_i_i := f_i_i:
m_i_i := m_i_i:
tau_B := tau_B:
tau_J := tau_J:
# Mit diesem Arbeitsblatt werden die Vorwärtsrekursive für Fixed-Base Modelle generiert. Erkenne welche Basis-Modellierung aktiv ist
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):
end if:
# Die kinetischen und potentiellen Energien aus (2) und (3) stehen ab hier durch T_fixb und U_fixb zur Verfügung.
# Der Parametervektor 'PV2_vec' aus (13) wurde in 'robot_tree_floatb_twist_definitions.mw' aufgestellt.

# Parameterlinearisierung
# Parameterlinearisierung auf Basis von [HRL_IDR] (14) und (15)
# Linearisierung_Gelenkmomente
tauJ_regressor := Matrix(NQJ, 10*(NL-1)):

for i from 1 to NQJ do
  for j from 1 to 10*(NL-1) do
    tauJ_regressor[i,j] :=  diff(tau_J(i,1),PV2_vec[10+j,1]):
  end do:
end do:
# Terme vereinfachen
if use_simplify>=1 then
  tmp_t1:=time():
  tmp_l1 := length(tauJ_regressor):
  tauJ_regressor := simplify2(tauJ_regressor):
  tmp_t2:=time():
  tmp_l2 := length(tauJ_regressor):
  printf("%s: Gelenkmoment-Regressor vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t1):
end if:

# Linearisierung_Gelenkmomente(base)
#
tauB_regressor := Matrix(6,10*NL):
for i from 1 to 6 do
  for j from 1 to 10*NL do
    tauB_regressor[i,j] := diff(tau_B(i,1),PV2_vec[j,1]):
  end do:
end do:
# Terme vereinfachen
if use_simplify>=1 then
  tmp_t1:=time():
  tmp_l1 := length(tauB_regressor):
  tauB_regressor := simplify2(tauB_regressor):
  tmp_t2:=time():
  tmp_l2 := length(tauB_regressor):
  printf("%s: Basis-Kraft/Moment-Regressor vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t1):
end if:

# Linearisierung_Schnittmomente
m_regressor := Matrix(3*NL,10*NL):

m_i_i_vec := Matrix(3*(NL),1):
for i from 1 to NL do
  for j from 1 to 3 do
    m_i_i_vec[3*(i-1)+j] := m_i_i[j,i]:
  end do:
end do:

for i from 1 to 3*NL do
  for j from 1 to 10*NL do
    m_regressor[i,j] := diff~(m_i_i_vec(i,1),PV2_vec[j,1]):
  end do:
end do:

# Terme vereinfachen
if use_simplify>=1 then
  tmp_t1:=time():
  tmp_l1 := length(m_regressor):
  m_regressor := simplify2(m_regressor):
  tmp_t2:=time():
  tmp_l2 := length(m_regressor):
  printf("%s: Schnittmoment-Regressor vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t1):
end if:

# Linearisierung_Schnittkräfte
f_regressor := Matrix(3*NL,10*NL):
f_i_i_vec := Matrix(3*(NL),1):
for i from 1 to NL do
  for j from 1 to 3 do
    f_i_i_vec[3*(i-1)+j] := f_i_i[j,i]:
  end do:
end do:

for i from 1 to 3*NL do
  for j from 1 to 10*NL do
    f_regressor[i,j] := diff~(f_i_i_vec(i,1),PV2_vec[j,1]):
  end do:
end do:
# Terme vereinfachen
if use_simplify>=1 then
  tmp_t1:=time():
  tmp_l1 := length(f_regressor):
  f_regressor := simplify2(f_regressor):
  tmp_t2:=time():
  tmp_l2 := length(f_regressor):
  printf("%s: Schnittkraft-Regressor vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t1):
end if:

# Export
# Maple Export
save tauJ_regressor, sprintf("../codeexport/%s/tmp/fixb_NewtonEuler_tauJ_regressor_maple.m", robot_name):
save tauB_regressor, sprintf("../codeexport/%s/tmp/fixb_NewtonEuler_tauB_regressor_maple.m", robot_name):
save m_regressor, sprintf("../codeexport/%s/tmp/fixb_NewtonEuler_m_regressor_maple.m", robot_name):
save f_regressor, sprintf("../codeexport/%s/tmp/fixb_NewtonEuler_f_regressor_maple.m", robot_name):
printf("Maple-Ausdrücke exportiert. %s\n", FormatTime("%Y-%m-%d %H:%M:%S")):
# Matlab Export
if codegen_act then
  MatlabExport(convert_t_s(tauJ_regressor), sprintf("../codeexport/%s/tmp/invdyn_%s_NewtonEuler_tauJ_regressor_matlab.m", robot_name, expstring), codegen_opt):
  MatlabExport(convert_t_s(tauB_regressor), sprintf("../codeexport/%s/tmp/invdyn_%s_NewtonEuler_tauB_regressor_matlab.m", robot_name, expstring), codegen_opt):
  MatlabExport(convert_t_s(m_regressor), sprintf("../codeexport/%s/tmp/invdyn_%s_NewtonEuler_m_regressor_matlab.m", robot_name, expstring), codegen_opt):
  MatlabExport(convert_t_s(f_regressor), sprintf("../codeexport/%s/tmp/invdyn_%s_NewtonEuler_f_regressor_matlab.m", robot_name, expstring), codegen_opt):
end if
;

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


text:
# Maple Code Generator

##
## Copyright (C) 2015-2017 Maplesoft
## Authors:    Behzad Samadi <bsamadi@maplesoft.com>
## Created:    September 2017
## Version:    0.2
## Keywords:   Maple, Code Generation, Model Predictive Control (MPC)
##
## Procedures:
##             CGMRES
##

MapleCodeGenerator := module()
description "Maple code generation for numerical optimization";
option package;

export DotProduct,
       GMRES,
       GradientDescent,
       LMGMRES,
       Minimize1D,
       NewName,
       NewtonsMethod,
       NewtonsMethodGMRES,
       nfGMRES,
       Norm2;

$include "DotProduct.mm"
$include "GMRES.mm"
$include "GradientDescent.mm"
$include "LMGMRES.mm"
$include "Minimize1D.mm"
$include "NewName.mm"
$include "NewtonsMethod.mm"
$include "NewtonsMethodGMRES.mm"
$include "nfGMRES.mm"
$include "Norm2.mm"

end module; # MCG

ext: mpl

max_stars_repo_name: SymbolicNMPC/MapleCodeGenerator

max_issues_repo_name: SymbolicNMPC/MapleCodeGenerator

max_forks_repo_name: SymbolicNMPC/MapleCodeGenerator


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

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

ff     := z -> params_a_BLOC_a + params_a_BLOC_b*z:
mkappa := (x, t) -> params_a_kappa:

$include "tpss_x.mpl"

(* Equation (5) *)

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

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
ctfm_expr :=  exp(x) - exp(lv) + ((1/(1-x/rv) - 1/(1-lv/rv)))^(rv-x+1) ;
ctfmp_expr := ((diff(ctfm_expr,x)));
octfmp_expr := -(diff(1/ctfmp_expr,x));

read ("/home/arman/FD/FD.mpl"):
MFD(); grid_functions:={};
Gen_Eval_Code(ctfm_expr,input="c",proc_name="init_ctfm2");
Gen_Eval_Code(ctfmp_expr,input="c",proc_name="init_ctfmp2");
Gen_Eval_Code(octfmp_expr,input="c",proc_name="init_octfmp2");


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

a2 := 146/2025:
a3 := -73/405:
a4 := 0.131957187845257783631757384393: (* DD + 100/(81*81*kappa) *)

kappa := 0.804:

xx := (p, qt) -> MU_GE*p + a2*qt*qt + a3*qt*p + a4*p*p:

f := (rs, x, t, u) -> 1 + kappa - \
  kappa^2/(kappa + xx(X2S^2*x^2, 6*X2S^2*t - 9/20 - X2S^2*x^2/12)):


text:
text.usetex : True
text.latex.preamble: \usepackage{bm}
axes.grid : False
figure.figsize : 6, 5.25
figure.dpi  : 150
font.size : 16
xtick.major.width : 1.5
xtick.minor.width : 1.5
ytick.major.width : 1.5
ytick.minor.width : 1.5

figure.subplot.left    : 0.15  # 0.12  # the left side of the subplots of the figure
figure.subplot.right   : 0.85  # 0.88  # the right side of the subplots of the figure
figure.subplot.bottom  : 0.15  # 0.14  # the bottom of the subplots of the figure
figure.subplot.top     : 0.95  # 0.86  # the top of the subplots of the figure

max_stars_repo_name: warrickball/bluebell

max_issues_repo_name: warrickball/bluebell

max_forks_repo_name: warrickball/bluebell


text:
repeat process ACV_cubes {
/execute @e[type=armor_stand,name=ACV_Cube] ~ ~ ~ detect ~ ~1 ~ minecraft:lava * kill @e[type=armor_stand,name=ACV_Cube,r=0,c=1]

/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID=0] ~ ~ ~ detect ~ ~1 ~ minecraft:air * tp @e[type=armor_stand,name=ACV_Cube,r=0,c=1] ~ ~-0.25 ~
/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID=0] ~ ~ ~ detect ~ ~1 ~ minecraft:stone_pressure_plate * tp @e[type=armor_stand,name=ACV_Cube,r=0,c=1] ~ ~-0.25 ~
/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID=0] ~ ~ ~ detect ~ ~1 ~ minecraft:carpet * tp @e[type=armor_stand,name=ACV_Cube,r=0,c=1] ~ ~-0.25 ~
/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID=0] ~ ~ ~ detect ~ ~1 ~ minecraft:tallgrass * tp @e[type=armor_stand,name=ACV_Cube,r=0,c=1] ~ ~-0.25 ~

/scoreboard players add @a ACV_CUBE_ID 0
/scoreboard players add @e[type=armor_stand,name=ACV_Cube] ACV_CUBE_ID 0
/scoreboard players tag @a remove ACV_HoldingArrow
/scoreboard players tag @a add ACV_HoldingArrow {SelectedItem:{id:"minecraft:arrow",tag:{display:{Name:"Picking up Cubes"}}}}
/scoreboard players set @a[score_ACV_CUBE_ID_min=1,tag=!ACV_HoldingArrow] ACV_CUBE_ID 0
/execute @p[score_ACV_CUBE_ID=0,tag=ACV_HoldingArrow] ~ ~ ~ scoreboard players operation ACV_Last_ID ACV_CUBE_ID %= 5 ACV_CONSTANT
/execute @p[score_ACV_CUBE_ID=0,tag=ACV_HoldingArrow] ~ ~ ~ scoreboard players add ACV_Last_ID ACV_CUBE_ID 1
/scoreboard players operation @r[type=player,score_ACV_CUBE_ID=0,tag=ACV_HoldingArrow] ACV_CUBE_ID += ACV_Last_ID ACV_CUBE_ID

/scoreboard players tag @a remove ACV_CubeLinked
/scoreboard players tag @e[type=armor_stand,name=ACV_Cube] remove ACV_CubeLinked
/execute @p[score_ACV_CUBE_ID_min=1,score_ACV_CUBE_ID=1] ~ ~ ~ execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,score_ACV_CUBE_ID=1] ~ ~ ~ scoreboard players tag @e[score_ACV_CUBE_ID_min=1,score_ACV_CUBE_ID=1] add ACV_CubeLinked
/execute @p[score_ACV_CUBE_ID_min=2,score_ACV_CUBE_ID=2] ~ ~ ~ execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=2,score_ACV_CUBE_ID=2] ~ ~ ~ scoreboard players tag @e[score_ACV_CUBE_ID_min=2,score_ACV_CUBE_ID=2] add ACV_CubeLinked
/execute @p[score_ACV_CUBE_ID_min=3,score_ACV_CUBE_ID=3] ~ ~ ~ execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=3,score_ACV_CUBE_ID=3] ~ ~ ~ scoreboard players tag @e[score_ACV_CUBE_ID_min=3,score_ACV_CUBE_ID=3] add ACV_CubeLinked
/execute @p[score_ACV_CUBE_ID_min=4,score_ACV_CUBE_ID=4] ~ ~ ~ execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=4,score_ACV_CUBE_ID=4] ~ ~ ~ scoreboard players tag @e[score_ACV_CUBE_ID_min=4,score_ACV_CUBE_ID=4] add ACV_CubeLinked
/execute @p[score_ACV_CUBE_ID_min=5,score_ACV_CUBE_ID=5] ~ ~ ~ execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=5,score_ACV_CUBE_ID=5] ~ ~ ~ scoreboard players tag @e[score_ACV_CUBE_ID_min=5,score_ACV_CUBE_ID=5] add ACV_CubeLinked

/scoreboard players set @e[type=armor_stand,name=ACV_Cube,tag=!ACV_CubeLinked] ACV_CUBE_ID 0
/execute @a[tag=!ACV_CubeLinked] ~ ~ ~ /scoreboard players operation @e[type=armor_stand,name=ACV_Cube,r=3,c=1] ACV_CUBE_ID += @p[tag=!ACV_CubeLinked] ACV_CUBE_ID

/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,score_ACV_CUBE_ID=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,c=1] @p[score_ACV_CUBE_ID_min=1,score_ACV_CUBE_ID=1]
/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=2,score_ACV_CUBE_ID=2] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,c=1] @p[score_ACV_CUBE_ID_min=2,score_ACV_CUBE_ID=2]
/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=3,score_ACV_CUBE_ID=3] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,c=1] @p[score_ACV_CUBE_ID_min=3,score_ACV_CUBE_ID=3]
/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=4,score_ACV_CUBE_ID=4] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,c=1] @p[score_ACV_CUBE_ID_min=4,score_ACV_CUBE_ID=4]
/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=5,score_ACV_CUBE_ID=5] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,c=1] @p[score_ACV_CUBE_ID_min=5,score_ACV_CUBE_ID=5]

// r=1 wegen https://bugs.mojang.com/browse/MC-90967
/execute @a[rym=169,ry=-168,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~0 ~-0.25 ~-2
/execute @a[rym=-167,ry=-146,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~0.7653668647301795434569199680608 ~-0.25 ~-1.8477590650225735122563663787936
/execute @a[rym=-145,ry=-123,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~1.4142135623730950488016887242097 ~-0.25 ~-1.4142135623730950488016887242097
/execute @a[rym=-122,ry=-101,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~1.8477590650225735122563663787936 ~-0.25 ~-0.7653668647301795434569199680608
/execute @a[rym=-100,ry=-78,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~2 ~-0.25 ~0
/execute @a[rym=-77,ry=-56,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~1.8477590650225735122563663787936 ~-0.25 ~0.7653668647301795434569199680608
/execute @a[rym=-55,ry=-33,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~1.4142135623730950488016887242097 ~-0.25 ~1.4142135623730950488016887242097
/execute @a[rym=-32,ry=-11,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~0.7653668647301795434569199680608 ~-0.25 ~1.8477590650225735122563663787936
/execute @a[rym=-10,ry=11,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~0 ~-0.25 ~2
/execute @a[rym=12,ry=33,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~-0.7653668647301795434569199680608 ~-0.25 ~1.8477590650225735122563663787936
/execute @a[rym=34,ry=56,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~-1.4142135623730950488016887242097 ~-0.25 ~1.4142135623730950488016887242097
/execute @a[rym=57,ry=78,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~-1.8477590650225735122563663787936 ~-0.25 ~0.7653668647301795434569199680608
/execute @a[rym=79,ry=101,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~-2 ~-0.25 ~0
/execute @a[rym=102,ry=123,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~-1.8477590650225735122563663787936 ~-0.25 ~-0.7653668647301795434569199680608
/execute @a[rym=124,ry=146,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~-1.4142135623730950488016887242097 ~-0.25 ~-1.4142135623730950488016887242097
/execute @a[rym=147,ry=168,score_ACV_CUBE_ID_min=1] ~ ~ ~ tp @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID_min=1,r=1,c=1] ~-0.7653668647301795434569199680608 ~-0.25 ~-1.8477590650225735122563663787936

/kill @e[type=armor_stand,name=cube_update]
/execute @e[type=armor_stand,name=ACV_Cube,score_ACV_CUBE_ID=0] ~ ~ ~ summon armor_stand ~ ~1.25 ~ {CustomName:"cube_update",Invisible:1b,Invulnerable:1b,Marker:1b}
}


text:
with(ListTools):

####################
f1 := 2*x^2-x*y+4:
f2 := x*y-2*y^2+4:
sub := [x=x1+I*x2, y=y1+I*y2]:
f1 := expand(subs( sub, f1)):
f2 := expand(subs( sub, f2)):
p1 := Re(f1) assuming real:
p2 := Im(f1) assuming real:
p3 := Re(f2) assuming real:
p4 := Im(f2) assuming real:

################
nld := proc(d, nb_vars)
local l, i, j, p, x, vars;
l := [];
vars := [];
for i from 1 to nb_vars do
    p := -1;
    for j from 1 to nb_vars do
        if i=j then
           p := p+x||j^d;
        else
           p := p+x||j;
        end if;
    end do;
    l := [op(l), p];
    vars := [op(vars), x||i];
end do;
return l, vars;
end proc:

################
nqleven := proc(nb_vars, even_d)
local x, vars, eqs, i;
vars := [ x||1 ];
eqs := [ x||1^even_d-2 ];
for i from 2 to nb_vars do
    vars := [ x||i, op(vars)];
    eqs := [op(eqs), x||i^even_d + x||i^(even_d/2)- x||(i-1)];
od;
return eqs, vars;
end proc:

################
simplenql := proc(nbvar, d)
local x, eqs, i, vars;
vars := [ x||1 ];
eqs := [ x||1^d-2 ];
for i from 2 to nbvar do
    vars := [ x||i, op(vars) ];
    eqs := [op(eqs), x||i^d - x||(i-1)];
end do;
return eqs, vars;
end proc:

####################
cm :=
[
    [   #(1) 4-body-homog ((1)-(47) are from changbo's folder BCLM-cm14)
        [-2*p^3+2*p^3*phi^3-4*phi^3*s*p^2+5*phi^3*s^3*p-phi^3*s^5,
         -2*s*p^3-2*phi^3*s^2+phi^3*s^4-3*phi^3*s^2*p+2*phi^3*p,
         -2*s^2+s^4-4*s^2*p+phi^2+1+4*p
        ],
        [phi,s,p]
    ],
    [
        #(2) 5-body-homog
        [-6*p^3+4*p^3*phi^3+15*phi^3*s^3*p-3*phi^3*s^5-12*phi^3*s*p^2-3*phi^3*s*p+phi^3*s^3,
        -9*phi^3*s^2*p-5*phi^3*s^2-6*s*p^3+3*phi^3*s^4+5*phi^3*p,
        -12*s^2*p-6*s^2+3*s^4+4*phi^2+3+12*p
        ],
        [phi,s,p]
    ],
    [
        #(3) Arnborg-Lazard
        [x^2*y*z + x*y^2*z + x*y*z^2 + x*y*z + x*y + x*z + y*z,
        x^2*y^2*z + x*y^2*z^2 + x^2*y*z + x*y*z + y*z + x + z,
        x^2*y^2*z^2 + x^2*y^2*z + x*y^2*z + x*y*z + x*z + z + 1
        ],
        [z,y,x]
    ],
    [
        #(4) Arnborg-Lazard-rev
        [x^2*y*z + x*y^2*z + x*y*z^2 + x*y*z + x*y + x*z + y*z,
        x^2*y^2*z + x*y^2*z^2 + x^2*y*z + x*y*z + y*z + x + z,
        x^2*y^2*z ^2 + x^2*y^2*z + x*y^2*z + x*y*z + x*z + z + 1
        ],
        [x,y,z]
    ],
    [
        #(5) Barry
        [-x^5+y^5-3*y-1,
         5*y^4-3,
         -20*x+y-z
        ],
        [z,y,x]
    ],
    [
        #(6) Caprasse
        [y^2*z+2*x*y*t-2*x-z,
         -x^3*z+4*x*y^2*z+4*x^2*y*t+2*y^3*t+4*x^2-10*y^2+4*x*z-10*y*t+2,
         2*y*z*t+x*t^2-x-2*z,
         -x*z^3+4*y*z^2*t+4*x*z*t^2+2*y*t^3+4*x*z+4*z^2-10*y*t-10*t^2+2
        ],
        [t,z,y,x]
    ],
    [
        #(7) Caprasse-Li
        [y^2*z+2*x*y*t-2*x-z,
         -x^3*z+4*x*y^2*z+4*x^2*y*t+2*y^3*t+4*x^2-10*y^2+4*x*z-10*y*t+2,
         2*y*z*t+x*t^2-x-2*z,
         -x*z^3+4*y*z^2*t+4*x*z*t^2+2*y*t^3+4*x*z+4*z^2-10*y*t-10*t^2+2
        ],
        [x,y,z,t]
    ],
    [
        #(8) chemical-reaction
        [2 - 7*x1 + x1^2*x2 - 1/2*(x3 - x1),
        6*x1 - x1^2*x2 - 5*(x4 - x2 ),
        2 - 7*x3 + x3^2*x4 - 1/2*(x1 - x3 ),
        6*x3 - x3^2*x4 + 1 + 1/2*(x2 - x4 )
        ],
        [x1,x2,x3,x4]
    ],
    [
        #(9) circles
        [product( (x-'i'), 'i'=1..10) +1/100,
         product ( (x-'i')^2+(y-'i')^2-2, 'i'=1..5)+1/1000
        ],
        [x,y]
    ],
    [
        #(10) cyclic-5
        [a + b + c + d + e,
         a*b + b*c + c*d + d*e + e*a,
         a*b*c + b*c*d + c*d*e + d*e*a + e*a*b,
         a*b*c*d + b*c*d*e + c*d*e*a + d*e*a*b  + e*a*b*c,
         a*b*c*d*e - 1
        ],
        [e,d,c,b,a]
    ],
    [
        #(11) Czapor-Geddes-Wang
        [b-2,
         c*c*q*q-2*d*c*q*q+d*d*q*q-2*c*c*q+2*d*p*c*q+2*b*p*c*q+2*d*c*q-2*d*d*p*q+2*b*p*q-2*p*q-2*b*d*d*q-2*b*q+2*q+c*c-2*d*p*c-2*b*d*c+b*b*d*d*p*p-2*b*d*d*p*p+d*d*p*p+2*b*b*p*p+2*p*p-2*b*b*d*d*p+4*b*d*d*p-4*b*b*p+4*b*p+b*b*d*d+2*b*b-4*b+2,
         -b*p*c*q-p*c*q+b*c*q-2*c*q+2*d*p*q+d*q+p*c+b*b*c-b*c+2*c+b*b*d*p*p+b*d*p*p-2*d*p*p-2*b*b*d*p+b*d*p-d*p+b*b*d-2*b*d,
         -2*p*q-2*q-b*b*p*p+2*p*p+2*b*b*p-b*b+2,
         3*c*c*q*q-3*d*d*q*q-4*q*q-6*c*c*q+6*b*d*p*c*q+6*b*d*c*q+6*d*d*p*q+3*c*c+6*b*d*p*c-6*b*d*c+3*b*b*d*d*p*p-3*d*d*p*p+b*b*p*p-6*b*b*d*d*p-2*b*b*p+4*p+3*b*b*d*d+b*b
        ],
        [b,d,p,c,q]
    ],
    [
        #(12) fabfaux
        [30752*x^3 + 46128*y*x^2 - 216256*x + 46128*x*y^2 + 141980 + 30752*y^3- 216256*y,
        30752*x^3 + 46128*z*x^2 - 216256*x + 46128*x*z^2 + 141980 + 30752*z^3+ 216256*z,
        46128*y^3 + 46128*x*y^2 + 46128*z*y^2 + 46128*y*x^2 - 432512*y +46128*x*z*y + 46128*z^2*y + 46128*x^3 - 432512*x + 425940 +46128*z*x^2 - 432512*z + 46128*x*z^2 + 46128*z^3
        ],
        [z,y,x]
    ],
    [
        #(13) geometric-constraints
        [1/100 - 4*s*(s - 1)*(s - b)*(s - c),
         1/5 - b*c,
         2*s - 1 - b - c
        ],
        [s,b,c]
    ],
    [
        #(14) GonzalezGonzalez
        [x^3 + y*x^2 + x + z,
         x^2*y +3*x +2,
         x^2-y*x+z
        ],
        [z,y,x]
    ],
    [
        #(15) Katsura-4
        [2*x^2 + 2*y^2 + 2*z^2 + 2*t^2 + u^2 - u ,
        2*x*y + 2*y*z + 2*z*t + 2*t*u - t ,
        2*x*z + 2*y*t + t^2 + 2*z*u - z ,
        2*x*t + 2*z*t + 2*y*u - y ,
        2*x + 2*y + 2*z + 2*t + u - 1
        ],
        [u,t,z,y,x]
    ],
    [
        #(16) lhlp1
        [2*x^7+2-3*x^2-x^4,
         x^2+2*y^2-5,
         x*z-1
        ],
        [x,y,z]
    ],
    [
        #(17) lhlp2
        [y*(320+1600*y^4-240*y^5-471*y^6+36*y^7-48*y^2+36*y^8),
        -40*y^2+3*y^3+6*y^4+8*x,
        -8*y*z-8-40*y^4+3*y^5+6*y^6
        ],
        [x,y,z]
    ],
    [
        #(18) lhlp3
        [2450*x^6-1241*x^4+196*x^2-49,
         86*x^2+35*y*x^2-77*y-14,
         175*x^4+70*x^2*z-76*x^2-154*z+21
        ],
        [x,y,z]
    ],
    [
        #(19) lhlp4
        [x^4+y^4-1,
         x^5*y^2-4*x^3*y^3+x^2*y^5-1
        ],
        [x,y]
    ],
    [
        #(20) lhlp5
        [-7*x*y*z+6*y*z-14*x*z+9*z-3*x*y-12*y-x+1,
         2*x*y*z-y*z+14*z+15*x*y+14*y-15*x,
         -8*x*y*z+11*y*z-12*x*z-5*z+15*x*y+2*y+10*x-14
        ],
        [x,y,z]
    ],
    [
        #(21) lhlp6
        [p1,p2,p3,p4
        ],
        [y1,x1,y2,x2]
    ],
    [
        #(22) neural-network
        [1 - c*x - x*y^2 - x*z^2,
         1 - c*y - y*x^2 - y*z^2,
         1 - c*z - z*x^2 - z*y^2,
         8*c^6 + 378*c^3 - 27
        ],
        [z,c,y,x]
    ],
    [
        #(23) nld-3-4
        [x^3 + y + z + t- 1,
         x + y^3 + z + t -1,
         x + y + z^3 + t-1,
         x + y + z + t^3 -1
        ],
        [t,z,y,x]
    ],
    [
        #(24) nld-3-5
        [x^3 + y + z + t + u- 1,
         x + y^3 + z + t + u-1,
         x + y + z^3 + t + u-1,
         x + y + z + t^3 + u-1,
         x + y + z + t + u^3 -1
        ],
        [u,t,z,y,x]
    ],
    [
        #(25) nld-4-5
        [x^4 + y + z + t + u- 1,
         x + y^4 + z + t + u-1,
         x + y + z^4 + t + u-1,
         x + y + z + t^4 + u-1,
         x + y + z + t + u^4 -1
        ],
        [x,y,z,t,u]
    ],
    [
        #(26) nld-7-3
        [x^7 + y + z - 1,
         x + y^7 + z - 1,
         x + y + z^7 - 1
        ],
        [z,y,x]
    ],
    [
        #(27) nld-8-3
        [x^8 + y + z - 1,
         x + y^8 + z - 1,
         x + y + z^8 - 1
        ],
        [z,y,x]
    ],
    [
        #(28) nld-9-3
        [x^9 + y + z - 1,
         x + y^9 + z - 1,
         x + y + z^9 - 1
        ],
        [z,y,x]
    ],
    [
        #(29) nld-10-3
        [x^10 + y + z - 1,
         x + y^10 + z - 1,
         x + y + z^10 - 1
        ],
        [z,y,x]
    ],
    [
        #(30) nql-5-4
        nqleven(5,4)
    ],
    [
        #(31) nql-10-2
        nqleven(10,2)
    ],
    [
        #(32) nql-10-4
        nqleven(10,4)
    ],
    [
        #(33) nql-15-2
        nqleven(15,2)
    ],
    [
        #(34) p3p-special
        [x^2 + y^2 - x*y - 1,
         y^2 + z^2 - y*z - a^2,
         z^2 + x^2 - z*x - b^2,
         a^2 - 1 + b - b^2,
         3*b^6 + 56*b^4 - 122*b^3 + 56*b^2 + 3
        ],
        [b,a,x,y,z]
    ],
    [
        #(35) PlateForme2d-easy
        [(x1+5)^2+ (y1-0)^2 - 1,
         (x2-5)^2+ (y2-0)^2 - 1,
         (x3)^2  + (y3-5)^2 - 1,
         (x1-x2)^2 + (y1-y2)^2 - 3,
         (x1-x3)^2 + (y1-y3)^2 - 3,
         (x2-x3)^2 + (y2-y3)^2 - 3
        ],
        [y3,x3,y2,x2,y1,x1]
    ],
    [
        #(36) r-5
        [a^2 + a,
         a*b + b + a*b^2 +a,
         b^2*c + c + b*c^3 + b,
         c^3*d + d + c*d^4 + c,
         d^4*e + e + d*e^5 + d
        ],
        [e,d,c,b,a]
    ],
    [
        #(37) r-6
        [a^2 + a,
         a*b + b + a*b^2 +a,
         b^2*c + c + b*c^3 + b,
         c^3*d + d + c*d^4 + c,
         d^4*e + e + d*e^5 + d,
         e^5*f + f + e*f^6 + e
        ],
        [f,e,d,c,b,a]
    ],
    [
        #(38) Reif
        [x4*x13+x5*x14+x6*(1-x13-x14),
         x4*x15+x5*x16-x6*(x15+x16) ,
         x7*x13+x8*x14+x9*(1-x13-x14) ,
         x7*x15+x8*x16-x9*(x15+x16)-1 ,
         x10*x13+x11*x14+x12*(1-x13-x14) ,
         x10*x15+x11*x16-x12*(x15+x16) ,
         x1*x13+x2*x14+x3*(1-x13-x14) ,
         x1*x15+x2*x16-x3*(x15+x16) ,
         x1*x4*x13+x2*x5*x14+x3*x6*(1-x13-x14)-1 ,
         x1*x4*x15+x2*x5*x16-x3*x6*(x15+x16) ,
         x1*x7*x13+x2*x8*x14+x3*x9*(1-x13-x14) ,
         x1*x7*x15+x2*x8*x16-x3*x9*(x15+x16) ,
         x1*x10*x13+x2*x11*x14+x3*x12*(1-x13-x14) ,
         x1*x10*x15+x2*x11*x16-x3*x12*(x15+x16)-1
        ],
        [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16]
    ],
    [
        #(39) Rose
        [7*y^4 - 20*x^2,
         2160*x*x*z^4+1512*x*z^4+315*z^4-4000*x*x-2800*x-490,
         -10080000*x^4*z^3-28224000*x^3*z^3-15288000*x*x*z^3-1978032*x*z^3-180075*z^3-23520000*x^4*y*z*z-41395200*x^3*y*z*z-26726560*x*x*y*z*z-7727104*x*y*z*z-852355*y*z*z+40320000*x^6*y*y*z+28800000*x^5*y*y*z+21168000*x^3*y*y*z+4939200*x*x*y*y*z+347508*x*y*y*z+67200000*x^5*y^3+94080000*x^4*y^3+40924800*x^3*y^3+2634240*x*x*y^3-2300844*x*y^3-432180*y^3
        ],
        [x,y,z]
    ],
    [
        #(40) simple-nql-20-30
        simplenql(20, 30)
    ],
    [
        #(41) Takeuchi-Lu
        [2*x1*(2 - x1 - y1 ) + x2 - x1,
         2*x2*(2 - x2 - y2 ) + x1 - x2,
         2*y1*(5 - x1 - 2*y1 ) + y2 - y1,
         y2*(3 - 2*x2 - 4*y2 ) + y1 - y2
        ],
        [x1,x2,y1,y2]
    ],
    [
        #(42) Themos-net
        [10000*d - 323536,
         (x-48)^2+(y-89)^2+(28)^2-(48*a+89*b+28*c)^2-d^2,
         (x-77)^2+(y-3)^2+(37)^2-(77*a+3*b+37*c)^2-d^2,
         (x-49)^2+(y-23)^2+(57)^2-(49*a+23*b+57*c)^2-d^2,
         a*x+b*y,
         a^2+b^2+c^2-1
        ],
        [d,c,b,a,y,x]
    ],
    [
        #(43) Trinks-2
        [45*y+35*u-165*v-36 ,
         35*y+25*z+40*t-27*u ,
         25*y*u-165*v^2+15*x-18*z+30*t ,
         15*y*z+20*t*u-9*x ,
         -11*v^3+x*y+2*z*t ,
         -11*u*v+3*v^2+99*x ,
         10000*v^3+6600*v+2673
        ],
        [v,u,z,t,y,x]
    ],
    [
        #(44) Trinks-difficult
        [45*y+35*u-165*v-36 ,
         35*y+25*z+40*t-27*u ,
         25*y*u-165*v^2+15*x-18*z+30*t ,
         15*y*z+20*t*u-9*x ,
         -11*v^3+x*y+2*z*t ,
         -11*u*v+3*v^2+99*x
        ],
        [v,u,z,t,y,x]
    ],
    [
        #(45) Uteshev-Bikker
        [x^2+x*y+y^2-2*x*z-4*y*z+3*z^2-3*x*t+2*y*t+t^2-3*x-2*y+3*z-2*t-2 ,
        2*x^2-x*y+y^2-x*z-y*z-6*z^2-x*t+y*t-5*z*t-3*t^2-5*x+y+5*z+2*t+5 ,
        -3-3*x*y+2*x*z+x*t^2-5*x*z^2-5*z^2*t-3*x*t-2*z*t+x*y*z+x*y*t-x^2*z+x^2-y^2+2*z^2+11*z-2*t-x+y+x^3+y^3-3*z^3+2*t^3-3*t^2-5*y^2*z+7*y*z^2 ,
        -15+2*x*y+11*x*t^2+5*x*z^2-z*t-4*x*y*z+6*x*y*t-x^2*z+3*x^2+2*y^2-z^2+4*z-10*t-35*x-14*y-x^3+6*y^3+15*z^3+4*t^3+5*t^2+6*y^2*z+4*y*z^2-x*z*t+6*x^2*y-12*x*y^2-7*y^2*t+2*y*t
        ],
        [t,z,y,x]
    ],
    [
        #(46) wilkinson20
        [mul( (x-i),i=0..20 )+1/100],
        [x]
    ],
    [
        #(47) wilkinsonxy
        [mul( (x-i),i=1..5 )+1/100,
	     mul( (y-x-i),i=1..5)
        ],
        [x,y]
    ],
    [
        #(48) trivial-5 ((48) is from changbo's folder RCtest-zerodim)
        [x^2 -1,
         y^2 - 1,
         (x-1)*(y-1)
        ],
        [y,x]
    ],
    [
        #(49) nld-4-4
        nld(4, 4)
    ],
    [
        #(50) nld-5-4
        nld(5, 4)
    ],
    [
        #(51) nld-6-4
        nld(6, 4)
    ],
    [
        #(52) nld-5-3
        nld(5, 3)
    ]
]:

max_stars_repo_name: lihaokun/StrongSfTriDec

max_issues_repo_name: lihaokun/StrongSfTriDec

max_forks_repo_name: lihaokun/StrongSfTriDec


text:
# The function brent_fsolve uses Brent's algorithm to find a root of
# a real valued function of one variable.   This code is based on the
# Matlab implementation by John Burkardt, which is in turn based on the
# original Fortran implementation by Richard Brent.  However, this
# version is designed for the case where we want to retain a lot
# of additional information that is generated in the process of
# evaluating the relevant function.
#
# The first argument F should be a function which accepts a single
# real argument t and returns a table F(t).  This should contain an
# entry F(t)["err"], which we want to be zero.  The second and third
# arguments should be numbers a and b such that F(a)["err"] and
# F(b)["err"] have opposite signs.  Other entries in F(t) can contain
# auxiliary information generated during the calculation.
#
# If F(a) and/or F(b) have already been computed, then they can be
# supplied as the arguments Fa_ and Fb_.  Otherwise, these arguments
# should be set to false.

#@ brent_fsolve
brent_fsolve := proc(F,a,b,Fa_,Fb_,tol_,epsilon_)
 local tol,tol1,epsilon,c,sa,sb,sc,Fa,Fb,Fc,fa,fb,fc,d,e,m,p,q,r,s,steps,ret;

 tol := `if`(nargs > 5,tol_,10.^(-20));
 epsilon := `if`(nargs > 6,epsilon_,10.^(2 - Digits));

 sa := a;
 sb := b;
 if nargs < 4 or Fa_ = false then Fa := eval(F(sa)); else Fa := eval(Fa_); fi;
 if nargs < 5 or Fb_ = false then Fb := eval(F(sb)); else Fb := eval(Fb_); fi;
 fa := Fa["err"];
 fb := Fb["err"];

 if not((fa <= 0 and 0 <= fb) or
        (fb <= 0 and 0 <= fa)) then
  error("function has the same sign at endpoints");
 fi;

 c := sa;
 fc := fa;
 Fc := eval(Fa);
 e := sb - sa;
 d := e;

 steps := 0;

 while true do
  steps := steps+1;

  if ( abs ( fc ) < abs ( fb ) ) then
   sa := sb;
   sb := c;
   c := sa;
   fa := fb;
   fb := fc;
   fc := fa;
   Fa := eval(Fb);
   Fb := eval(Fc);
   Fc := eval(Fa);
  fi;

  tol1 := 2.0 * epsilon * abs ( sb ) + tol;
  m := 0.5 * ( c - sb );

  if ( abs ( m ) <= tol or fb = 0.0 ) then
   break;
  fi;

  if ( abs ( e ) < tol or abs ( fa ) <= abs ( fb ) ) then
   e := m;
   d := e;
  else
   s := fb / fa;
   if ( sa = c ) then
    p := 2.0 * m * s;
    q := 1.0 - s;
   else
    q := fa / fc;
    r := fb / fc;
    p := s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
    q := ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
   fi;

   if ( 0.0 < p ) then
    q := - q;
   else
    p := - p;
   fi;

   s := e;
   e := d;

   if ( 2.0 * p < 3.0 * m * q - abs ( tol * q ) and p < abs ( 0.5 * s * q ) ) then
    d := p / q;
   else
    e := m;
    d := e;
   fi;
  fi;

  sa := sb;
  fa := fb;
  Fa := eval(Fb);

  if ( tol < abs ( d ) ) then
   sb := sb + d;
  elif ( 0.0 < m ) then
   sb := sb + tol;
  else
   sb := sb - tol;
  fi;

  Fb := eval(F(sb));
  fb := Fb["err"];

  if ( ( 0.0 < fb and 0.0 < fc ) or ( fb <= 0.0 and fc <= 0.0 ) ) then
   c := sa;
   fc := fa;
   Fc := eval(Fa);
   e := sb - sa;
   d := e;
  fi;
 od;

 ret := eval(Fb);
 ret["steps"] := steps;

 return [sb,eval(ret)];
end:


text:
$ifndef _UTILS_
$define _UTILS_

tpop:=proc(t,k)
	local r;
	r:=t[k][1];
	t[k]:=t[k] minus {r};
	return r;
end proc:

# 表达式按照复杂度升序排序
sortByComplexity:=proc(_s::list,{index::boolean:=false})
    local s,t,i,n;
    if index then
        s:=_s;
        t:=table();
        n:=numelems(s);
        for i from 1 to n do
            tappend(t,s[i],i);
        end do;
        s:=ListTools[Reverse](SolveTools[SortByComplexity](s));
        return map[2](tpop,t,s);
    else
        return ListTools[Reverse](SolveTools[SortByComplexity](_s));
    end if;
end proc:

# 按照集合拓展table键值
tappend:=proc(t,k,v)
    if assigned(t[k]) then
        t[k]:=t[k] union {v};
    else
        t[k]:={v};
    end if;
end proc:

# 对象按键值分类
# 推荐对键做convert/global处理，以消除局部变量相等的问题。
collectObj:=proc(s,key,{output:=[val]::{[ind],[val],[ind,val]}})
    local t,v,res;
    t:=table();
    for v in s do
        tappend(t,key(v),v);
    end do;
    res:=();
    if (ind in output) then
        res:=res,[indices(t,nolist)];
    end if;
    if (val in output) then
        res:=res,[entries(t,nolist)];
    end if;
    return res;
end proc:

# 对象按键值唯一化
# 推荐对键做convert/global处理，以消除局部变量相等的问题
uniqueObj:=proc(s,key,{index::boolean:=false})
    local t,i,n;
    t:=table();
    n:=numelems(s);
    for i from 1 to n do
        t[key(s[i])]:=`if`(index,i,s[i]);
    end do;
    return [entries(t,nolist)];
end proc:

# 选择目标函数值最小的元素
MinSelect:=proc(v::list,fun)
    return v[min[index]([seq(fun(x),x in v)])];
end proc:

# 寻找特定值的下标
find:=proc(v::list,k)
    return map(i->`if`(v[i]=k,i,NULL),[seq(1..numelems(v))]);
end proc:

$endif


text:
restart;
writeto(terminal);
NULL;
S := readdata("s.txt", float, 3)[1];
X1 := S[1];
X2 := S[2];
T := S[3];
Y_str := parse(readdata("y.txt", string, 3)[1][1], statement);
y := unapply(Y_str, x1, x2, t);
y(0, 0, 42);
L_str := parse(readdata("L_eq.txt", string, 100)[1][1], statement);
G_str := parse(readdata("G.txt", string, 100)[1][1], statement);
c := 1;
L := unapply(L_str, x1, x2, t);
u := (x1, x2, t) -> L(x1, x2, t);
G := unapply(G_str, x1, x2, t);
;
u(0.5, 0.5, 0.5);
G(0, 1/5, 0.3);
y(0, 0, 42);

plotsetup(maplet);
with(plots);
animate3d(y(x1, x2, t), x1 = -X1 .. X1, x2 = -X2 .. X2, t = 0 .. 1, color = sin(x1*x2));
Y0rl := Matrix([[1/2, 1/2, 0], [1/3, 1/3, 0], [1/4, 1/4, 0], [1/5, 1/5, 0], [1/6, 1/6, 0]]);
L0 := 5;
Y0i := Vector(5);
for i to 5 do
    Y0i[i] := y(Y0rl[i, 1], Y0rl[i, 2], Y0rl[i, 3]);
end do;
YGrl := Matrix([[1/2, 1, 1/7], [0, 0, 1/2], [0, 1/3, 1/5], [1/2, 0, 1/4], [0, 1/2, 1/8]]);
LG := 5;
YGi := Vector(5);
for i to 5 do
    YGi[i] := evalf(y(YGrl[i, 1], YGrl[i, 2], YGrl[i, 3]));
end do;
Sm0 := Matrix([[1/2, 1/2, -1], [1/3, 1/3, -2]]);
M0 := 2;
SmG := Matrix([[1.1, 1.5, 1/2], [-0.5, 2, 1/2]]);
MG := 2;

yI := 0;
for i1 to 10 do
    for i2 to 10 do for i3 to 10 do yI1 := piecewise(c(i3/10)^4 + (-(i1/10)^2 - (i2/10)^2) <> 0, 0.001*1/(2*Pi*c*(c(t - i3/10)^4 + (-(x1 - i1/10)^2 - (x2 - i2/10)^2)))*L(i1/10, i2/10, i3/10), 0); yI := yI + yI1; end do; end do;
end do;
yInf := unapply(yI, x1, x2, t);
evalf(yInf(0, 0.65, 0.3));
Y0 := Vector(L0);
for l to L0 do
    Y0[l] := Y0i[l] - evalf(yInf(Y0rl[l, 1], Y0rl[l, 2], Y0rl[l, 3]));
end do;
YG := Vector(LG);
for l to LG do
    YG[l] := YGi[l] - evalf(yInf(YGrl[l, 1], YGrl[l, 2], YGrl[l, 3]));
end do;
Ybar := Vector(L0 + LG);
for i to L0 do
    Ybar[i] := Y0[i];
end do;
for i to LG do
    Ybar[L0 + i] := YG[i];
end do;
A11 := Matrix(L0, M0);
A := Matrix(L0 + LG, M0 + MG);
for i to L0 do
    for j to M0 do A11[i, j] := evalf(G(Y0rl[i, 1] - Sm0[j, 1], Y0rl[i, 2] - Sm0[j, 2], Y0rl[i, 3] - Sm0[j, 3])); A[i, j] := A11[i, j]; end do;
end do;
A11;
A12 := Matrix(L0, MG);
for i to L0 do
    for j to MG do A12[i, j] := evalf(G(Y0rl[i, 1] - SmG[j, 1], Y0rl[i, 2] - SmG[j, 2], Y0rl[i, 3] - SmG[j, 3])); A[i, j + M0] := A12[i, j]; end do;
end do;
i := 2;
j := 1;
A12[i, j] := evalf(G(Y0rl[i, 1] - SmG[j, 1], Y0rl[i, 2] - SmG[j, 2], Y0rl[i, 3] - SmG[j, 3]));
A12;
A21 := Matrix(LG, M0);
for i to LG do
    for j to M0 do A21[i, j] := evalf(G(YGrl[i, 1] - Sm0[j, 1], YGrl[i, 2] - Sm0[j, 2], YGrl[i, 3] - Sm0[j, 3])); A[L0 + i, j] := A21[i, j]; end do;
end do;
A21;
A22 := Matrix(LG, MG);
for i to LG do
    for j to MG do A22[i, j] := evalf(G(YGrl[i, 1] - SmG[j, 1], YGrl[i, 2] - SmG[j, 2], YGrl[i, 3] - SmG[j, 3])); A[i + L0, j + M0] := A22[i, j]; end do;
end do;
A22;
A;
# Псевдообернена:
with(LinearAlgebra);
At := Multiply(Transpose(A), MatrixInverse(Multiply(A, Transpose(A))));
ufinal := Multiply(At, Ybar);
u0 := Vector(M0);
y0 := 0;
for i to M0 do
    u0[i] := ufinal[i];
end do;
for i to M0 do
    y0 := y0 + G(x1 - Sm0[i, 1], x2 - Sm0[i, 2], t - Sm0[i, 3])*u0[i];
end do;
uG := Vector(MG);
for i to MG do
    uG[i] := ufinal[M0 + i];
end do;
yG := 0;
for i to MG do
    yG := yG + G(x1 - SmG[i, 1], x2 - SmG[i, 2], t - SmG[i, 3])*uG[i];
end do;
yfinal := unapply(evalf(yInf(x1, x2, t) + y0 + yG), x1, x2, t);

save yfinal, "y_final.log";
animate3d(yfinal(x1, x2, t), x1 = -X1 .. X1, x2 = -X2 .. X2, t = 0 .. T, color = sin(x1*x2));


text:
#ifndef SERVER_STUB_H__
#define SERVER_STUB_H__

#include <ccf/session_request.h>
#include <ccf/server.h>
#include <string>
#include "server/proto.h"

#define SVR_IMPL(NAME, req, zone) \
	void svr_##NAME(ccf::session_request<NAME/*, ccf::address*/> req, \
			ccf::auto_zone& zone)

namespace server {

%self.each do |msg|
%next unless msg.id
SVR_IMPL([%msg.name%], req, zone);
%end

static inline void dispatch(ccf::method_t method, ccf::msgobj param,
		ccf::session_responder response,
		/*const ccf::address& from, */ccf::auto_zone& z)
{
	switch(method) {
	%self.each do |msg|
	%next unless msg.id
	case [%msg.name%]::method:
		svr_[%msg.name%](ccf::session_request<[%msg.name%]/*, ccf::address*/>(
				param, response/*, from*/), z);
		break;
	%end
	default:
		throw std::runtime_error("unknown method");  // FIXME
	}
}

}  // namespace server

#endif /* server/stub.h */
%# vim: syntax=mplex

max_stars_repo_name: frsyuki/ccf

max_issues_repo_name: frsyuki/ccf

max_forks_repo_name: frsyuki/ccf


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

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

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

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

params_a_gamma := (1 - log(2))/Pi^2:
params_a_BB    := 1:
$include "gga_c_pbe.mpl"

ff := (z, t) -> mphi(z)^(params_a_alpha*t^3):

f  := (rs, z, xt, xs0, xs1) ->
  f_pw(rs, z) + ff(z, tp(rs, z, xt))*fH(rs, z, tp(rs, z, xt)):


text:
with(Groebner):

################ proc PolylistToNormal
# input:
# Polylist: a nonempty list of nonzero polynomials which are different from each other (list)
# VarList: a list of variables which appear in Polylist (list)
# output: If Polylist has finite solutions, then the output is a strong triangular decomposition (STD) of Polylist.
#         (a row vector of some row vectors)
#         Else FAIL.
PolylistToNormal := proc(PolyList, VarList)
local i, k, ans, TriList, j;
for i from 1 to nops(PolyList) do
    if type(PolyList[i], constant)=true then
       return Vector[row]();
    end if;
end do;
TriList := IsTri(PolyList, VarList);
if TriList[1]=true then
    if IniConstant(PolyList, TriList[2])=0 then
        if nops(PolyList)=nops(VarList) then
            ans := Vector[row](nops(VarList));
            for i from 1 to nops(TriList[3]) do
                k := TriList[3][i];
                ans[-i] := PolyList[k];
            end do;
            for i from 2 to nops(VarList) do
                ans[i] := NormalForm(ans[i], [seq(ans[j], j=1..i-1)], plex(op(VarList)));
            end do;
            return Vector[row](1, fill = ans);
        else
            error "FAIL1: This is not a zero-dimensional system.";
        end if;
    else
        return SubPolylistToNormal(PolyList, VarList, 1, 1);
    end if;
else
    return SubPolylistToNormal(PolyList, VarList, 2, 1);
end if;
end proc:

################ proc SubPolylistToNormal
# GBmethod: 0: fgb+convert
#           1: buchberger
#           2: the first time 0, else 1 (def)
# SATGBmethod: 0: fgb+convert
#              1: buchberger
SubPolylistToNormal := proc(PolyList, VarList, GBmethod := 2, SATGBmethod := 1)
local ComputedSet, GBmethod1, SplitSet, P, GB_Grad, GB_Lex, wCharList, AbnormalNum, IniNum, GB, AddF, i, GB_mainvar, lindex, lnum;
ComputedSet := Vector[row]();   # ComputedSet is a row vector of some row vectors
SplitSet := Vector[row]();      # SplitSet is a row vector of some lists
SplitSet(1) := PolyList;
GBmethod1 := GBmethod;
lindex := 1;
lnum := 1;
while lindex<=lnum do
    P := SplitSet[lindex];
    lindex := lindex+1;
    if GBmethod1=1 then
        GB_Lex := Basis(P, plex(op(VarList)), method = Buchberger);
    else
        GB_Grad := Basis(P, tdeg(op(VarList)), method = fgb);
        GB_Lex  := Basis(GB_Grad, plex(op(VarList)), method = convert);
        if GBmethod1=2 then
            GBmethod1 := 1;
        end if;
    end if;
    if GB_Lex <> [1] then
       GB_mainvar := [seq(MainVar(GB_Lex[i], VarList), i=1..nops(GB_Lex))];
       if lindex=2 then
           if not(IsZeroDim(GB_Lex, VarList)) then
              error "FAIL2: This is not a zero-dimensional system.";
           end if;
       end if;
       wCharList := wChar(GB_Lex, VarList, GB_mainvar); # wCharList is a row vector
       AbnormalNum := IsNormal(wCharList, VarList);
       if AbnormalNum = 0 then
            ComputedSet(numelems(ComputedSet)+1) := wCharList;
       else
            IniNum := SquareFree(primpart(lcoeff(wCharList[AbnormalNum], VarList[-AbnormalNum])));
            GB := SaturIdeal(IniNum, [seq(wCharList[i], i=1..AbnormalNum-1)], VarList, SATGBmethod); # GB is a row vector
            AddF := Vector[row]();
            for i from 1 to numelems(GB) do
                if NormalForm(GB[i], [seq(wCharList[i], i=1..AbnormalNum-1)], plex(op(VarList)))<>0 then
                    AddF(numelems(AddF)+1) := GB[i];
                end if;
            end do;
            SplitSet(++lnum):= [op(GB_Lex), seq(AddF[i], i=1..numelems(AddF))];
            SplitSet(++lnum):= [op(GB_Lex), IniNum];
       end if;
    end if;
end do;
return ComputedSet;
end proc:

################ proc NormalToSquarefree
# input:
# NormalChain: a strong chain (a row vector)
# VarList: a list of variables which appear in Polylist (list)
# output:
# ans: a strong square-free triangular decomposition (SSFTD) of NormalChain
#         (a row vector of some row vectors)
NormalToSquarefree := proc(NormalChain, VarList, SATGBmethod:=1)
local NormalSet, ans, P, AbNum, temp, tempSet, f, GB, AddF, i, r1, r2, lindex, lnum;
NormalSet := Vector[row]();   # NormalSet is a vector of some row vectors
NormalSet(1) := NormalChain;
ans := Vector[row]();         # ans is a vector of some row vectors
lindex := 1;
lnum := 1;
while lindex<=lnum do
   P := NormalSet[lindex];           # P is a row vector
   lindex := lindex+1;
   AbNum := IsSquareFree(P, VarList);
   if AbNum=0 then
        ans(numelems(ans)+1) := P;
   elif AbNum=1 then
        temp := factors(P[1]);
        tempSet := [seq(temp[2][i][1], i = 1 .. nops(temp[2]))];
        for i from 1 to nops(tempSet) do
            P(1) := tempSet[i];
            NormalSet(++lnum) := Vector[row](P);
        end do;
    else
        f := SquareFree(diff(P[AbNum], VarList[-AbNum]));
        GB := SaturIdeal(f, [seq(P[i], i=1..AbNum)], VarList, SATGBmethod); # GB is a row vector
        AddF := Vector[row]();
        for i from 1 to numelems(GB) do
            if NormalForm(GB[i], [seq(P[i], i=1..AbNum)], plex(op(VarList)))<>0 then
                AddF(numelems(AddF)+1) := GB[i];
            end if;
        end do;
        r1 := SubPolylistToNormal([seq(P[i], i=1..numelems(P)), f], VarList); # r1 is a row vector of some row vectors
        r2 := SubPolylistToNormal([seq(P[i], i=1..numelems(P)), seq(AddF[i], i=1..numelems(AddF))], VarList);
        # r2 is a row vector of some row vectors
        for i from 1 to numelems(r1) do
            NormalSet(++lnum) := r1[i];
        end do;
        for i from 1 to numelems(r2) do
            NormalSet(++lnum) := r2[i];
        end do;
    end if;
end do;
return ans;
end proc:

################ proc PolysetToSfNormal
# input:
# Polylist: a nonempty list of nonzero polynomials which are different from each other (list)
# VarList: a nonempty list of variables
# output:
# output: If Polylist has finite solutions, then the output is a strong square-free triangular decomposition (SSFTD) of Polylist.
#         (a row vector of some row vectors)
#         Else FAIL.
PolylistToSfNormal := proc(PolySet, VarList, SATGBmethod:=1)
local NormalDec, ans, i, r, j;
NormalDec := PolylistToNormal(PolySet, VarList);   # NormalDec is a row vector of some row vectors
if NormalDec=FAIL then
   return FAIL;
else
   ans := Vector[row]();     # ans is a row vector of some row vectors
   for i from 1 to numelems(NormalDec) do
       r := NormalToSquarefree(NormalDec[i], VarList, SATGBmethod);    # r is a row vector of some row vectors
       for j from 1 to numelems(r) do
           ans(numelems(ans)+1) := r[j];
       end do;
   end do;
   return ans;
end if;
end proc:

################ proc IsTri
IsTri := proc(PolyList, VarList)
local JudgeList, i, k, mainvarList;
JudgeList := [seq(-1, i=1..nops(VarList))];
mainvarList := [seq(MainVar(PolyList[i], VarList), i=1..nops(PolyList))];
for i from 1 to nops(PolyList) do
    member(mainvarList[i], VarList, 'k');
    if JudgeList[k]=-1 then
        JudgeList[k] := i;
    else
        return [false];
    end if;
end do;
return [true, mainvarList, JudgeList];
end proc:

################ proc wChar
# input:
# GroeBasis: the reduced lex Groebner basis w.r.t. VarList[1]>...>VarList[n]
# VarList: the variable list of GroeBasis
# mainvar: the main variable list of GroeBasis such that the main variable of GroeBasis[i] is mainvar[i]
# output:
# wCharList: the W-characteristic list of GroeBasis (row vector)
wChar := proc(GroeBasis, VarList, mainvar)
local k, i, PolySet, wCharList, wCharPrep, j;
PolySet := Vector[row](nops(VarList));
for k from 1 to nops(VarList) do
    PolySet[k] := Vector[row]();
end do;
for i from 1 to nops(GroeBasis) do
    member(mainvar[i], VarList, 'k');
    PolySet[k](numelems(PolySet[k])+1) := GroeBasis[i];
end do;
wCharList := Vector[row](nops(VarList));
for i from 1 to nops(VarList) do
    wCharPrep := PolySet[i][1];
    for j from 2 to numelems(PolySet[i]) do
        if TestOrder(LeadingMonomial(PolySet[i][j], plex(op(VarList))), LeadingMonomial(wCharPrep, plex(op(VarList))), plex(op(VarList))) then
           wCharPrep := PolySet[i][j];
        end if;
    end do;
    wCharList[-i] := wCharPrep;
end do;
return wCharList;
end proc:

################ proc IsNormal
# input:
# TriSet: a zero-dim triangle set T=[T1,...,Tn] such that 0<lv(T1)<...<lv(Tn) (row vector)
# VarList: the list of main variables of TriSet such that the main variable of TriSet[i] is VarList[-i] (list)
# output:
# NormNum: if Triset is normal, then return 0
# else return the smallest integer (>1) such that [T1,...,Tk] is abnormal
IsNormal := proc(TriSet, VarList)
local NormNum, i;
NormNum := 0;
for i from 1 to numelems(TriSet) do    # Note that Triset[1] is a univariate polynomial.
    if type(lcoeff(TriSet[i], VarList[-i]), constant)=false then
       NormNum := i;
       break;
    end if;
end do;
return NormNum;
end proc:

################ proc IsNormalSpe
# input:
# PolyList: a polynomial list which is a out-of-order triangle set (a list)
# mainvarList: lv(PolyList[i]) = mainvarList[i]
# output:
# NormNum: if all initials of PolyList are constants, then return 0; else return 1
IniConstant := proc(PolyList, mainvarList)
local i;
for i from 1 to nops(PolyList) do    # Note that Triset[1] is a univariate polynomial.
    if type(lcoeff(PolyList[i], mainvarList[i]), constant)=false then
       return 1;
    end if;
end do;
return 0;
end proc:

################ proc IsSquareFree
# input:
# TriSet: a triangle set T=[T1,...,Tn] such that 0<lv(T1)<...<lv(Tn) (a row vector)
# (a list)
# VarList: the list of main variables of TriSet such that the main variable of TriSet[i] is VarList[-i]
# output:
# NormNum: if Triset is square-free, then return 0
# else return the smallest integer such that [T1,...,Tk] is not squarefree
IsSquareFree := proc(TriSet, VarList)
local NormNum, j, i, dis, Tri;
NormNum := 0;
if primpart(discrim(TriSet[1], VarList[-1]))=0 then
   return 1;
else
   for i from 2 to numelems(TriSet) do
        dis[i] := primpart(discrim(TriSet[i], VarList[-i]));
        Tri := [seq(TriSet[j], j=1..i-1)];
        if ResTri(dis[i], Tri, VarList)=0 then
           NormNum := i; break;
        end if;
        if NormNum <> 0 then break; end if;
    end do;
    return NormNum;
end if;
end proc:

################ proc ResTri
# input:
# P : a polynomial
# TriSet: a triangle set T=[T1,...,Tn] such that 0<lv(T1)<...<lv(Tn)   (a list)
# VarList: the list of main variables of TriSet such that the main variable of TriSet[i] is VarList[-i]
# output:
# F: Res(P, TriSet)
ResTri := proc(P, TriSet, VarList)
local F, n,i,j,temp,tempSet1,tempSet2;
F := P;
for n from nops(TriSet) by -1 to 1 do
    F := primpart(resultant(F, TriSet[n], VarList[-n]));
end do;
return F;
end proc:

################ proc SaturIdeal
# input:
# Poly : a polynomial
# PolyList: a list of polynomials (a list)
# SATGBmthod: 0: fgb+convert; 1: buchberger
# output:
# ans: the lex reduced Groebner basis of <op(PolyList)>:Poly^{infinity} (a vector)
SaturIdeal := proc(Poly, PolyList, VarList, SATGBmethod:=0)
local Z, GB_Grad, GB_Lex, i, ans;
ans := Vector[row]();
if SATGBmethod=0 then
    GB_Grad := Basis([op(PolyList), 1-Poly*Z], tdeg(Z, op(VarList)), method = fgb);
    GB_Lex  := Basis(GB_Grad, plex(Z, op(VarList)), method = convert);
else
    GB_Lex  := Basis([op(PolyList), 1-Poly*Z], plex(Z, op(VarList)), method = buchberger);
end if;
for i from 1 to nops(GB_Lex) do
    if has(GB_Lex[i], Z)=false then
       ans(numelems(ans)+1) := GB_Lex[i];
    end if;
end do;
return ans;
end proc:

########## proc Class, proc MainVar, proc WuInitial
# input:
# f: a polynomial
# ord: a list of variables of f such that ord[1]>ord[2]>...
########## proc Class
Class := proc(f::polynom, ord::list)
local i;
option cache;
for i from 1 to nops(ord) do
	if has(f,ord[i]) then
		return i;
	end if;
end do;
return 0;
end proc:

########## proc MainVar
MainVar := proc(f::polynom, ord::list)
local i;
option cache;
for i from 1 to nops(ord) do
	if has(f, ord[i]) then
		return ord[i];
	end if;
end do;
return 0;
end proc:

########## proc WuInitial
WuInitial := proc(f::polynom, ord::list)
option cache;
if Class(f, ord)=0 then
	return 1;
else
	return primpart(lcoeff(f, MainVar(f,ord)));
end if;
end proc:

########## proc SquareFree
SquareFree := proc(poly)
option cache;
local FactorList,i;
FactorList := sqrfree(poly);
return mul(FactorList[2][i][1], i = 1 .. nops(FactorList[2]));
end proc:

########## proc IsZeroDim
IsZeroDim := proc(GBlex, VarList)
local JudgeList,i,IndSet,j;
JudgeList := [seq(0, i=1..nops(VarList))];
for i from 1 to nops(GBlex) do
   IndSet := indets(LeadingMonomial(GBlex[i], plex(op(VarList))));
   if nops(IndSet)=1 then
       for j from 1 to nops(VarList) do
           if IndSet={VarList[j]} then
               JudgeList[j] := 1;
               break;
           end if;
       end do;
   end if;
end do;
for i from 1 to nops(JudgeList) do
   if JudgeList[i]=0 then
       return false;
   end if;
end do;
return true;
end proc:


text:
program brainfuck;
var
    size: integer;
    cells: array[30000] of char;
    ptr: integer;
    src: array[30000] of char;
    cursol: integer;

procedure load;
var
    i, jump: integer;
begin
    read(size);
    i := 0;
    jump := 0;
    while i < size do
    begin
        read(src[i]);
        if src[i] = '[' then
            jump := jump + 1
        else if (src[i] = ']') and (jump > 0) then
            jump := jump - 1
        else
        begin
            writeln('mismatching square parenthesis');
            return
        end;
        i := i + 1
    end;
    if jump > 0 then
    begin
        writeln('unclosed square parenthesis');
        return
    end;
end;

begin
    call load;
    cursol := 0;
    ptr := 0;
    while cursol < size do
    begin
        if src[cursol] = '>' then
            ptr := ptr + 1
        else if src[cursol] = '<' then
            ptr := ptr - 1
        else if src[cursol] = '+' then
            cells[ptr] := char(integer(cells[ptr]) + 1)
        else if src[cursol] = '-' then
            cells[ptr] := char(integer(cells[ptr]) - 1)
        else if src[cursol] = '.' then
            write(cells[ptr])
        else if src[cursol] = ',' then
            read(cells[ptr])
        else if (src[cursol] = '[') and (integer(cells[ptr]) = 0) then
            while src[cursol] <> ']' do
                cursol := cursol + 1
        else if (src[cursol] = ']') and (integer(cells[ptr]) <> 0) then
            while src[cursol] <> '[' do
                cursol := cursol - 1;
        cursol := cursol + 1;
    end
end.

max_stars_repo_name: shouth/LanguageProcessing

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: shouth/LanguageProcessing


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

sys := [
diff(N(t), t) = -delta_NE*N(t)*P(t) - mu_N*N(t),
diff(S(t), t) = -mu_LE*S(t)*E(t) + delta_EL*S(t) - S(t)^2*mu_LL - S(t)*delta_LM,
diff(M(t), t) = S(t)*delta_LM - M(t)*mu_M,
diff(P(t), t) = rho_P*P(t)^2 - S(t)*P(t)*mu_PL - E(t)*mu_PE*P(t) - P(t)*mu_P,
diff(E(t), t) = -mu_EE*E(t)^2 + delta_NE*N(t)*P(t) - delta_EL*E(t) + E(t)*P(t)*rho_E,
y1(t) = N(t),
y3(t) = M(t),
y2(t) = S(t) + E(t)
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


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

s1 := 0.6:
s2 := 2.6:

(* The equations to solve in order to obtain the coeficients cc are
  G(s1) = 0
  G(s2) = 1
 G'(s1) = 0
 G'(s2) = 0
G''(s1) = 0
G''(s2) = 0
*)

cc0 :=  s1^3*(s1^2 - 5*s1*s2 + 10*s2^2)/(s1 - s2)^5:
cc1 := -30*s1^2*s2^2/(s1 - s2)^5:
cc2 :=  30*s1*s2*(s1 + s2)/(s1 - s2)^5:
cc3 := -10*(s1^2 + 4*s1*s2 + s2^2)/(s1 - s2)^5:
cc4 :=  15*(s1 + s2)/(s1 - s2)^5:
cc5 := -6/(s1 - s2)^5:

$define gga_x_rpbe_params
$include "gga_x_rpbe.mpl"
$include "gga_x_wc.mpl"

g := s -> cc0 + cc1*s + cc2*s^2 + cc3*s^3 + cc4*s^4 + cc5*s^5:

f0 := s -> convert(piecewise(
   s < s1, f0_wc(s),
   s > s2, f0_rpbe(s),
   g(s)*f0_rpbe(s) + (1 - g(s))*f0_wc(s)
), 'Heaviside'):

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


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

# Not that the files have to be included in this specific order
$define gga_x_pw91_params
$include "gga_x_pw91.mpl"

$include "op.mpl"

op_qab         := 2.3706:
op_enhancement := xs -> pw91_f(xs):


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

cases := [
    [
        [[a / b, b], [a, b]], true
    ],
    [
        [[a * b, a + b], [a, b]], false
    ],
    [
        [[a + b^2 - c^2, b * c, b / c], [a, b^2, c^2, b / c]], true
    ]
];

num_passed := 0:
num_failed := 0:

for case in cases do
    input := case[1]:
    correct := case[2]:
    if CompareFields(op(input)) = correct then
        printf("PASSED");
        num_passed := num_passed + 1:
    else
        printf("FAILED");
        num_failed := num_failed + 1:
    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("est_pi.mpl"):
printf("\nDoing calculation using loop:\n\n"):
est_pi(n):
printf("\nDoing calculation using sum:\n\n"):
native_est_pi(n):


text:
input := FileTools:-Text:-ReadFile("AoC-2021-5-input.txt" ):
lines:=map(l->[[parse(l[1])],[parse(l[-1])]],
    {op(map(StringTools:-Split, StringTools:-Split(input,"\n"), "->"))}):

pointcounts1 := table(sparse=0):
pointcounts2 := table(sparse=0):

for l in lines do
    step := signum~(l[2] - l[1]);
    ls := l[1];
    pointcounts2[op(ls)] += 1;
    if 0 in step then
        pointcounts1[op(ls)] += 1;
    end if;
    while ls <> l[2] do
        ls := ls + step;
        pointcounts2[op(ls)] += 1;
        if 0 in step then
            pointcounts1[op(ls)] += 1;
        end if;
    end do;
end do:

answer1:=nops(select(i->i>1,[entries(pointcounts1, nolist)]));
answer2:=nops(select(i->i>1,[entries(pointcounts2, nolist)]));

# Visualization

# just the lines
plots:-display(seq(plottools:-line(l[], thickness=0.01),l in lines));

# fancy hotspots
rng := [min,max](entries(pointcounts2,nolist));
colors := [ColorTools:-Color("CVD 8"),ColorTools:-Gradient("CVD 2".."CVD 9", number=rng[2]-rng[1]-2)[]];
epc := [seq(Array(map([lhs], select(p->rhs(p)=i,[entries(pointcounts2,pairs)]))),i=rng[1]..rng[2])]:

plots:-display(
   seq(plottools:-point(epc[i], symbol=solidcircle, color=colors[i], symbolsize=i),i=rng[1]..rng[2] ),
   size=[5000,5000], axes=none, background="Black");


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

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

(* type: lda_exc *)
(* prefix:
  lda_c_chachiyo_params *params;

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

e0 := rs -> params_a_ap*log(1 + params_a_bp/rs + params_a_bp/rs^2):
e1 := rs -> params_a_af*log(1 + params_a_bf/rs + params_a_bf/rs^2):

f_chachiyo := (rs, zeta) -> e0(rs) + (e1(rs) - e0(rs))*f_zeta(zeta):
f := (rs, zeta) -> f_chachiyo(rs, zeta):


text:
flavor 1
srclang 3
id 65535
numfuncs 7
import "InterfaceTest.mplt"
import "/home/bravewtz/Desktop/openarkcompiler/libjava-core/libjava-core.mplt"
entryfunc &LInterfaceTest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V
fileinfo {
  @INFO_filename "InterfaceTest.jar"}
srcfileinfo {
  1 "InterfaceTest.java"}
javaclass $LBase_3B <$LBase_3B>
javaclass $LDerived_3B <$LDerived_3B>
javainterface $LInter_3B <$LInter_3B> abstract
javaclass $LInterfaceTest_3B <$LInterfaceTest_3B> public
func &LBase_3B_7C_3Cinit_3E_7C_28_29V constructor (var %_this <* <$LBase_3B>>) void
func &LBase_3B_7Cfoo_7C_28_29V public virtual (var %_this <* <$LBase_3B>>) void
func &LDerived_3B_7C_3Cinit_3E_7C_28_29V constructor (var %_this <* <$LDerived_3B>>) void
func &LDerived_3B_7Cfoo_7C_28_29V public virtual (var %_this <* <$LDerived_3B>>) void
func &LInter_3B_7Cfoo_7C_28_29V public virtual (var %_this <* <$LInter_3B>>) void
func &LInterfaceTest_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$LInterfaceTest_3B>>) void
func &LInterfaceTest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg5_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void
var $_C_STR_477e2dc723978d7245344fc51fc13e2b fstatic <[4] u64> readonly = [0, 0x1500000000, 0x65736142062a270a, 0x29286f6f662e]
var $_C_STR_d820ddedd204b763a52374e860ce318e fstatic <[5] u64> readonly = [0, 0x1b00000000, 0x69726544357c6116, 0x286f6f662e646576, 41]
var $_C_STR_a3cc511f7ef4b5079c824f99850ea7ef fstatic <[4] u64> readonly = [0, 0x1700000000, 0x65746e497a755c95, 0x29286f6f662e72]
var $Ljava_2Flang_2FSystem_3B_7Cout extern <* <$Ljava_2Fio_2FPrintStream_3B>> final public static
var $__cinf_Ljava_2Flang_2FString_3B <$__class_meta__>
func &MCC_GetOrInsertLiteral () <* <$Ljava_2Flang_2FString_3B>>
func &LBase_3B_7C_3Cinit_3E_7C_28_29V constructor (var %_this <* <$LBase_3B>>) void {
  funcid 48153
  var %Reg1_R43697 <* <$LBase_3B>>
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>>

  dassign %Reg1_R43697 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43697))
  callassigned &Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  return ()
}
func &LBase_3B_7Cfoo_7C_28_29V public virtual (var %_this <* <$LBase_3B>>) void {
  funcid 48154
  var %Reg2_R43697 <* <$LBase_3B>>
  var %Reg0_R562 <* <$Ljava_2Fio_2FPrintStream_3B>>
  var %Reg1_R43 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161333 <* <$Ljava_2Flang_2FString_3B>>

  dassign %Reg2_R43697 0 (dread ref %_this)
  #INSTIDX : 0||0000:  getstatic
  intrinsiccallwithtype <$Ljava_2Flang_2FSystem_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R562 0 (dread ref $Ljava_2Flang_2FSystem_3B_7Cout)
  #INSTIDX : 3||0003:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_477e2dc723978d7245344fc51fc13e2b) { dassign %L_STR_161333 0 }
  dassign %Reg1_R43 0 (dread ptr %L_STR_161333)
  #INSTIDX : 5||0005:  invokevirtual
  virtualcallassigned &Ljava_2Fio_2FPrintStream_3B_7Cprintln_7C_28Ljava_2Flang_2FString_3B_29V (dread ref %Reg0_R562, dread ref %Reg1_R43) {}
  #INSTIDX : 8||0008:  return
  return ()
}
func &LDerived_3B_7C_3Cinit_3E_7C_28_29V constructor (var %_this <* <$LDerived_3B>>) void {
  funcid 48155
  var %Reg1_R43699 <* <$LDerived_3B>>
  var %Reg1_R43697 <* <$LBase_3B>>

  dassign %Reg1_R43699 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  dassign %Reg1_R43697 0 (retype ref <* <$LBase_3B>> (dread ref %Reg1_R43699))
  callassigned &LBase_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg1_R43697) {}
  #INSTIDX : 4||0004:  return
  return ()
}
func &LDerived_3B_7Cfoo_7C_28_29V public virtual (var %_this <* <$LDerived_3B>>) void {
  funcid 48156
  var %Reg2_R43699 <* <$LDerived_3B>>
  var %Reg0_R562 <* <$Ljava_2Fio_2FPrintStream_3B>>
  var %Reg1_R43 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161344 <* <$Ljava_2Flang_2FString_3B>>

  dassign %Reg2_R43699 0 (dread ref %_this)
  #INSTIDX : 0||0000:  getstatic
  intrinsiccallwithtype <$Ljava_2Flang_2FSystem_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R562 0 (dread ref $Ljava_2Flang_2FSystem_3B_7Cout)
  #INSTIDX : 3||0003:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_d820ddedd204b763a52374e860ce318e) { dassign %L_STR_161344 0 }
  dassign %Reg1_R43 0 (dread ptr %L_STR_161344)
  #INSTIDX : 5||0005:  invokevirtual
  virtualcallassigned &Ljava_2Fio_2FPrintStream_3B_7Cprintln_7C_28Ljava_2Flang_2FString_3B_29V (dread ref %Reg0_R562, dread ref %Reg1_R43) {}
  #INSTIDX : 8||0008:  return
  return ()
}
func &LInter_3B_7Cfoo_7C_28_29V public virtual (var %_this <* <$LInter_3B>>) void {
  funcid 48157
  var %Reg2_R43701 <* <$LInter_3B>>
  var %Reg0_R562 <* <$Ljava_2Fio_2FPrintStream_3B>>
  var %Reg1_R43 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161347 <* <$Ljava_2Flang_2FString_3B>>

  dassign %Reg2_R43701 0 (dread ref %_this)
  #INSTIDX : 0||0000:  getstatic
  intrinsiccallwithtype <$Ljava_2Flang_2FSystem_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R562 0 (dread ref $Ljava_2Flang_2FSystem_3B_7Cout)
  #INSTIDX : 3||0003:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_a3cc511f7ef4b5079c824f99850ea7ef) { dassign %L_STR_161347 0 }
  dassign %Reg1_R43 0 (dread ptr %L_STR_161347)
  #INSTIDX : 5||0005:  invokevirtual
  virtualcallassigned &Ljava_2Fio_2FPrintStream_3B_7Cprintln_7C_28Ljava_2Flang_2FString_3B_29V (dread ref %Reg0_R562, dread ref %Reg1_R43) {}
  #INSTIDX : 8||0008:  return
  return ()
}
func &LInterfaceTest_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$LInterfaceTest_3B>>) void {
  funcid 48158
  var %Reg1_R43703 <* <$LInterfaceTest_3B>>
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>>

  dassign %Reg1_R43703 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43703))
  callassigned &Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  return ()
}
func &LInterfaceTest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg5_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void {
  funcid 48159
  var %Reg0_R43699 <* <$LDerived_3B>>
  var %Reg2_R43699 <* <$LDerived_3B>>
  var %Reg3_R43699 <* <$LDerived_3B>>
  var %Reg3_R43697 <* <$LBase_3B>>
  var %Reg4_R43699 <* <$LDerived_3B>>
  var %Reg4_R43701 <* <$LInter_3B>>

  intrinsiccallwithtype <$LInterfaceTest_3B> JAVA_CLINIT_CHECK ()
  #INSTIDX : 0||0000:  new
  intrinsiccallwithtype <$LDerived_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R43699 0 (gcmalloc ref <$LDerived_3B>)
  #INSTIDX : 3||0003:  dup
  #INSTIDX : 4||0004:  invokespecial
  callassigned &LDerived_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg0_R43699) {}
  #INSTIDX : 7||0007:  astore_1
  dassign %Reg2_R43699 0 (dread ref %Reg0_R43699)
  #INSTIDX : 8||0008:  aload_1
  #INSTIDX : 9||0009:  invokevirtual
  virtualcallassigned &LDerived_3B_7Cfoo_7C_28_29V (dread ref %Reg2_R43699) {}
  #INSTIDX : 12||000c:  new
  intrinsiccallwithtype <$LDerived_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R43699 0 (gcmalloc ref <$LDerived_3B>)
  #INSTIDX : 15||000f:  dup
  #INSTIDX : 16||0010:  invokespecial
  callassigned &LDerived_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg0_R43699) {}
  #INSTIDX : 19||0013:  astore_2
  dassign %Reg3_R43699 0 (dread ref %Reg0_R43699)
  #INSTIDX : 20||0014:  aload_2
  #INSTIDX : 21||0015:  invokevirtual
  dassign %Reg3_R43697 0 (retype ref <* <$LBase_3B>> (dread ref %Reg3_R43699))
  virtualcallassigned &LBase_3B_7Cfoo_7C_28_29V (dread ref %Reg3_R43697) {}
  #INSTIDX : 24||0018:  new
  intrinsiccallwithtype <$LDerived_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R43699 0 (gcmalloc ref <$LDerived_3B>)
  #INSTIDX : 27||001b:  dup
  #INSTIDX : 28||001c:  invokespecial
  callassigned &LDerived_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg0_R43699) {}
  #INSTIDX : 31||001f:  astore_3
  dassign %Reg4_R43699 0 (dread ref %Reg0_R43699)
  #INSTIDX : 32||0020:  aload_3
  #INSTIDX : 33||0021:  invokeinterface
  dassign %Reg4_R43701 0 (retype ref <* <$LInter_3B>> (dread ref %Reg4_R43699))
  interfacecallassigned &LInter_3B_7Cfoo_7C_28_29V (dread ref %Reg4_R43701) {}
  #INSTIDX : 38||0026:  return
  return ()
}

max_stars_repo_name: Bravewtz/OpenArkCompiler

max_issues_repo_name: Bravewtz/OpenArkCompiler

max_forks_repo_name: Bravewtz/OpenArkCompiler


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

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

(* type: mgga_exc *)

$include "mgga_x_tm.mpl"
$include "lda_x_erf.mpl"

(* See text after eq 7  *)
js18_H := x -> (2*tm_lambda - 1)^2 * tm_p(x):
(* See text after eq 7 *)
js18_G := (x, t) ->
  (3*(tm_lambda^2 - tm_lambda + 1/2)*(t - K_FACTOR_C - x^2/72) - (t - K_FACTOR_C)
   + 7/18*(2*tm_lambda - 1)^2*x^2)/K_FACTOR_C:

(* not to run into problems with erfs and exp *)
js18_A := (rs, z, x) -> m_max(1e-10, a_cnst*rs/(tm_f0(x)*opz_pow_n(z,1/3))):

js18_DME_SR := (rs, z, x, t) ->
  + attenuation_erf   (js18_A(rs, z, x))/tm_f0(x)^2
  + attenuation_erf_f2(js18_A(rs, z, x))*7*js18_G(x, t)/(9*tm_f0(x)^4)
  + attenuation_erf_f3(js18_A(rs, z, x))*245*js18_H(x)/(54*tm_f0(x)^4):

pjs18_f := (rs, z, x, u, t) -> js18_DME_SR(rs, z, x, t):

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


text:
with(Groebner):
F:={2*x1 + 5*x2 + x3,2*x1 + x2 + x3,9*x1 + 6*x2 + x3,x1^2 + x2^2 + x3^2 - 1}:
m_ord:=tdeg(x1,x2,x3):
G:=Basis(F,m_ord):
NormalSet(G,m_ord)[1]

max_stars_repo_name: marcusvaltonen/alggeompoly

max_issues_repo_name: marcusvaltonen/alggeompoly

max_forks_repo_name: marcusvaltonen/alggeompoly


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

# A standard Stasheff tree for n is a Stasheff tree for the set
# {1,...,n} with its standard ordering.

`is_element/standard_stasheff_trees` := (n::posint) -> proc(TT)
 local A,T,i;
 global reason;

 A := {seq(i,i=1..n)};

 if not(`is_element/full_trees`(A)(TT)) then
  reason := ["is_element/standard_stasheff_trees","TT is not a full tree",TT,reason];
  return false;
 fi;

 for T in TT do
  if nops(T) <> max(T) - min(T) + 1 then
   reason := ["is_element/standard_stasheff_trees","T is not an interval",T];
   return false;
  fi;
 od;

 return true;
end;

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

`is_equal/standard_stasheff_trees` := (n::posint) -> proc(TT,UU)
 local A,i;

 A := {seq(i,i=1..n)};
 return `is_equal/trees`(A)(TT,UU);
end;

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

`is_leq/standard_stasheff_trees` := (n::posint) -> proc(TT,UU)
 local A,i;

 A := {seq(i,i=1..n)};
 return `is_leq/trees`(A)(TT,UU);
end;

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

`random_element/standard_stasheff_trees` := (n::posint) -> proc()
 local N,pi,TT,UU,U,i;

 if n = 0 then return FAIL; fi;
 if n = 1 then return {{1}}; fi;

 N := {seq(i,i=1..n)};

 pi := {N};
 while nops(pi) = 1 do
  pi := `random_element/interval_partitions`(n)();
 od;

 TT := N;

 for U in pi do
  UU := `random_element/standard_stasheff_trees`(nops(U))();

  UU := map(V -> map(v -> v+min(U)-1,V),UU);
  TT := TT,op(UU);
 od;

 TT := {TT};

 return TT;
end;

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

`list_elements/standard_stasheff_trees` := proc(n::posint)
 option remember;

 local TTT,UUU,TT,UU,N,IP,pi,pi_list,P,B,b,i;
 if n = 0 then return []; fi;
 if n = 1 then return [{{1}}]; fi;

 N := {seq(i,i=1..n)};
 TTT := NULL;
 IP := select(pi -> pi <> {N},`list_elements/interval_partitions`(n));

 for pi in IP do
  pi_list := [op(pi)];
  P := NULL;
  for B in pi_list do
   b := min(op(B));
   UUU := `list_elements/standard_stasheff_trees`(nops(B));
   UUU := map(UU -> map(V -> map(v -> b-1+v,V),UU),UUU);
   P := P,UUU;
  od;
  P := `cartesian_product/lists`(P);
  TTT := TTT,op(map(UU -> {N,op(map(op,UU))},P));
 od:
 return([TTT]);
end:

######################################################################
# This is A001003 in OEIS

`count_elements/standard_stasheff_trees` := proc(n)
 local f,L,i;
 if n = 0 then
  return 0;
 elif n = 1 then
  return 1;
 else
  f := (k) -> op([seq(i,i=0..k+1),seq(k-i,i=0..k)]);
  L := [0];
  for i from 1 to n-2 do
   L := map(f,L);
  od;
  return nops(L);
 fi;
end;

# Alternative formula for n > 1:
# add(2^k*3^(n-2-2*k)*binomial(n-2,2*k)*binomial(2*k,k)/(k+1),k=0..floor((n-2)/2));

`dim/standard_stasheff_trees` := (n) -> proc(TT)
 2*n - 1 - nops(TT);
end:

`count_by_dim/standard_stasheff_trees` := proc(n,k)
 if k < 0 or k > n-2 then return 0; fi;

 return binomial(n-2,n-k-2) * binomial(2*n-k-2,n) / (n - k - 1);
end:

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

`draw/standard_stasheff_trees` := (n::posint) -> proc(TT)
 local p,P,T,U,m,i;
 p := table();
 P := NULL;
 for T in TT do
  p[T] := [`+`(op(T))/nops(T),nops(T)];
 od;

 m := parent_map({seq(i,i=1..n)})(TT);

 for T in TT do
  if nops(T) < n then
   U := m[T];
   P := P,line(p[T],p[U]);
  fi;
 od;
 for T in TT do
  P := P,point(p[T],colour=red,symbol=solidcircle,symbolsize=20);
 od;

 P := display(P,axes=none,scaling=constrained);
 return P;
end:

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

`root_children/standard_stasheff_tree` := (n::posint) -> proc(TT)
 local C,A,i,XX,X,m;

 if n = 1 then return []; fi;

 C := NULL;
 A := {seq(i,i=1..n)};
 while A <> {} do
  i := min(op(A));
  XX := select(X -> member(i,X) and nops(X) < n,TT);
  m := max(op(map(nops,XX)));
  XX := select(X -> nops(X) = m,XX);
  X := XX[1];
  C := C,X;
  A := A minus X;
 od:
 return [C];
end:

######################################################################
# A Stasheff tree TT on A is binary if every non-minimal set in TT has
# precisely two children.  This holds iff |TT| = 2|A|-1.

`is_binary/standard_stasheff_trees` := (n::posint) -> proc(TT)
 return evalb(nops(TT) = 2*n-1);
end:

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

`e/loday` := (n::posint) -> proc(pq)
 local p,q;
 p,q := op(pq);
 return [0 $ (p-1), 1 $ (q-p), 0 $ (n-q)];
end:

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

`m/loday` := (n::posint) -> pq -> (pq[2]-pq[1]+1)*(pq[2]-pq[1])/2;

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

`p/loday` := (n::posint) -> proc(pq,u)
 local i;
 return add(u[i],i=pq[1]..pq[2]-1) - `m/loday`(n)(pq);
end:

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

`is_element/PK` := (n::posint) -> proc(u)
 global reason;
 local p,q;

 if not `is_element/R`(n-1)(u) then
  reason := [convert(procname,string),"u is not a list of length n-1",u,n-1];
  return false;
 fi;

 if `p/loday`(n)([1,n],u) <> 0 then
  reason := [convert(procname,string),"p_A(u) <> 0",u];
  return false;
 fi;

 for p from 1 to n do
  for q from p to n do
   if `p/loday`(n)([p,q],u) < 0 then
    reason := [convert(procname,string),"p_J(u) < 0",u,[p,q]];
    return false;
   fi;
  od;
 od;

 return true;
end;

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

`tau/loday` := (n::posint) -> proc(u)
 local TT,p,q,i;
 TT := NULL;
 for p from 1 to n do
  for q from p to n do
   if `p/loday`(n)([p,q],u) = 0 then
    TT := TT,{seq(i,i=p..q)};
   fi;
  od;
 od;
 TT := {TT};

 return TT;
end;

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

`is_element/Sigma0_loday` := (n::posint) -> proc(sigma)
 global reason;
 local JJ,J,g,i,j,k;

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

 JJ := {seq(seq({seq(k,k=i..j)},j=i+1..n),i=1..n)};
 if map(op,{indices(sigma)}) <> JJ then
  reason := [convert(procname,string),"sigma is not indexed by the appropriate set of intervals",sigma,JJ];
  return false;
 fi;

 for J in JJ do
  g := sigma[J];
  if not (type(g,set) and nops(g) = 2 and max(op(g)) - min(op(g)) = 1) then
   reason := [convert(procname,string),"sigma(J) is not a gap",sigma,J,g];
   return false;
  fi;

  if g minus J <> {} then
   reason := [convert(procname,string),"sigma(J) is not a gap in J",sigma,J,g];
   return false;
  fi;
 od;

 return true;
end;

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

`w0/loday` := (n::posint) -> proc(sigma)
 local v,JJ,J,i,j,k;

 v := Vector(n-1);
 JJ := {seq(seq({seq(k,k=i..j)},j=i+1..n),i=1..n)};
 for J in JJ do
  i := min(op(sigma[J]));
  v[i] := v[i] + 1;
 end:

 return convert(v,list);
end;

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

`pi/loday` := (n::posint) -> (TT) -> proc(J)
 local UU,k;

 UU := select(T -> J minus T = {},TT);
 k := min(map(nops,UU));
 UU := select(U -> nops(U) = k,UU);
 return(UU[1]);
end;

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

`w/loday` := (n::posint) -> proc(TT)
 local v,JJ,J,P,G,T,S,g,m,i,j,k;

 v := Vector(n-1);
 JJ := {seq(seq({seq(k,k=i..j)},j=i+1..n),i=1..n)};
 P := table();
 for J in JJ do
  P[J] := `pi/loday`(n)(TT)(J);
 od;
 G := {seq({i,i+1},i=1..n-1)};
 for J in JJ do
  T := P[J];
  S[J] := NULL;
  for j from min(J) to max(J)-1 do
   g := {j,j+1};
   if P[g] = T then
    S[J] := S[J],g;
   fi;
  od;
  S[J] := {S[J]};
  m := nops(S[J]);
  for g in S[J] do
   v[min(g)] := v[min(g)] + 1/m;
  od;
 od;
 return convert(v,list);
end:

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

`w/loday/bin` := (n::posint) -> proc(TT)
 local C,T,v,i;
 C := children_map({seq(i,i=1..n)})(TT);
 v := Vector(n-1);
 for i from 1 to n-1 do
  T := `pi/loday`(n)(TT)({i,i+1});
  v[i] := `*`(op(map(nops,C[T])));
 od;
 return convert(v,list);
end:

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

`T/loday` := (n::posint) -> proc(i)
 local j,k;
 {seq({j},j=1..n),{seq(j,j=1..n)},
  seq({seq(k,k=1..j)},j=2..i),
  seq({seq(k,k=j..n)},j=i+1..n-1)
 };
end:

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

`phi/standard_stasheff_trees/PK` := (n::posint) -> (TT) -> `w/loday`(n)(TT);

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

`phi/realisation/standard_stasheff_trees/PK` := (n::posint) -> proc(x)
 local C,v,TT;

 C := map(op,[indices(x)]);
 v := [0$(n-1)];
 for TT in C do
  v := v +~ x[TT] *~ `phi/standard_stasheff_trees/PK`(n)(TT);
 od;

 return v;
end:


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

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

(* type: mgga_exc *)

$include "mgga_x_tm.mpl"

(* revtm uses the tpss definition of qtilde *)
tm_b := 0.4:
tm_qtilde := (x, t) ->
     9/20 * (tm_alpha(x, t) - 1)/sqrt(1 + tm_b*tm_alpha(x, t)*(tm_alpha(x, t) - 1))
     + 2*tm_p(x)/3:


text:
# Note this is valid only because this file is technically a top-level file,
# i.e. it isn't included anywhere with `$include'
with(Domain_Type):

local Simplify_DConstrain := (can_simp, do_simp) -> module()
    uses Domain;

    export ModuleApply := proc(vs0 :: DomBound, sh :: DomShape, $)
        local vs := vs0, ctx;
        if _Env_HakaruSolve=false then
           return sh end if;
        vs  := Domain:-Bound:-varsOf(vs,"set");
        ctx := Domain:-Bound:-toConstraints(vs0, 'bound_types');
        subsindets(sh, DomConstrain
                  ,x->do_simpl_constraints(vs0, vs, ctx, x));
    end proc;

    local do_simpl_constraints := proc(vs0, vars, ctx_vs, x, $)
        local ctx1, ctx, ss, td, rest, d, in_vs;
        ss, ctx := selectremove(q->depends(q,indets(vars,Name)), x);
        td, rest := selectremove(can_simp(vs0), ss);
        ctx1 := { op(ctx), op(ctx_vs), op(rest) };
        d := map(x->do_simp(vs0,ctx1,x), td);
        DConstrain(op(d), op(ctx), op(rest));
    end proc;
end module;

Domain_try_improve_exp := module()
  uses Utilities;

  export SimplName  := "Try improve exp";
  export SimplOrder := 6 - 2/10;

  local can := proc(vs)
    c -> c::relation and (not(lhs(c) :: Name) and not(rhs(c) :: Name) and has(c,{ln,exp}));
  end proc;

  local `try` := proc(dbnd, ctx1, q0, $)
    local vars, vars_q, q := q0;
    vars := Domain:-Bound:-varsOf(dbnd,"set");
    vars_q := indets(q, name) intersect vars;
    if nops(vars_q)=0 then return q end if;

    vars_q := op(1,vars_q);
    try_improve_exp(q, vars_q, ctx1);
  end proc;

  export ModuleApply := Simplify_DConstrain(can, `try`);
end module;

constraints_about_vars := module()
  uses Utilities;
  export SimplName  := "Make constraints abouts vars";
  export SimplOrder := 6;

  local can_make_about := proc(vs)
    local vars := Domain:-Bound:-varsOf(vs,"set");
    c -> c::relation and (not(lhs(c) in vars) and not(rhs(c) in vars));
  end proc;
  local try_make_about := proc(dbnd, ctx1_, q0, $)
      local vars_q, vars, q_s, q_r
          , q := q0
          , ctx1 := ctx1_ ;
      vars := Domain:-Bound:-varsOf(dbnd,"set");
      vars_q := indets(q, name) intersect vars;
      if nops(vars_q)=0 then return q end if;

      vars_q := op(1,vars_q);
      q := try_improve_exp(q, vars_q, ctx1);
      q_r := `if`(has(q,{ln,exp}),q0,q);

      q, ctx1 := map(KB:-chill, [q,ctx1])[];
      # Unevaluated for debugging
      q_s := 'solve({q},[op(vars_q)], 'useassumptions'=true) assuming (op(ctx1))';
      q_s := eval(q_s);
      q, ctx1 := map(KB:-warm, [q,ctx1])[];
      if q_s::list and not has(q_s, signum) then
        if nops(q_s)=0 then return q_r end if;
        q_s := map(s->remove(c->c in ctx1 or `and`(c::relation,lhs(c)::name,lhs(c)=rhs(c)), s), q_s);
        q_s := remove(x->x=[],q_s);
        if nops(q_s)=0 then return q_r end if;
        userinfo(3, 'constraints_about_vars'
                ,printf("Found a solution to %a (with solve):\n\t%a\n", q0, op(1,q_s)));
        op(op(1,q_s)); # pick the first solution arbitrarily!
      else
        q_r
      end if;
  end proc;

  export ModuleApply := Simplify_DConstrain(can_make_about, try_make_about);
end module;

clamp_extraneous_constraints := module()
  uses Utilities;
  export SimplName  := "clamp_extraneous_constraints";
  export SimplOrder := 6 - 1/10;

  local can := proc(vs)
    local vars := Domain:-Bound:-varsOf(vs,"set");
    c-> c::relation and
        (ormap(s->s(c)::Name and not(c in vars),[lhs,rhs]) and has(c, vars));
  end proc;

  local `try` := proc(dbnd, ctx1, q0, $)
    local bnd, b_ty, b_rel, b_var, b_bnd, b_bnd_vs, vars, extremum, q := q0;
    if has(q, {ln,exp}) then
      vars := Domain:-Bound:-varsOf(dbnd, "set");
      bnd := classify_relation(q0, x -> type(x, Name) and not (x in vars));
      if bnd=FAIL then return q0 end if;
      b_ty, b_rel, b_var, b_bnd := op(bnd):
      b_bnd_vs := indets(b_bnd, satisfies(x -> x in vars));
      if b_bnd_vs <> {} then
        b_ty := `if`(b_ty=B_LO, 'minimize', 'maximize');
        extremum := b_ty(b_bnd, op(b_bnd_vs)) assuming(op(ctx1));
        if not (ext :: SymbolicInfinity or has(extremum, {maximize,minimize})) then
          q := q, b_rel(b_var,extremum);
        end if;
      end if;
    end if;
    q
  end proc;
  export ModuleApply := Simplify_DConstrain(can, `try`);
end module;

# Pushes constraints down, or pulls them up, when there are such constraints.
local do_ctx_dir := dir -> proc(vs :: DomBound, sh :: DomShape, $)
    local ctx := Domain:-Bound:-toConstraints(vs, 'no_infinity');
    ctx := remove(x->x::`::`,ctx);
    if ctx <> {} then
        dir(ctx)(sh);
    else
        sh;
    end if;
end proc;

push_ctx_down := module()
  uses Domain;
  export ModuleApply := do_ctx_dir(ctx->sh->subsindets(sh, DomConstrain, x->if nops(x)<>0 then DConstrain(op(x), op(ctx)) else DConstrain() end if));
  export SimplName  := "Push context down";
  export SimplOrder := 1;
end module;

pull_ctx_out := module()
  export ModuleApply :=
    do_ctx_dir(ctx->sh->subsindets(sh, DomConstrain
              ,x->remove(c->c in ctx or (is(c) assuming op(ctx)), x)));
  export SimplName  := "Pull context out";
  export SimplOrder := 20;
end module;

# Turns constraints(relations) into bounds
classify_DConstrains := module()
  uses Domain, Utilities;

  export SimplName  := "Classify constraints";
  export SimplOrder := 8;

  export ModuleApply := proc(dbnds::DomBound,dsh::DomShape,$)
    subsindets(dsh, DomConstrain, x->classify_Rels([op(x)], dbnds));
  end proc;

  # Split relations into those we can incorporate
  # into bounds and others; this process also includes
  # some intermediate computations about the new bound:
  #  [ DInto(var), orig, upd_bound, bound_kind, numVars ]
  #    var         :=   variable to be modified
  #    orig        :=   the original relation
  #    upd_bound   :=   how to modify (a function)
  #    bound_kind  :=   {LO/HI}
  #    numVars     :=   how many variables are mentioned by the bound
  local classify_Rel := proc(r0,dbnd,$)
    local ret, vars, mk_k, bnd_k, i, var_k, q, cl, lhs_r, rhs_r, r := r0;
    ret := [ DConstrain, r0 ];
    vars := Domain:-Bound:-varsOf(dbnd,"set");

    # extract a representation of the relation, and check that it is a bound
    cl := classify_relation(r, vars);
    if cl=FAIL or not (op(1,cl) in {B_LO,B_HI}) then return ret; end if;
    bnd_k,mk_k,lhs_r,rhs_r := op(cl);

    i := Domain:-Bound:-varIx_mb(dbnd, lhs_r);
    if i >= 0 then
      var_k := op([1,i,3], dbnd);
      q := Domain:-ExtBound[var_k]:-RecogBound(mk_k,rhs_r);
      if q <> NULL then
        ret := [ DInto(lhs_r), r0, q, bnd_k, countVs(vars minus {lhs_r},r0) ];
      end if;
    end if;
    ret;
  end proc;

  # Classifies a set of relations as potential bounds. The process
  # is as follows:
  #  - Classify each relation individually (see classify_Rel)
  #  - "Categorize" by the classification
  #  - Split relations which will be DConstrains, and those where we found
  #      more than one LO/HI bound or more than 2 bounds total.
  #  - Sort the relations which are to become bounds by the number of variables
  #      mentioned in the bound
  #  - Rebuild the continuations producing each new bound (see make_DInto)
  #  - Rebuild the leftover relations as a DConstrain.
  #  - Squash together ("foldr (.) id fs z") the continuations and leftover relations.
  local classify_Rels := proc(rs0,dbnd,$)
    local cs1, cs2, cs, rs := rs0;
    rs := map(r->classify_Rel(r,dbnd), rs);
    rs := [ListTools[Categorize]( ((a,b) -> op(1,a)=op(1,b) ), rs )];
    cs1, rs := selectremove(x->op([1,1],x)=DConstrain, rs);
    rs, cs2 := selectremove(x->nops(x)=1 or (nops(x)=2 and op([1,4],x)<>op([2,4],x)), rs);
    rs := sort(rs, key=solOrder);
    rs := map(r->make_DInto(r,dbnd),rs);
    cs := DConstrain(map(c->map(q->op(2,q),c)[], [op(cs1),op(cs2)])[]);
    foldr((f,g)->f@g,x->x,op(rs))(cs);
  end proc;

  # Rebuild a classified (by classify_Rel) relation; this
  # takes a list of such classifications, all of whose variable should
  # be identical, and for which the different bound kinds are all distinct.
  # None of this is checked!
  # (since there are 2 bound kinds, this list is at most length 2)
  local make_DInto := proc(r::list([specfunc(DInto),anything,anything,identical(B_LO,B_HI),nonnegint]), dbnd, $)
    local var_nm, var_ty;
    var_nm := op([1,1,1], r);
    var_ty := Domain:-Bound:-get(dbnd, var_nm)[1];
    var_ty := foldr( ((q,x)->op(3,q)(x)), var_ty, op(r) );
    (ctx -> DInto(var_nm, var_ty, ctx))
  end proc;
  local flip_rel := table([`=`=`=`,`<>`=`<>`,`<=`=`>=`,`<`=`>`,`>=`=`<=`,`>`=`<`]);
  local countVs := ((vs,c)-> nops(indets(c, DomBoundVar) intersect {op(vs)} ));
  local solOrder := (x->`+`(map(z->op(5,z),x)[]));
end module;

# todo; this should actually solve for a variable, then substitute
# that variable in. In most cases, it would probably be enough to
# leave that as it is; it would simplify later.
singular_pts := module()
  uses Domain, Domain_Type, Hakaru;

  export SimplName  := "Single_pts";
  export SimplOrder := 7;

  export ModuleApply := proc(bnds :: DomBound, sh :: DomShape, $)
    local kb, kb_rn, ns, sh1; sh1 := sh;
    kb, kb_rn := Domain:-Bound:-toKB(bnds)[];
    sh1 := subs(kb_rn, sh1);
    ns := op(1,bnds);
    ns := select(b->op(3,b) in {`Int`,`Ints`}, ns);
    ns := map(curry(op,1), ns);

    sh1 := Partition:-Simpl:-singular_pts(sh1, kb, _name_cands=ns);
    eval(sh1, [DSum=DSum,DConstrain=DConstrain]);
  end proc;
end module;

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
               2019 Susi Lehtola

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

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

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

$include "gga_c_pw91.mpl"

optc_f2 := (rs, z, xt, xs0, xs1) ->
  + f_pw91(rs*(2/(1 + z))^(1/3),  1, xs0, xs0, 0)*opz_pow_n( z,1)/2
  + f_pw91(rs*(2/(1 - z))^(1/3), -1, xs1, 0, xs1)*opz_pow_n(-z,1)/2:

f  := (rs, z, xt, xs0, xs1) ->
  + params_a_c1*f_pw91(rs, z, xt, xs0, xs1) + (params_a_c2 - params_a_c1)*optc_f2(rs, z, xt, xs0, xs1):