text:
#@ Not autoload

p := 3;   # the ambient prime
n := 2;   # height
d := 243; # degree of approximate FGL

# find the FGL
T := make_fgl_Morava(p,n,d):
F_E := eval(T["sum"]);
F_K := eval(T["sum0"]);

# Generators and relations for O_Y, where Y is the scheme defined in
# Section sec_Y.  Note that the quantities called a,b and c in the
# main document are aa,bb and cc here.  Note also that we build up
# the relations in two stages: we first do a preliminary version
# without the condition \psi^4(D) = D, then we add that condition
# as a second step.

vars_Y := plex(aa,bb,cc,c[1],c[3],c[2]):
rels_Y := {
 c[1] - (aa+bb+cc),
 c[2] - (aa*bb+bb*cc+cc*aa),
 c[3] - (aa*bb*cc),
 c[1]^9,c[2]^9,c[3]^9,
 aa + F_K(bb,cc)
}:
basis_Y := Basis(rels_Y,vars_Y,characteristic=p); # Preliminary version

# aa4 = [4](aa), bb4 = [4](bb), cc4 = [4](cc)
aa4 := mods(NormalForm(F_K(aa,aa^(p^n)),basis_Y,vars_Y,characteristic=p),p);
bb4 := mods(NormalForm(F_K(bb,bb^(p^n)),basis_Y,vars_Y,characteristic=p),p);
cc4 := mods(NormalForm(F_K(cc,cc^(p^n)),basis_Y,vars_Y,characteristic=p),p);

# To impose the condition psi^4(D) = D, the elementary symmetric functions
# in aa4, bb4, cc4 must be the same as for aa, bb and cc.
rels_Y := {
 op(basis_Y),
 c[1] - (aa4+bb4+cc4),
 c[2] - (aa4*bb4+bb4*cc4+cc4*aa4),
 c[3] - (aa4*bb4*cc4)
}:

basis_Y := Basis(rels_Y,vars_Y,characteristic=p);
NF_Y := (u) -> mods(NormalForm(u,basis_Y,vars_Y,characteristic=p),p):

# The above Groebner basis includes the variables aa, bb and cc which were
# adjoined to O_Y for convenience.  Because these were placed first in a pure
# lexicographic order, we can obtaina Groebner basis for O_Y itself by simply
# omitting any relations that involve aa, bb or cc.
vars_Y1 := plex(c[1],c[3],c[2]):
basis_Y1 := remove(has,basis_Y,{aa,bb,cc});
cobasis_Y1 := cobasis(basis_Y1,vars_Y1);
nops(cobasis_Y1);

# Definitions of some auxiliary elements

# gamma is used by default for the Euler-Mascheroni constant; we need
# to cancel this.
unprotect('gamma'):

alpha := c[3]^3-c[2]^3*c[3];
beta  := 1 + c[2]^4;
gamma := c[2]^3*c[3]^3;
dd := NF_Y(F_K(aa,-bb));
mm := aa^5*bb^5 + bb^5*cc^5 + cc^5*aa^5;
uu := c[2] + c[2]^2*c[3]^2;
for i from 0 to 18 do
 s[i] := NF_Y(aa^i + bb^i + cc^i);
od:

# Formulae in the Claim labelled claim-A
_ASSERT(
 map(NF_Y,{
  c[1]+c[3]^3,
  c[1]^3
 }) = {0},
 "claim-A"
):

# Formulae in the Claim labelled claim-B and its proof
_ASSERT(
 map(NF_Y,{
  c[2]^6,
  s[1]-c[1],
  s[2]-(c[1]^2+c[2]),
  s[3],
  s[4]-c[1]*s[3]+c[2]*s[2]-c[3]*s[1],
  s[12]-s[4]^3,
  aa^27,bb^27,cc^27,
  aa4 - (aa + aa^9),
  bb4 - (bb + bb^9),
  cc4 - (cc + cc^9),
  aa4^4 - (aa^4 + aa^12),
  bb4^4 - (bb^4 + bb^12),
  cc4^4 - (cc^4 + cc^12)
 }) = {0},
 "claim-B"
):

# Formulae in the Claim labelled claim-C and its proof
_ASSERT(
 map(NF_Y,{
  c[1]*c[2]                                   - c[2]^3 * (-c[2]*c[3]),
  c[3]^3*c[2]                                 - c[2]^3 * ( c[2]*c[3]),
  aa^3-c[1]*aa^2+c[2]*aa-c[3],
  aa^9 - c[3]^3                               - c[2]^3 * (-aa^3),
  aa^9 + c[1]                                 - c[2]^3 * (-aa^3),
  aa4 - (aa - c[1])                           - c[2]^3 * (-aa^3),
  aa4*bb4*cc4 - (aa-c[1])*(bb-c[1])*(cc-c[1]) - c[2]^3*(-c[2]*c[3]),
  aa4*bb4*cc4 - (c[3]-c[2]*c[1])              - c[2]^3*(-c[2]*c[3]),
  aa4*bb4*cc4 - c[3]
 }) = {0},
 "claim-C"
):

# Formulae in Corollary cor-C and its proof
_ASSERT(
 map(NF_Y,{
  c[3]^6*c[2]^2,
  c[3]^3*c[2]^4,
  c[3]^3*c[2]-c[2]^3*(c[2]*c[3]),
  c[2]^6
 }) = {0},
 "cor-C"
):

# Formulae in the Claim labelled claim-D and its proof
_ASSERT(
 map(NF_Y,{
  c[3]^3*c[2]-c[3]*c[2]^4,
  aa^3-c[1]*aa^2+c[2]*aa-c[3],
  aa^9+c[2]^3*aa^3-c[3]^3,
  aa^9+c[2]^3*aa^3-c[3]^3,
  aa^9-c[2]^3*c[3]^3*aa^2-c[2]^4*aa+c[2]^3*c[3]-c[3]^3,
  F_K(aa,aa^9) - (aa+aa^9),
  F_K(aa,aa^9) - (alpha + beta*aa + gamma*aa^2),
  alpha*beta - alpha,
  beta*gamma - gamma,
  gamma^2,
  alpha^3,
  alpha^2*gamma,
  beta^3 - 1,
  alpha^2*c[1],
  alpha*gamma*c[1]*c[2],
  gamma*c[1]*c[3],
  aa4*bb4*cc4 -
   ((alpha+beta*aa+gamma*aa^2) *
    (alpha+beta*bb+gamma*bb^2) *
    (alpha+beta*cc+gamma*cc^2)),
  aa4*bb4*cc4 - (alpha^2*c[1]+alpha*c[2]+alpha*gama*c[1]*c[2]+c[3]+gamma*c[1]*c[3]),
  alpha^2*c[1],
  -(c[3]^3-c[2]^3*c[3])^2*c[3]^3,
  alpha*gamma*c[1]*c[2],
  gamma*c[1]*c[3],
  aa4*bb4*cc4 - (c[3]+alpha*c[2]),
  alpha*c[2]
 }) = {0},
 "claim-D"
):

# Formulae in the Claim labelled claim-E and its proof
_ASSERT(
 map(NF_Y,{
  c[3]^6-c[2]^5,
  aa4^2 - (aa^2-aa^10+aa^18),
  aa4^2+bb4^2+cc4^2 - (s[2]-s[10]+s[18]),
  s[18],
  s[10],
  aa^9-(alpha+(beta-1)*aa+gamma*aa^2),
  gamma,
  aa^10 - (alpha*aa + (beta-1)*aa^2),
  s[10] - (alpha*s[1] + (beta-1)*s[2]),
  s[10] - (-c[3]^6+c[2]^5)
 }) = {0},
 "claim-E"
):

# Groebner basis for O_Y
_ASSERT(
 mods(basis_Y1,3) = [c[2]^6,c[3]^3*c[2]-c[3]*c[2]^4,c[3]^6-c[2]^5,c[1]+c[3]^3],
 "prop-Y"
):

# Formulae in the Claim labelled claim-F and its proof
_ASSERT(
 map(NF_Y,{
  dd - (aa-bb+aa^6*bb^3-aa^3*bb^6),
  aa^9+c[2]^3*aa^3-c[3]^3,
  bb^9+c[2]^3*bb^3-c[3]^3,
  (c[2]^3)^2,
  c[2]^3*c[3]^3,
  (c[3]^3)^2-c[2]^5,
  aa^9*bb^9 - c[2]^5,
  aa^3-(-c[3]^3*aa^2-c[2]*aa+c[3]),
  aa^12 * bb^9  - c[2]^5*c[3],
  aa^9  * bb^12 - c[2]^5*c[3],
  aa^9 * bb^9 * (bb^3 - aa^3)
 }) = {0},
 "claim-F"
):

# Formulae in the Claim labelled claim-G and its proof
_ASSERT(
 map(NF_Y,{
  c[2] - (-(aa^2+aa*bb+bb^2)-(aa+bb)*c[3]^3),
  c[3] - (-aa*bb*(aa+bb)-aa*bb*c[3]^3),
  aa+bb - (-c[3]^3-cc),
  aa*bb - (cc^2+c[3]^3*cc+c[2]),
  aa^9*bb^9 - c[2]^5,
  (aa-bb)^2 - (-c[2]-(aa+bb)*c[3]^3),
  aa^3-bb^3 - (-c[2]*(aa-bb)-c[3]^3*(aa^2-bb^2)),
  aa^9-bb^9 - c[2]^4*(aa-bb),
  aa+bb+cc-c[1],
  aa+bb+cc+c[3]^3,
  c[2] - (-(aa^2+aa*bb+bb^2)-(aa+bb)*c[3]^3),
  c[3] - (-aa*bb*(aa+bb)-aa*bb*c[3]^3),
  cc^2 + cc * c[3]^3 + c[2] - aa * bb,
  - aa^18 - aa^9 * bb^9 - bb^18,
  aa^18 + bb^18 + c[2]^5
 }) = {0},
 "claim-G"
):

# Formulae in the Claim labelled claim-H
_ASSERT(
 map(NF_Y,{
  s[ 1] - (-c[3]^3),
  s[ 2] - (c[2]+c[2]^5),
  s[ 3],
  s[ 4] - (-c[3]^4-c[2]^2),
  s[ 5] - (c[3]*c[2]),
  s[ 6] - (c[2]^3),
  s[ 7] - (-c[3]^5+c[3]*c[2]^2),
  s[ 8] - (-c[2]^4+c[3]^2*c[2]),
  seq(s[i],i=9..18)
 }) = {0},
 "claim-H"
);

# Formula in the Claim labelled claim-I
_ASSERT(
 NF_Y(
  dd - (aa-bb)*(1-c[3]^2*c[2]+c[3]^2*c[2]^5-c[2]^4-(c[3]*c[2]^2-c[3]^5)*cc-c[2]^3*cc^2)
 ) = 0,
 "claim-I"
);

# Formulae in the Claim labelled claim-J and its proof
_ASSERT(
 map(NF_Y,{
 dd^3 - (aa-bb)^3,
 dd^3 - (-c[3]^3*(aa^2-bb^2)-c[2]*(aa-bb)),
 c[2]^4*dd - c[2]^4*(aa-bb)*(1-c[2]*c[3]^2),
 c[2]^4*dd^4,
 c[2]^4*(aa-bb)^4,
 c[2]^5*dd - c[2]^5*(aa-bb),
 c[2]^5*dd^2,
 dd-(aa-bb) - aa^3*bb^3*(aa^3-bb^3),
 aa^9*bb^9*(aa^9-bb^9),
 c[2]^4*c[3]^3,
 aa^3-bb^3+c[3]^3*(aa^2-bb^2)+c[2]*(aa-bb),
 c[2]^4*(aa^3-bb^3)+c[2]^5*(aa-bb),
 (c[3]^3)                      * c[2]^5,
 (c[2])                        * c[2]^5,
 (aa*bb-cc^2)                  * c[2]^5,
 ((aa*bb)^3-cc^6)              * c[2]^5,
 (cc^3-c[3])                   * c[2]^5,
 (cc^6-c[3]^2)                 * c[2]^5,
 c[2]^4*(aa^3-bb^3)*(aa*bb)^3 + c[2]^5*c[3]^2*(aa-bb),
 (aa-bb)^2+c[2]+(aa+bb)*c[3]^3
 }) = {0},
 "claim-J"
);

points_E := map(NF_Y,[
 F_K(aa,-bb), F_K(bb,-cc), F_K(cc,-aa),
 F_K(bb,-aa), F_K(cc,-bb), F_K(aa,-cc)
]):
f_E := NF_Y(expand(mul(t - u,u in points_E)));

vars_Z := plex(aa,bb,cc,c[1],z,c[3],c[2]):

az := rem(NF_Y(F_K(aa,z)),z^9,z);
bz := rem(NF_Y(F_K(bb,z)),z^9,z);
cz := rem(NF_Y(F_K(cc,z)),z^9,z);

rels_Z := {
 op(basis_Y),
 z^9,
 c[1] - (az+bz+cz),
 c[2] - (az*bz+bz*cz+cz*az),
 c[3] - (az*bz*cz)
}:
basis_Z := Basis(rels_Z,vars_Z,characteristic=p):
NF_Z := (u) -> mods(NormalForm(u,basis_Z,vars_Z,characteristic=p),p):

vars_Z1 := plex(c[1],z,c[3],c[2]):
basis_Z1 := remove(has,basis_Z,{aa,bb,cc});
cobasis_Z1 := cobasis(basis_Z1,vars_Z1);

# Formulae in the Claim labelled claim-K and its proof
_ASSERT(
 map(NF_Z,{
  t^3 * f_E - (t^9 - c[2]^5*t^7 + c[2]^3 * t^3),
  F_K(aa,-bb) - dd,
  F_K(bb,-aa) - (-dd),
  F_K(aa,-cc) - (-dd+aa^9),
  F_K(cc,-aa) - ( dd-aa^9),
  F_K(cc,-bb) - (-dd-bb^9),
  F_K(bb,-cc) - ( dd+bb^9),
  f_E - (t^2-dd^2)*(t^2-(dd-aa^9)^2)*(t^2-(dd+bb^9)^2),
  (dd-aa^9)^2*(dd+bb^9)^2 - dd^4,
  coeff(f_E,t,4) - (-dd^2-(dd-aa^9)^2-(dd+bb^9)^2),
  coeff(f_E,t,4) - (-c[2]^4*(aa-bb)*dd+c[2]^5),
  c[2]^4*dd-c[2]^4*(aa-bb)*(1-c[2]*c[3]^2),
  (aa-bb)^2+c[2]+(aa+bb)*c[3]^3,
  coeff(f_E,t,4) + c[2]^5,
  coeff(f_E,t,2) - (dd^2*(dd-aa^9)^2+dd^2*(dd+bb^9)^2+(dd-aa^9)^2*(dd+bb^9)^2),
  coeff(f_E,t,2),
  coeff(f_E,t,0) + dd^6
 }) = {0},
 "claim-K"
);

# Formulae in the Claim labelled claim-L and its proof
_ASSERT(
 map(NF_Z,{
  z^9,
  z^3+(c[2]+c[2]^2*c[3]^2)*z,
  (c[3]^3-c[3]*c[2]^3)*z,
  az - z - aa*(1 - z^3*aa^5 - z^6*aa^2),
  (aa-z)*(bb-z)*(cc-z) -
   c[3]*(1-z^3*aa^5-z^6*aa^2)*
        (1-z^3*bb^5-z^6*bb^2)*
        (1-z^3*cc^5-z^6*cc^2),
  c[2]^3*z^3,
  c[3]-c[2]*z+c[1]*z^2-z^3-(c[3]-c[3]*s[5]*z^3-c[3]*s[2]*z^6-c[3]*mm*z^6),
  s[5]-c[3]*c[2],
  c[2]*z+z^3*(1-c[3]^2*c[2]),
  c[3]^4*c[2]^2*z^3,
  c[3]^2*c[2]^5*z^3,
  (c[2]+c[3]^2*c[2]^2)*z+z^3,
  az^2 - (aa^2-z*aa+z^2+z^3*aa^7+z^4*aa^6+z^6*(aa^12+aa^4)+z^7*aa^3),
  az^2 + bz^2 + cz^2 - (s[2]-z*s[1]+z^3*s[7]+z^6*(s[12]+s[4])+z^7*s[3]),
  -z*s[1]+z^3*s[7]+z^4*s[6]+z^6*(s[12]+s[4])+z^7*s[3],
  -z*s[1]+z^3*s[7]+z^7*s[3],
  s[1]+c[3]^3,
  s[3],
  s[7]-(-c[3]^5+c[3]*c[2]^2),
  z^3+(c[3]^2*c[2]^2+c[2])*z,
  (c[3]^3-c[3]*c[2]^3)*z
 }) = {0},
 "claim-L"
):

vars_ZW := plex(aa,bb,cc,c[1],z,w,c[3],c[2]):
rels_ZW := {
 op(basis_Z),
 op(subs(z = w,basis_Z))
}:
basis_ZW := Basis(rels_ZW,vars_ZW,characteristic=p):
NF_ZW := (u) -> mods(NormalForm(u,basis_ZW,vars_ZW,characteristic=p),p):

vars_ZW1 := plex(c[1],z,w,c[3],c[2]):
basis_ZW1 := remove(has,basis_ZW,{aa,bb,cc});
cobasis_ZW1 := cobasis(basis_ZW1,vars_ZW1);

points_ZW := map(NF_ZW,[
 0,z,-z,w,-w,F_K(z,w),F_K(z,-w),F_K(-z,w),F_K(-z,-w)
]):

f_ZW := t^9 + uu^3*z^2*w^2*t^7 - (uu*z^2*w^2+uu^2*(z^2+w^2))*t^3;

# Formulae in the Claim labelled claim-M and its proof

### UNFINISHED
_ASSERT(
 map(NF_ZW,{
  f_ZW - mul(t-m,m in points_ZW),
  z^3 + uu*z,
  z^9,
  uu^4*z,
  F_K(z, w) - (z+w)*(1+uu^3*z*w),
  F_K(z,-w) - (z-w)*(1-uu^3*z*w),
  F_K(z, w)^2 - (z+w)^2 - uu^3*z^2*w^2,
  F_K(z,-w)^2 - (z-w)^2 - uu^3*z^2*w^2,
  uu      * z^2*w^2*uu^3,
  z       * z^2*w^2*uu^3,
  w       * z^2*w^2*uu^3
 }) = {0},
 "claim-M"
):

vars_C := vars_ZW:
rels_C := {op(basis_ZW),coeffs(mods(expand(f_ZW - t^3*f_E),p),t)}:
basis_C := Basis(rels_C,vars_C,characteristic=p);
NF_C := (u) -> mods(NormalForm(u,basis_C,vars_C,characteristic=p),p);

vars_C1 := plex(c[1],z,w,c[3],c[2]);
basis_C1 := mods(remove(has,basis_C,{aa,bb,cc}),p);
cobasis_C1 := cobasis(basis_C1,vars_C1);

_ASSERT(nops(cobasis_C1) = 108,"C(K,G) has rank 108");

rels_C1_alt := [
c[3]^6-c[2]^5,
c[3]^3*c[2]-c[3]*c[2]^4,
c[2]^6,
c[3]^3*z-c[3]*c[2]^3*z,
c[2]^4*z,
z^3+c[2]*z+c[3]^2*c[2]^2*z,
c[3]^3*w-c[3]*c[2]^3*w,
c[2]^4*w,
w^3+c[2]*w+c[3]^2*c[2]^2*w,
uu*z^2*w^2+uu^3+uu^2*(z^2+w^2)
]:

map(NF_C,rels_C1_alt);

cza[1] := mods(expand(rem(expand(a + F_E(a,z) + F_E(a,-z)),z^8+3,z)),p^6):
cza[2] := mods(expand(rem(expand(a * F_E(a,z) + a * F_E(a,-z) + F_E(a,z)*F_E(a,-z)),z^8+3,z)),p^6):
cza[3] := mods(expand(rem(expand(a * F_E(a,z) * F_E(a,-z)),z^8+3,z)),p^6):

# Series for [3](x) in terms of y = x(x+_Fz)(x-_Fz) assuming <3>(z)=0
nn := (z,y) ->
   z^6*y+
28    *y^3-
27*z^2*y^5-
39*z^4*y^7+
 9*z^2*y^13-
27    *y^19+
27*z^4*y^23-
 9*z^6*y^25+
27*z^4*y^31-
27*z^6*y^49-
27*z^6*y^73;

# Series for c_2({x,x+_Fz,x-_Fz}) in terms of y = x(x+_Fz)(x-_Fz) assuming <3>(z)=0
mm := (z,y) ->
-z^2+
38*z^4*y^2+
35*z^6*y^4+
27    *y^6-
21*z^2*y^8-
 6*z^4*y^10+
 9*z^6*y^12+
18*z^2*y^16+
18*z^6*y^20+
 9    *y^22-
24*z^2*y^24-
15*z^4*y^26+
39*z^6*y^28+
36*z^2*y^32-
 9*z^4*y^34-
27*z^6*y^36+
27*z^2*y^40-
27*z^4*y^42+
27*z^6*y^44+
27    *y^46+
18*z^2*y^48-
36*z^4*y^50+
18*z^6*y^52-
27*z^2*y^56+
27*z^4*y^58;

_ASSERT(
 mods(expand(rem(T["p_series"](a)-nn(z,cza[3]),z^8+3,z)),p^4) = 0,
 "Defining property of nn"
);

# The assertion below actually works mod (p^4,p^3z^6) but for simplicity
# we just check mod p^3.
_ASSERT(
 mods(expand(rem(cza[2]-mm(z,cza[3]),z^8+3,z)),p^3) = 0,
 "Defining property of nn"
);

mu := 1+3^8+3^16+3^18;

# 3-adic square root of -2
rr := 5-3^3+3^5-3^6-3^7+3^11+3^12-3^14-3^15-3^17-3^18+3^20-3^22-3^23-3^24;
nu := mods((rr-2)/(3*mu),p^25);
lm := 1 + 3 * mu * nu;

_ASSERT(
 modp(rem(T["p_series"](x),x^9+3*mu*x,x),p^25) = 0,
 "Defining property of mu"
);

_ASSERT(
 modp(rr^2+2,p^25) = 0,
 "Defining property of rr = sqrt(-2)"
);

_ASSERT(
 mods(rem(rem(nn(z,y),y^3+y*z^6,y),z^8+3,z),p^4) = 0,
 "Consistency check for nn"
);

_ASSERT(
 mods(rem(rem(mm(z,y)+z^2+nu*z^4*y^2,y^3+y*z^6,y),z^8+3,z),p^4) = 0,
 "Consistency check for mm"
);

char0 := proc(t)
 local u,v;
 u := subs(LL = 1+3*MU*NU,t);
 v :=
 [rem(subs({z=0,w=0,c[2]=3*a^2,c[3]=a^3},u),a^9+3*MU*a,a),
  rem(rem(subs({w= 0,c[2]=-NU*z^4*y^2-z^2,c[3]=y},u),y^3+y*z^6,y),z^8+3*MU,z),
  rem(rem(subs({w= z,c[2]=-NU*z^4*y^2-z^2,c[3]=y},u),y^3+y*z^6,y),z^8+3*MU,z),
  rem(rem(subs({w=-z,c[2]=-NU*z^4*y^2-z^2,c[3]=y},u),y^3+y*z^6,y),z^8+3*MU,z),
  rem(rem(subs({z= 0,c[2]=-NU*w^4*y^2-w^2,c[3]=y},u),y^3+y*w^6,y),w^8+3*MU,w)
 ];
 v := simplify(subs(NU = (sqrt(-2)-2)/(3*MU),v));
 return v;
end:

char := proc(t)
 local u;
 u := char0(t);
 [seq(coeff(u[1],a,i),i=0..8),
  seq(seq(seq(coeff(coeff(u[i],y,k),z,j),k=0..2),j=0..7),i=2..4),
      seq(seq(coeff(coeff(u[5],y,k),w,j),k=0..2),j=0..7)];
end:

SS :=
 (3*MU^2*NU^2+2*MU*NU-3)/MU/((1+3*MU*NU)^3-27) * c[3]^2*(c[2]^3-27*c[3]^2) +
 (3/MU) * c[3]^4 +
 c[2]^2 +
 w^2*z^2 +
 (w^2+z^2)*c[2] +
 (NU/(1+3*MU*NU)^2) * (z^2+w^2)*c[2]^2*c[3]^2;

rels_E_lift := [
 (LL+6)*c[2]^6+(-36*LL-192)*c[3]^2*c[2]^3+(3*LL*MU+18*MU)*c[2]^2+(252*LL+864)*c[3]^4,
 (-LL*MU-3*MU)*c[2]^5+(16*LL*MU+68*MU)*c[3]^2*c[2]^2+(-3*LL*MU^2-9*MU^2)*c[2]+(-32*LL-36)*c[3]^6,
 -c[2]^4*c[3]+(4*LL+24)*c[3]^3*c[2]+(36*LL*MU-27*MU)*c[3],
 z*c[3]*(-c[2]^3+(LL+6)*c[3]^2),
 w*c[3]*(-c[2]^3+(LL+6)*c[3]^2),
 z*((c[2]^4+3*MU)-(4*LL+20)/3*c[2]*c[3]^2),
 w*((c[2]^4+3*MU)-(4*LL+20)/3*c[2]*c[3]^2),
 (LL^2*(c[2]+z^2) + (LL-1)/(3*MU)*c[2]^2*c[3]^2)*z,
 (LL^2*(c[2]+w^2) + (LL-1)/(3*MU)*c[2]^2*c[3]^2)*w,
 SS
];

ss := mods(subs({MU = mu,NU = nu},SS),p^25);

_ASSERT(
 {rem(subs({z=0,w=0,c[2]=3*a^2,c[3]=a^3},SS),a^8+3*MU,a),
  rem(rem(subs({w= 0,c[2]=-NU*z^4*y^2-z^2,c[3]=y},SS),y^3+y*z^6,y),z^8+3*MU,z),
  rem(rem(subs({w= z,c[2]=-NU*z^4*y^2-z^2,c[3]=y},SS),y^3+y*z^6,y),z^8+3*MU,z),
  rem(rem(subs({w=-z,c[2]=-NU*z^4*y^2-z^2,c[3]=y},SS),y^3+y*z^6,y),z^8+3*MU,z),
  rem(rem(subs({z= 0,c[2]=-NU*w^4*y^2-w^2,c[3]=y},SS),y^3+y*w^6,y),w^8+3*MU,w)
 } = {0},
 "Character theory relation"
);

vars_E := vars_C1:
rels_E := {op(basis_C1),mods(ss,p)}:
basis_E := Basis(rels_E,vars_E,characteristic=p);
NF_E := (u) -> mods(NormalForm(u,basis_E,vars_E,characteristic=p),p);
cobasis_E := cobasis(basis_E,vars_E);

_ASSERT(nops(cobasis_E) = 105,"E^0(BG) has rank 105");


text:
make_complex_projective_space := proc(n::nonnegint)
 local T,x,i;

 T := table([]):
 T["name"] := sprintf("CP^%d",n);
 T["dimension"] := 2*n;
 T["cohomology_gens"] := tdeg(x);
 T["cohomology_degrees"] := table([x = 2]);
 T["cohomology_rels"] := [x^(n+1)];
 T["cohomology_basis"] := [seq(x^i,i=0..n)];
 T["cohomology_graded_basis"] := table([]):
 for i from 0 to 2*n do
  T["cohomology_graded_basis"][i] :=
   `if`(modp(i,2)=0,[x^(i/2)],[]);
 od:

 return eval(T);
end:

ext: mpl


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

`is_element/P2` := (N::posint) -> (A::set) -> proc(u)
 local n,A2,ab,ba,v,i;
 global reason;

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

 n := nops(A);

 A2 := `list_elements/pairs`(A);
 if {indices(u)} <> {op(A2)} then
  reason := [convert(procname,string),"u is not indexed by pairs in A",u,A2];
  return false;
 fi;

 for ab in A2 do
  v := u[op(ab)];
  if not(`is_element/R`(N)(v)) then
   reason := [convert(procname,string),"u[a,b] is not in R^N",op(ab),v,N];
   return false;
  fi;
  if simplify(add(v[i]^2,i=1..N) - 1) <> 0 then
   reason := [convert(procname,string),"u[a,b] is not in S^(N-1)",op(ab),v,N];
   return false;
  fi;
 od;

 for ab in A2 do
  ba := [ab[2],ab[1]];
  if u[op(ab)] +~ u[op(ba)] <> [0$N] then
   reason := [convert(procname,string),"u[a,b]+u[b,a] <> 0",op(ab),u[op(ab)],u[op(ba)],N];
   return false;
  fi;
 od;

 return true;
end:

`phi2/Fbar/P2` := (N::posint) -> (A::set) -> proc(x)
 local A2,u,ab,a,b,f;

 A2 := `list_elements/pairs`(A);
 u := table();

 for ab in A2 do
  a,b := op(ab);
  f := `normalise_2/F`(N)({a,b})(x[{op(ab)}]);

  u[op(ab)] := combine(sqrt(2) *~ f[b]);
 od:

 return(eval(u));
end:


Default

songs:
Baby I Love Your Way.mp3
Lost in your eyes.mp3

last: Lost in your eyes.mp3


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

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

f0 := s -> s^2/(1 + params_a_a*s^2):
f := x -> 1
  + params_a_c1*f0(X2S*x)
  + params_a_c2*f0(X2S*x)^2
  + params_a_c3*f0(X2S*x)^3:

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
program TeaBreak;
begin
;;;;;;;;;;
writeln
(
'Say "BREAK!!" '
)
;
;;;;;;;;;;
break
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing


text:
# A full flag of subspaces in (Z/p)^d is represented by a list
# [w[1],...,w[d]] where the subspace of dimension i is
# spanned by [w[1],...,w[i]] and the last nonzero entry in
# w[i] is a 1 in position h[i] say, and w[i+1],...,w[d]
# all have zero in position h[i].

`count_elements/vectors` := (p,d) -> p^d;

`list_elements/vectors` := proc(p,d)
 local L,i,j,u;

 L := [[]];
 for i from 1 to d do
  L := [seq(seq([op(u),j],j=0..p-1),u in L)];
 od:
 return L;
end:

`height/vectors` := (p,d) -> proc(u)
 local h;

 h := d;
 while h > 0 and u[h] = 0 do h := h - 1; od;
 return h;
end:

`count_elements/monic_vectors` := (p,d) -> (p^d - 1)/(p - 1);

`list_elements/monic_vectors` := proc(p,d)
 local i,u;

 [seq(seq([op(u),1,0$(d-1-i)],u in `list_elements/vectors`(p,i)),i=0..d-1)];
end:

`count_elements/full_flags` := proc(p,d) local i; mul((p^i-1)/(p-1),i=1..d); end:

`list_elements/full_flags` := proc(p,d)
 local L,L1,nL,H,V,i,j,k,u,v,q;
 L := [[]];
 H := [[]];

 for i from 1 to d do
  nL := nops(L);
  L1 := L;
  L := [];
  V := `list_elements/monic_vectors`(p,d - i + 1);
  H := map(W -> sort(map(`height/vectors`(p,d),W)),L1);
  for j from 1 to nL do
   for v in V do
    u := [];
    q := 1;
    for k from 1 to d do
     if member(k,H[j]) then
      u := [op(u),0];
     else
      u := [op(u),v[q]];
      q := q + 1;
     fi;
    od:
    L := [op(L),[op(L1[j]),u]];
   od;
  od:
 od:

 return L;
end:

`random_element/full_flags` := (p,d) -> proc()
 local V,i,v;
 V := table():
 V[0] := `bot/subspaces`(p,d);
 for i from 1 to d do
  v := [0$d];
  while `is_member/subspaces`(p,d)(v,V[i-1]) do
   v := [seq(rand(0..p-1)(),i=1..d)];
  od;
  V[i] := `span/subspaces`(p,d)([op(V[i-1][1]),v]);
 end:
 return `from_subspace_list/full_flags`(p,d)([seq(V[i],i=1..d)]);
end:

`is_element/full_flags` := (p,d) -> proc(W)
 local h,i,j;

 if not(type(W,list(list(integer)))) then
  return false;
 fi;

 if map(nops,W) <> [d$d] then
  return false;
 fi;

 if W <> modp(W,p) then
  return false;
 fi;

 h := table():
 for i from 1 to d do
  if W[i] = [0$d] then return false; fi;
  j := d;
  while W[i][j] = 0 do j := j - 1; od;
  if W[i][j] <> 1 then
   return false;
  fi;
  h[i] := j;
 od:

 for i from 1 to d do
  for j from i+1 to d do
   if W[j][h[i]] <> 0 then return false; fi;
  od;
 od;

 return true;
end:

`count_elements/full_flags` := proc(p,d) local i; mul((p^i-1)/(p-1),i = 1..d); end:

`basepoint/full_flags` := (p,d) -> convert(IdentityMatrix(d),listlist);

`to_subspace_list/full_flags` := (p,d) -> proc(W)
 local i;
 [seq(`span/subspaces`(p,d)([op(1..i,W)]),i=1..d)];
end:

`from_subspace_list/full_flags` := (p,d) -> proc(V)
 local i,j,h,v,W;

 v := V[1][1][1];
 j := d; while j > 0 and v[j] = 0 do j := j - 1; od;
 v := modp(v /~ v[j],p);
 W := [v];
 h := [j];
 for i from 2 to d do
  j := 1;
  while j < i and `is_member/subspaces`(p,d)(V[i][1][j],V[i-1]) do
   j := j + 1;
  od;
  v := V[i][1][j];
  for j from 1 to i-1 do
   v := modp(v -~ v[h[j]] *~ W[j],p);
  od:
  j := d; while j > 0 and v[j] = 0 do j := j - 1; od;
  v := modp(v /~ v[j],p);
  W := [op(W),v];
  h := [op(h),j];
 od:
 return W;
end:

# It would be better to do this directly and use it in the
# `from_subspace/full_flags` procedure.
`from_generators/full_flags` := (p,d) -> proc(U)
 local V,i;
 V := [seq(`span/subspaces`(p,d)([op(1..i,U)]),i=1..d)];
 return `from_subspace_list/full_flags`(p,d)(V);
end:

`matrix_act/full_flags` := (p,d) -> proc(g,W)
 local gW,V,i;

 gW := convert(Transpose(Matrix(g) . Transpose(Matrix(W))),listlist);
 V := [seq(`span/subspaces`(p,d)([op(1..i,gW)]),i=1..d)];
 return `from_subspace_list/full_flags`(p,d)(V);
end:

`jordan_permutation/full_flags` := (p,d) -> proc(U,V)
 local U0,V0,UV,UV1,i,j,w;
 U0 := `to_subspace_list/full_flags`(p,d)(U);
 U0 := table([0 = `bot/subspaces`(p,d),seq(i = U0[i],i=1..d)]);
 V0 := `to_subspace_list/full_flags`(p,d)(V);
 V0 := table([0 = `bot/subspaces`(p,d),seq(i = V0[i],i=1..d)]);
 UV := table():
 for i from 0 to d do
  for j from 0 to d do
   UV[i,j] := `inter/subspaces`(p,d)(U0[i],V0[j]);
  od;
 od;
 UV1 := table();
 w := table():
 for i from 1 to d do
  for j from 1 to d do
   UV1[i,j] := `sum/subspaces`(p,d)(UV[i-1,j],UV[i,j-1]);
   if `dim/subspaces`(p,d)(UV[i,j]) > `dim/subspaces`(p,d)(UV1[i,j]) then
    w[i] := j;
   fi;
  od;
 od;
 return [seq(w[i],i=1..d)];
end:


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

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

RigidMotionsParameterSpaceDecompostion := module()
  option package;
  uses   RigidMotionsParameterSpaceDecompostionRecursive, RigidMotionsParameterSpaceCommon;

  local  GetQuadric, IsMonotonic, ComputeSetOfQuadrics,
         ComputeEventsATypeGrid, ComputeEventsBTypeGrid, ComputeEventsCTypeGrid,
         ComputeAsymptoticABEventsGrid, ComputeAsymptoticAAEvents,
         ComputeEventsFromAlgebraicNumbers, IsAsymptotic;

         # Only the procedures with prefix Launch* should be called! Other are exported only for
         # the grid framework.
  export IsAsymptoticIntersection, ComputeSamplePoints, ParallelComputeSamplePoints,
         LaunchComputeEvents, LaunchComputeSamplePoints;

  #Variables shared by grid nodes;
  global Q, events, vars, dbPath, skipped;

# Procedure: GetQuadric
#   Compute a quadric.
#
# Parameters:
#   R          - a 3x3 rotation matrix obtained from Cayley transform
#   neighbor   - a vector which points on a neighbor
#   hGridPlane - a half-grid plane i.e. k_dim + 1/2
#   axis       - a given row of R matrix
#
# Output:
#   Multivariate polynomial of degree 2 related to the changed of configuration of image patch and
#   given in [Quadrics:2016] as equation (7).
GetQuadric := proc(R::~Matrix, neighbor::~Vector, hGridPlane::~Vector, axis::integer)
  local r := R . neighbor - hGridPlane;
  if axis > 0 and axis < 4 then
    return simplify( ( 2 * denom( R[1][1] ) ) * r[axis] );
  else
    error "Wrong dimension: dim must be in [1,3].";
  end if:
end proc:

# Procedure: IsMonotonic
#   Check if given polynomial of degree 2 is non-positive or non-negative
#
# Parameters:
#   x      - a polynomial
#
# Output:
#   true if polynomial of degree 2 is non-positive and false otherwise
IsMonotonic := proc( x::~polynom )
  local homo, hessian, signmap, clean:
  homo := Groebner:-Homogenize( x, v ):
  hessian := 1 / 2 * VectorCalculus:-Hessian( homo, [ op( indets( x ) ), v ] ):
  signmap := map( signum, LinearAlgebra:-Eigenvalues( hessian ) ):
  signmap := convert( signmap, list ):
  clean := remove( `=`, signmap , 0 ):
  return ( andmap( `=`, clean, 1 ) or andmap( `=`, clean, -1 ) ):
end proc:

# Procedure: ComputeSetOfQuadrics
#   Compute a set of quadrics which are reduced by duplicated and these ones
#   which are strictly positive or negative.
#
# Parameters:
#   R          - a rotation matrix obtained from Cayley transform
#   nType      - size of neighborhood i.e. N_1, N_2, N_3.
#   axis       - a given axis i.e. 1 = x, 2 = y and 3 = z
#   kRange     - a range of planes to consider
#
# Output:
#   List of quadrics of form f - g reduced by duplicated---up to constant
#   factor---and these ones which are strictly positive or negative.
ComputeSetOfQuadrics := proc( R::~Matrix,
                               nType::string,
                               axis::~integer,
                               kRange::~list )
  local neighborhood := GetNeighborhood( nType ):
  local quadrics := Array([]):
  local f::polynom, g::polynom, quadric::polynom:
  local indexes := []:
  local T := combinat:-cartprod( [ neighborhood, neighborhood, kRange, kRange ] ):

  if LinearAlgebra:-Determinant(R) <> 1 or [upperbound(R)] <> [3,3] or nops(indets(R)) = 0 then
    error "Used matrix is not a correct 3 x 3 rotation matrix!";
  fi;

  if axis < 1 or axis > 3 then
    error "Axis is not from range [1,3].";
  fi:

  while not T[ finished ] do
    indexes := T[ nextvalue ]():
    f := GetQuadric( R, Vector( indexes[1] ), Vector( 3, indexes[3] + 1/2 ), axis ):
    g := GetQuadric( R, Vector( indexes[2] ), Vector( 3, indexes[4] + 1/2 ), axis ):
    quadric := f - g:
    if not type(quadric, constant) then
      quadric := quadric / lcoeff( quadric ):
      ArrayTools:-Append(quadrics,quadric):
    end if:
  end do:
  quadrics := convert(quadrics,set):
  quadrics := remove( IsMonotonic, quadrics ):
  quadrics := Threads:-Map(proc(x) normal(x) * denom(x) end proc, quadrics):
  return quadrics:
end proc:

# Procedure: IsAsymptotic
#   Checks if given quadric has an asymptotic critical value. For this moment a direction is fixed
#   to a.
#
# Parameters:
#   x          - a quadric in three variables
#   vars       - a list of variables
#
# Output:
#   List of solution for which partial derivatives in vars[2] and vars[3] are collinear.
IsAsymptotic := proc(x::polynom, vars::list)
  local vec, Vb, Vc, VV, sols;
  if nops(vars) < 3 then
    error "Expected at least three variables!";
  fi;
  vec := VectorCalculus:-Gradient(x, vars);
  Vb := Vector(3, [coeff(vec[2],vars[2]),coeff(vec[2],vars[3]),eval(vec[2],[vars[-2]=0,vars[-1]=0])]);
  Vc := Vector(3, [coeff(vec[3],vars[2]),coeff(vec[3],vars[3]),eval(vec[3],[vars[-2]=0,vars[-1]=0])]);
  VV := LinearAlgebra:-CrossProduct(Vb, Vc);
  if norm(VV,1) = 0 then
    return {{vars[1]=0}};
  fi;
  sols := solve({VV[1] = 0, VV[2] = 0,VV[3] = 0,vars[1]>=0}, [vars[1]]);
  if sols = NULL then
    return {};
  else
    return sols;
  fi;
end proc:

# Procedure: IsAsymptoticIntersection
#   Checks if intersection of two quadrics has an asymptotic critical value. For this moment a
#   direction is fixed to a.
#
# Parameters:
#   p          - a quadric in three variables
#   q          - a quadric in three variables
#   vars       - a list of variables
#
# Output:
#   List of solution for which intersection of p and q has an asymptotic intersection
# Comment:
#   - Since Groebner package seems to have memory leak I should rather replace
#     PolynomialIdeals:-EliminationIdeal by resultant elimination similarly to what I did with
#     univariate polynomials.
# TODO:
#   - Allow user to chose a direction.
#   - if there is no intersection between quadrics then skip it.
IsAsymptoticIntersection := proc( p::polynom, q::polynom, vars::list )
  local J := PolynomialIdeals:-`<,>`(p,q);
  local Pb, Pc, Cb, Cc, sols, univ;
  if nops(vars) < 3 then
    error "Expected at least three variables!";
  fi;
  Pb := PolynomialIdeals:-EliminationIdeal(J,{op(vars[1..2])}):
  Pb := PolynomialIdeals:-IdealInfo:-Generators(Pb)[1];
  Pc := PolynomialIdeals:-EliminationIdeal(J,{op(vars[[1,3]])});
  Pc := PolynomialIdeals:-IdealInfo:-Generators(Pc)[1];
  Cb := lcoeff(Pb, vars[2]);
  Cc := lcoeff(Pc, vars[3]);
  return lcm(Cb,Cc);
end proc:

# Procedure: ComputeEventsATypeGrid
#   Compute events such that a sweep plane is tangent to a quadric.
#
# Parameters:
#   Q          - a set of quadrics or conics
#   dim        - a list of indexes of variables used to calculate partial derivatives
#   vars       - a list of variables
#
# Output:
#   List of ranges which contains roots of a system(q, d/db q, d/dc q).
#
# Comment:
#  - only the first direction is supported
ComputeEventsATypeGrid := proc( Q, dim::list, vars::list )
  local s, result := Array([]);
  if nops(vars) < 3 then
    error "Expected at least three variables!";
  fi;
  s := proc(i::integer, vars::list)
   local sys, univ;
   local q := Q[i];
   sys := [q, diff( q, vars[ dim[1] ] ), diff( q, vars[ dim[2] ] )];
   univ := UnivariatePolynomial(sys, vars);
   return SerializeEvents(GenerateEvents(univ, [i]));
  end proc:
  map[inplace](proc(x) ArrayTools:-Extend(result, x, inplace=true) end proc,
                                             [Grid:-Seq(s(i, vars), i=1..nops(Q))]);
  return ReconstructEvents(result);
end proc:

# Procedure: ComputeEventsBType
#   Compute events such that intersection of two quadrics is tangent to sweeping plane.
#
# Parameters:
#   dir        - a direction of a gradient product it should
#                be the same as director of sweep
#   Q          - a set of conics
#   vars       - a list of variables
#
# Output:
#   Indexes of quadrics which intersect and a component of a vector product of
#   their gradients in given direction have a common root.
ComputeEventsBTypeGrid := proc( Q, dir::integer, vars::list )
  local s, result := Array([]);
  if nops(vars) < 3 then
    error "Expected at least three a variables!";
  fi;
  s := proc(i, j, vars::list)
    local p, prod, univ, sys:
    prod := LinearAlgebra:-CrossProduct( VectorCalculus:-Gradient( Q[i], vars ),
                                    VectorCalculus:-Gradient( Q[j], vars ) )[dir]:
    sys := [Q[i], Q[j], prod];
    univ := UnivariatePolynomial(sys, vars);
    return SerializeEvents(GenerateEvents(univ, [i, j]));
  end proc:
  map[inplace](proc(x) ArrayTools:-Extend(result, x, inplace=true) end proc,
                                   [Grid:-Seq(seq(s(i, j, vars), j=i+1..nops(Q)), i=1..nops(Q))]);
  return ReconstructEvents(result);
end proc;

# Procedure: ComputeEventsCTypeGrid
#   Compute events such that three quadrics intersects in a point.
#
# Parameters:
#   Q          - a set of quadrics
#   vars       - a list of variables
#
# Output:
#   Indexes of quadrics which intersect in a point.
ComputeEventsCTypeGrid := proc( Q, vars::list )
  uses ArrayTools;
  local result := Array([]), s;
  if nops(vars) < 3 then
    error "Expected at least three a variables!";
  fi;
  s := proc (i, j, k, vars::list)
    local univ, sys;
    sys := [Q[i], Q[j], Q[k]];
    univ := UnivariatePolynomial(sys, vars);
    return SerializeEvents(GenerateEvents(univ, [i, j, k]));
  end proc;
  map[inplace](proc(x) Extend(result, x, inplace=true) end proc, [Grid:-Seq(seq(seq(s(i, j, k, vars),
                                       k=j+1..nops(Q)),j=i+1..nops(Q)),i=1..nops(Q))]);
  return ReconstructEvents(result);
end proc:

# Procedure: ComputeAsymptoticAAEvents
#   Compute real algebraic numbers which corresponds to
#   asymptotic cases given by one quadrics.
#
# Parameters:
#   Q          - a set of quadrics
#   vars       - a list of variables
#
# Output:
#   A list of real algebraic numbers and indexes of quadrics -- events.
ComputeAsymptoticAAEvents:=proc(Q, vars::list)
  uses ArrayTools;
  local result := Array([]), s;
  if nops(vars) < 3 then
    error "Expected at least three a variables!";
  fi;
  s:=proc(i::integer, vars::list)
    local events := Array([]), sol, rootsF, tmp:
    rootsF := IsAsymptotic(Q[i], vars):
    for sol in rootsF do
      sol := op(sol);
      if not type(rhs(sol), rational) then
        error "Irrational asymptotic case! Are you sure the input is a set of quadrics?"
      fi:
      ArrayTools:-Append(events, EventType(RealAlgebraicNumber(lhs(sol) * denom(rhs(sol)) -
      numer(rhs(sol)), rhs(sol), rhs(sol)), [i]));
    od:
    return events;
  end proc;
  map[inplace](proc(x) Extend(result, x, inplace=true) end proc, [seq(s(i, vars), i=1..nops(Q))]);
  return result;
end proc;

# Procedure: ComputeAsymptoticABEvents
#   Compute real algebraic numbers which corresponds to
#   asymptotic cases given by one quadrics.
#
# Parameters:
#   Q          - a set of quadrics
#   vars       - a list of variables
#
# Output:
#   A list of real algebraic numbers and indexes of quadrics -- events.
ComputeAsymptoticABEventsGrid:=proc(Q, vars::list)
  uses ArrayTools;
  local result := Array([]), s;
  if nops(vars) < 3 then
    error "Expected at least three a variables!";
  fi;
  s:=proc(i::integer, j::integer, vars::list)
   local poly;
   poly := RigidMotionsParameterSpaceDecompostion:-IsAsymptoticIntersection(Q[i], Q[j], vars);
   if poly = NULL or nops(poly) = 0 then
     return [];
   fi;
   return SerializeEvents(GenerateEvents(poly, [i, j]));
  end proc;
  map[inplace](proc(x) Extend(result, x, inplace=true) end proc, [Grid:-Seq(seq(s(i, j, vars),
  j=i+1..nops(Q)), i=1..nops(Q))]);
  return ReconstructEvents(result);
end proc;

# Procedure: ComputeEventsFromAlgebraicNumbers
#   Compute and sort events
#
# Parameters:
#   Q     - set of quadrics or conics
#   vars       - a list of variables
# Output:
#   Sorted set of events
ComputeEventsFromAlgebraicNumbers := proc( Q, vars::list )
  local events := Array([]);
  ArrayTools:-Extend(events, ComputeEventsATypeGrid( Q, [2, 3], vars ), inplace=true);
  ArrayTools:-Extend(events, ComputeEventsBTypeGrid( Q, 1, vars ), inplace=true);
  ArrayTools:-Extend(events, ComputeEventsCTypeGrid( Q, vars ), inplace=true);
  ArrayTools:-Extend(events, ComputeAsymptoticAAEvents(Q, vars), inplace=true);
  ArrayTools:-Extend(events, ComputeAsymptoticABEventsGrid(Q, vars), inplace=true);
  return AlgebraicSort(events);
end proc:

# Procedure: LaunchComputeEvents
#   Computes sorted list of events for rotational part of rigid motions using the grid framework
#
#
# Parameters:
#   variables     - list of variables in which the problem is expressed
#   databasePath  - a path to a copy of the database CompRegister.db
#   nType         - neighborhood type: N1, N2 or N3.
#   kRange        - range of grid lines passed as a list
# Output:
#   It populates a database given by databasePath.
LaunchComputeEvents := proc(variables::list, databasePath::string, nType::string, kRange::list)
  local Q, Q1, Q2, Q3, R, i, db := Object(ComputationRegister, databasePath);
  vars := variables;
  dbPath:=databasePath;
  R := CayleyTransform(variables);
  Q1 := ComputeSetOfQuadrics(R, nType, 1, kRange);
  Q2 := map2(subs, [a = b, b = c, c = a], Q1);
  Q3 := map2(subs, [a = b, b = c, c = a], Q2);
  Q := [op(Q1), op(Q2), op(Q3), op(variables)];
  for i from 1 to nops(Q) do
    InsertQuadric(db, i, Q[i]);
  od;
  SynchronizeQuadrics(db);
  events := ComputeEventsFromAlgebraicNumbers(Q, variables);
  events := select[flatten](proc(x) evalb(GetInterval(GetRealAlgebraicNumber(x))[2] >= 0) end proc,
                                                                                events);
  events := ReduceEvents(events);
  AdjustEvents(events, upperbound(Q), variables);
  #Insert events into the register
  for i from 1 to upperbound(events) do
    InsertEvent(db, i, events[i]);
  od;
  SynchronizeEvents(db);
  Close(db);
end proc:

# Procedure: ComputeSamplePoints
#   Computes sample points for rotational part of rigid motions
#
#
# Parameters:
#   Q                  - list of quadrics
#   events             - an array of events
#   first              - integer value which indicates a first event to proceed.
#   last               - integer value which indicates a last event to proceed.
#   vars               - list of variables in which conics are expressed
#   db                 - an instance of the class ComputationRegister
#   skipped            - a list of the events' indices to be skipped
#
# Output:
#   It populates a database, given by databasePath, with sample points.
ComputeSamplePoints := proc (Q, events::Array, first::integer, last::integer,
                             vars::list, db::ComputationRegister, skipped::list:=[])
local i, midpoint, sys, samplePoints, disjointEvent:=[], ranumI, ranumJ;
  if first < 0 or last < 0 or last < first or upperbound(events) <= last then
    error "Bounds of the list of the events range are incorrect.":
  end if:
  for i from first to last do
    if i in skipped then
      next;
    fi:
    sys := Q[GetQuadrics(events[i])];
    ranumI := GetRealAlgebraicNumber(events[i]);
    ranumJ := GetRealAlgebraicNumber(events[i+1]);
    disjointEvent:=DisjointRanges(ranumI, ranumJ);
    midpoint := (GetInterval(disjointEvent[1])[2] + GetInterval(disjointEvent[2])[1])/2:
    # never call eval with sets!!
    sys := eval(sys, vars[1] = midpoint);
    LaunchComputeSamplePoints2D(sys, midpoint, 1, false, vars[2..], db);
    SynchronizeSamplePoints(db);
    InsertComputedNumber(db, i);
  end do;
end proc:

# Procedure: ParallelComputeSamplePoints
#   Computes sample points for rotational part of rigid motions. It should be call via Grid
#   framework.
ParallelComputeSamplePoints := proc()
  local me, n, events, skipped, first, last, offset, noEvents;
  local db:=Object(ComputationRegister, dbPath);
  me := Grid:-MyNode();
  skipped := FetchComputedNumbers(db);
  offset := FetchLowerEventID(db) - 1;
  noEvents := NumberOfEvents(db);
  # events-1 because the last event is a copy of events[-2]
  n := trunc((noEvents-1) / Grid:-NumNodes());
  first := offset + me * n + 1; last := offset + (me + 1) * n;
  if me = Grid:-NumNodes() - 1 then
    last := last + ((noEvents + offset) - (last + 1));
  fi;
  # recreate events
  events := FetchEvents(db, first, last + 1);
  RigidMotionsParameterSpaceDecompostion:-ComputeSamplePoints(Q, events, first , last, vars, db, skipped);
  Close(db);
  Grid:-Barrier();
end proc:

# Procedure: LaunchComputeSamplePoints
#   Resumes computations of sample points for rotational part of rigid motions using
#   the grid framework.
#
#
# Parameters:
#   variables     - list of variables in which the problem is expressed
#   databasePath  - a path to a database file. If file does not exist it will be crated.
#   nType         - neighborhood type: N1, N2 or N3.
#   nodes         - number of nodes used in the parallel computations
# Output:
#   It populates a database given by databasePath.
LaunchComputeSamplePoints := proc(variables::list, databasePath::string, nType::string,
                                 nodes:=kernelopts(numcpus))
  local db:=Object(ComputationRegister, databasePath);
  vars:=variables;
  dbPath:=databasePath;
  Q := FetchQuadrics(db);
  Close(db);
  Grid:-Setup("local");
  Grid:-Launch(RigidMotionsParameterSpaceDecompostion:-ParallelComputeSamplePoints,
               imports=['Q', 'vars', 'dbPath'], numnodes=nodes, allexternal=false);
end proc;

end module:

max_stars_repo_name: copyme/MapleTools

max_issues_repo_name: copyme/MapleTools

max_forks_repo_name: copyme/MapleTools


text:
program change;
var n, count : integer;
begin
    writeln('please input change');
    readln(n);
    count := 0;
    while n > 0 do begin
        while n >= 10 do begin
            while n >= 100 do begin
                while n >= 1000 do begin
                    while n >= 10000 do begin
                        count := count + 1;
                        n := n - 10000;
                    end;
                    if count > 0 then
                        writeln('10000 yen : ', count);
                    count := 0;
                    if n < 1000 then
                        break;
                    count := count + 1;
                    n := n - 1000;
                end;
                if count > 0 then
                    writeln(' 1000 yen : ', count);
                count := 0;
                if n < 100 then
                    break;
                count := count + 1;
                n := n - 100;
            end;
            if count > 0 then
                writeln('  100 yen : ', count);
            count := 0;
            if n < 10 then
                break;
            count := count + 1;
            n := n - 10;
        end;
        if count > 0 then
            writeln('   10 yen : ', count);
        count := 0;
        if n < 1 then
            break;
        count := count + 1;
        n := n - 1;
    end;
    if count > 0 then
        writeln('    1 yen : ', count);
end.

max_issues_repo_name: Hiroya-W/lang-processing


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

# make gensym global, so that it can be shared with other 'global' routines
gensym := module ()
  export ModuleApply, SymbolsToGen;
  local gs_counter, utf8, blocks, radix, unicode, bad_name, do_gensym, ModuleLoad;
  gs_counter := -1;
  utf8 := proc(n :: integer, $)
    local m;
    if n<128 then n
    elif n<2048 then 192+iquo(n,64,'m'), 128+m
    elif n<65536 then 224+iquo(n,4096,'m'), 128+iquo(m,64,'m'), 128+m
    elif n<2097152 then 240+iquo(n,262144,'m'), 128+iquo(m,4096,'m'), 128+iquo(m,64,'m'), 128+m
    elif n<67108864 then 248+iquo(n,16777216,'m'), 128+iquo(m,262144,'m'), 128+iquo(m,4096,'m'), 128+iquo(m,64,'m'), 128+m
    elif n<2147483648 then 248+iquo(n,1073741824,'m'), 128+iquo(m,16777216,'m'), 128+iquo(m,262144,'m'), 128+iquo(m,4096,'m'), 128+iquo(m,64,'m'), 128+m
    end if
  end proc;
  SymbolsToGen :=
             [ ["`", "_"]
             , ["0", "9"]
             , ["a", "z"]
             , ["A", "Z"] ];
             # e.g.
             # [[ 19968, 20950 ]] for unicode

  bad_name := n -> type(n, protected) or (x -> assigned(x))(n);
  blocks := proc($)
    map(((l,u)->`[]`(StringTools:-Ord(l), StringTools:-Ord(u) - StringTools:-Ord(l)+1))@op,SymbolsToGen)
  end proc;
  radix := proc($) `+`(op(map2(op, 2, blocks()))) / 2; end proc;
  unicode := proc(nn, $)
    local n, b;
    n := nn;
    for b in blocks() do
      if n < op(2,b) then return n + op(1,b) else n := n - op(2,b) end if
    end do
  end proc;
  ModuleApply := proc(x::name, $)
    local x1 := do_gensym(x);
    if bad_name(x1) then ModuleApply(x) else x1 end if;
  end proc;
  do_gensym := proc(x::name, $)
    gs_counter := gs_counter + 1;
    cat(x, op(map(StringTools:-Char, map(utf8 @ unicode, applyop(`+`, 1, map(`*`, convert(gs_counter, 'base', radix()), 2), 1)))))
  end proc;
  ModuleLoad := proc()
    unprotect(SymbolsToGen);
  end proc;
  ModuleLoad();
end module: # gensym

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

BindingTools := module ()
  option package;
  export generic_evalat, generic_evalatstar, load_types_from_table;

  load_types_from_table := proc(types, unload::truefalse:=true, $)
    local ty_nm;
    if unload then
      for ty_nm in [ indices(types, nolist) ] do
        if TypeTools[Exists](ty_nm) then TypeTools[RemoveType](ty_nm) end if;
      end do;
    end if;

    for ty_nm in [ indices(types, nolist) ] do
      TypeTools[AddType]( ty_nm, types[ty_nm] );
    end do;
  end proc;

  generic_evalat := proc(vv::{name,list(name)}, mm, eqs, $)
    local v, m, eqsRemain, eq, rename, funs;
    funs := map2(op, 0, indets(mm, 'function'));
    eqsRemain := select(
      (eq -> op(1,eq) <> op(2,eq)
             and not has(op(1,eq), vv)
             and (not (op(1,eq) :: 'name')
                  or depends(mm, op(1,eq))
                  or member(op(1,eq), funs))),
      eqs);
    if nops(eqsRemain) = 0 then return vv, mm end if;
    m := mm;
    rename := proc(v::name, $)
      local vRename;
      if depends(eqsRemain, v) then
        vRename := gensym(v);
        m := subs(v=vRename, m);
        vRename
      else
        v
      end if
    end proc;
    if vv :: name then
      v := rename(vv)
    else
      v := map(rename, vv);
    end if;
    v, eval(m, eqsRemain);
  end proc:

  generic_evalatstar := proc(body, bound::list, eqs, $)
    local indefinite, e, n, b, j;
    e := foldl(((x,i) -> `if`(i::`=`,
                              lam(lhs(i), rhs(i), x),
                              lam(i, indefinite, x))),
               body, op(bound));
    e := eval(e, eqs); # Piggyback on `eval/lam`
    n := nops(bound);
    b := table();
    for j from n to 1 by -1 do
      b[j] := `if`(op(2,e)=indefinite, op(1,e), `=`(op(1..2,e)));
      e := op(3,e);
    end do;
    e, [entries(b, 'nolist', 'indexorder')]
  end proc:
end module: # BindingTools

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
// Create a portal for each entity with the name ACV_remote_createPortals_PARAM
// The color is specified by tagging the entity with one of: ACV_Blue, ACV_Orange, ACV_Purple or ACV_Red
// Params: entity named ACV_remote_createPortals_PARAM
impulse process ACV_remote_createPortals {
/execute @e[name=ACV_remote_createPortals_PARAM] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_Up",Tags:[ACV_Direction],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals_PARAM] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_Down",Tags:[ACV_Direction],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals_PARAM] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_South",Tags:[ACV_Direction],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals_PARAM] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_East",Tags:[ACV_Direction],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals_PARAM] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_North",Tags:[ACV_Direction],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals_PARAM] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_West",Tags:[ACV_Direction],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals_PARAM] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_Main",Duration:2147483647}

/execute @e[name=ACV_remote_createPortals_PARAM,tag=ACV_Blue] ~ ~ ~ scoreboard players set @e[type=area_effect_cloud,name=ACV_Main,r=0] ACV_COLOR 0
/execute @e[name=ACV_remote_createPortals_PARAM,tag=ACV_Orange] ~ ~ ~ scoreboard players set @e[type=area_effect_cloud,name=ACV_Main,r=0] ACV_COLOR 1
/execute @e[name=ACV_remote_createPortals_PARAM,tag=ACV_Purple] ~ ~ ~ scoreboard players set @e[type=area_effect_cloud,name=ACV_Main,r=0] ACV_COLOR 2
/execute @e[name=ACV_remote_createPortals_PARAM,tag=ACV_Red] ~ ~ ~ scoreboard players set @e[type=area_effect_cloud,name=ACV_Main,r=0] ACV_COLOR 3
/kill @e[name=ACV_remote_createPortals_PARAM]
ACV_createPortals()
}


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

br89_a := [0, 1, 0, -2, 0, 1]:
fw := t -> mgga_series_w(br89_a, 6, t):

br89_gamma := 0.8:
br89_Q := (x, t, u) -> (u - 4*t*br89_gamma*Fermi_D(x, t))/6:
br89_y := (x, t, u) -> 2*Pi^(2/3)/(3*br89_Q(x, t, u)):

(* lower piece *)
pgk_a1 := 1.5255251812009530:
pgk_a2 := 0.4576575543602858:
pgk_a3 := 0.4292036732051034:

pgk_b  := [0.4771976183772063, -1.7799813494556270, 3.8433841862302150,
       -9.5912050880518490, 2.1730180285916720, -30.425133851603660]:

pgk_c  := [0.7566445420735584, -2.6363977871370960, 5.4745159964232880,
       -12.657308127108290, 4.1250584725121360, -30.425133957163840]:

pgk_d  := [0.00004435009886795587, 0.58128653604457910, 66.742764515940610,
       434.26780897229770, 824.7765766052239000, 1657.9652731582120]:

pgk_e  := [0.00003347285060926091, 0.47917931023971350, 62.392268338574240,
       463.14816427938120, 785.2360350104029000, 1657.962968223273000000]:

pgk_UB := 2.085749716493756:

pgk_x_lower := y -> (-arctan(pgk_a1*y + pgk_a2) + pgk_a3) *
            add(pgk_c[i]*y^(i-1), i=1..6)/add(pgk_b[i]*y^(i-1), i=1..6):

pgk_x_upper := y -> (arccsch(pgk_UB*y) + 2) *
            add(pgk_d[i]*y^(i-1), i=1..6)/add(pgk_e[i]*y^(i-1), i=1..6):

pgk_x := (x, t, u) -> convert(piecewise(
  br89_y(x, t, u) <= 0, pgk_x_lower(br89_y(x, t, u)),
  pgk_x_upper(br89_y(x, t, u))), 'Heaviside'):

br89_v_full := x -> exp(x/3.0)*(1.0 - exp(-x)*(1.0 + x/2.0))/x:
br89_v_expa := x -> 1/2 + x/6 - x^2/18:

br89_v := (x, t, u) ->
       -2*Pi^(1/3)/X_FACTOR_C * br89_v_full(pgk_x(x, t, u)):

f   := (rs, x, t, u) -> -br89_v(x, t, u)/2:


text:
#---------------------------------------------------------------------
# ieeedouble converts a number to IEEE double extended format.
# returns sign (-1 or 1), exponent between -16383 and 16383, mantissa as a fraction between 0.5 and 1.

# TODO : use JMM procedure; check subnormals etc
ieeedoubleExt:=proc(xx)
local x, sign, logabsx, exponent, mantissa, infmantissa;
    x:=evalf(xx):
    if (x=0) then
        sign,exponent,mantissa := 0,0,0;
    else
        if (x<0) then sign:=-1:
        else sign:=1:
        fi:
        exponent := floor(log2(sign*x));
        if (exponent>16383) then mantissa:=infinity: exponent:=16383:
        elif (exponent< -16382) then
            # denorm
            exponent := -16383
        fi:
        infmantissa := sign*x*2^(63-exponent);
        if frac(infmantissa) <> 0.5 then mantissa := round(infmantissa)
        else
            mantissa := floor(infmantissa);
            if type(mantissa,odd) then mantissa := mantissa+1 fi;
        fi;
        mantissa := mantissa*2^(-63);
    fi;
    sign,exponent,mantissa;
end:

nearestExt := proc(x)
    local sign, exponent, mantissa;

    sign, exponent, mantissa := ieeedoubleExt(x);
    sign*mantissa*2^(exponent);
end:

ieeehexaExt:= proc(x)
local resultat, sign, exponent, mantissa, t;

    if(x=0) then resultat:=["0000","0000","0000","0000","0000"];
    elif(x=-0) then resultat:=["8000","0000","0000","0000","0000"];
    else
        sign,exponent,mantissa := ieeedoubleExt(x);
        t := 2**80 +  (exponent+16383)*2^64 + mantissa*2^63;
        if (sign=-1) then
            t := t + 2**79;
        fi:
   t := convert(t, hex);
   t:=convert(t, string):

   resultat:=[substring(t, 2..5),
              substring(t, 6..9 ),  substring(t, 10..13 ),
              substring(t, 14..17 ),  substring(t, 18..21 )];

    end if:
    resultat;
end proc:

printDoubleAsShort:=proc(x)
    local ss;
    ss:=ieeehexa(x);
    cat( "DOUBLE_HEX(",
         substring(ss[1], 1..4),  ", " ,
         substring(ss[1], 5..8),  ", " ,
         substring(ss[2], 1..4),  ", " ,
         substring(ss[2], 5..8)) ;
end proc:

printDoubleExtAsShort:=proc(x)
    local ss;
    ss:=ieeehexaExt(x);
    cat( "LDOUBLE_HEX(", ss[1],  ", ", ss[2],  ", ",ss[3],  ", ", ss[4],  ", ", ss[5],  ")");
end proc:

printDoubleAsULL:=proc(x)
    local ss;
    ss:=ieeehexa(x);
    cat( "ULL(", ss[1], ss[2],  ")");
end proc:

printDoubleAsHexInt:=proc(x)
    local ss;
    ss:=ieeehexa(x);
    cat(ss[1], ss[2]);
end proc:

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

#---------------------------------------------------------------------
# Like poly_exact, but the n first coefficients are exactly representable as the sum of two doubles.
# (to actually get the two doubles, use procedure hi_lo)

polyExact2Ext:=proc(P,n)
local deg,i, coef, coef_hi, 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_lo:=hiloExt(coef):
        Q:= Q + coef_hi*x^i:
        if(i<n) then
            Q := Q + coef_lo*x^i:
        fi:
    od:
    return(Q);
end:

printPolyExt := proc(fd,P,n, name_of_poly)
local deg,i, coef, coef_hi, coef_lo;
    convert(Q, polynom):
    deg:=degree(P,x):
    fprintf(fd, " static const long double %s[%d][2] = {\n", name_of_poly, deg+1);
    for i from 0 to deg do
        coef :=coeff(P,x,i):
        coef_hi, coef_lo:=hiloExt(coef):

        fprintf(fd,"{ %1.50eL, ",coef_hi);

        if(i<n) then
            fprintf(fd," %1.50eL},\n",coef_lo);
        else
            fprintf(fd,"0},\n");
        fi:
    od:
    fprintf(fd,"}; \n");
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: work_mgga_c *)
(* prefix:
  mgga_x_mn12_params *params;

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

$define lda_x_params
$include "lda_x.mpl"

omega_x := 2.5:
gamma_x := 0.004:

vx := (rs, z) -> 1/(1 + rs/(omega_x*RS_FACTOR)*(2/(1 + z))^(1/3)):
ux := x -> gamma_x*x^2/(1 + gamma_x*x^2):
wx := t -> (K_FACTOR_C - t)/(K_FACTOR_C + t):

pol1  := t-> params_a_c[ 1] + params_a_c[ 2]*wx(t) + params_a_c[ 3]*wx(t)^2 + params_a_c[ 4]*wx(t)^3
  + params_a_c[ 5]*wx(t)^4 + params_a_c[ 6]*wx(t)^5:
pol2  := t-> params_a_c[ 7] + params_a_c[ 8]*wx(t) + params_a_c[ 9]*wx(t)^2 + params_a_c[10]*wx(t)^3
   + params_a_c[11]*wx(t)^4:
pol3  := t-> params_a_c[12] + params_a_c[13]*wx(t) + params_a_c[14]*wx(t)^2 + params_a_c[15]*wx(t)^3:
pol4  := t-> params_a_c[16] + params_a_c[17]*wx(t) + params_a_c[18]*wx(t)^2:
pol5  := t-> params_a_c[19] + params_a_c[20]*wx(t) + params_a_c[21]*wx(t)^2 + params_a_c[22]*wx(t)^3
   + params_a_c[23]*wx(t)^4:
pol6  := t-> params_a_c[24] + params_a_c[25]*wx(t) + params_a_c[26]*wx(t)^2 + params_a_c[27]*wx(t)^3:
pol7  := t-> params_a_c[28] + params_a_c[29]*wx(t) + params_a_c[30]*wx(t)^2:
pol8  := t-> params_a_c[31] + params_a_c[32]*wx(t) + params_a_c[33]*wx(t)^2 + params_a_c[34]*wx(t)^3:
pol9  := t-> params_a_c[35] + params_a_c[36]*wx(t) + params_a_c[37]*wx(t)^2:
pol10 := t-> params_a_c[38] + params_a_c[39]*wx(t) + params_a_c[40]*wx(t)^2:

FMN12 := (rs, z, x, t) ->
  + pol1(t)
  + pol2(t)*ux(x)
  + pol3(t)*ux(x)^2
  + pol4(t)*ux(x)^3
  + pol5(t)*vx(rs, z)
  + pol6(t)*ux(x)*vx(rs, z)
  + pol7(t)*ux(x)^2*vx(rs, z)
  + pol8(t)*vx(rs, z)^2
  + pol9(t)*ux(x)*vx(rs, z)^2
  + pol10(t)*vx(rs, z)^3:

f_spin := (rs, z, x, t) ->
  lda_x_ax*(1 + z)^(4/3)/rs*FMN12(rs, z, x, t):

f_mn12 := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  f_spin(rs, z, xs0, ts0) + f_spin(rs, -z, xs1, ts1):

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


text:
#
#  Currently, these are not 'tests' in the usual sense, in that they do
# not have correct answers encoded.  They are just meant to all be
# runnable (and are part of the stuff serving as input for the array
# simplification paper we're writing).

kernelopts(assertlevel=2): # be strict on all assertions while testing
kernelopts(opaquemodules=false): # allow testing of internal routines
if not (NewSLO :: `module`) then
  WARNING("loading NewSLO failed");
  `quit`(3);
end if;

with(Hakaru):
with(NewSLO):
with(Loop):

kb := KB:-assert(n :: nonnegint, KB:-empty):
i, kb := KB:-genType(i, HInt(closed_bounds(1..n)), kb):

t1 := Product(idx(x,j),j=1..n):
r1 := unproduct(t1, x, i=1..n, [], `*`, kb, kb);

t2 := 2:
r2 := unproduct(t2, x, i=1..n, [], `*`, kb, kb);

t3 := 2^n:
r3 := unproduct(t3, x, i=1..n, [], `*`, kb, kb);

t4 := Product(2*idx(x,j),j=1..n);
r4 := unproduct(t4, x, i=1..n, [], `*`, kb, kb);

t5 := Product(2+idx(x,j),j=1..n);
r5 := unproduct(t5, x, i=1..n, [], `*`, kb, kb);

t6 := Product(f(idx(x,j)),j=1..n);
r6 := unproduct(t6, x, i=1..n, [], `*`, kb, kb);

t7 := exp(Sum(g(idx(x,k)),k=1..n));
r7 := unproduct(t7, x, i=1..n, [], `*`, kb, kb);

t8 := t6*t7;
r8 := unproduct(t8, x, i=1..n, [], `*`, kb, kb);

t9 := Product(f(idx(x,j+1)),j=0..n-1);
r9 := unproduct(t9, x, i=1..n, [], `*`, kb, kb);

t10 := Product(f(idx(x,j)),j=1..n-1);
r10 := unproduct(t10, x, i=1..n, [], `*`, kb, kb);

# oops, this ones doesn't !!  Going under (above) is
# ok, but over does not.  Neither of the following:
# t10a := Product(f(idx(x,j)),j=1..n+1);
# r10a := unproduct(t10, x, i=1..n, [], `*`, kb, kb);
# t10b := Product(f(idx(x,j)),j=0..n);
# r10b := unproduct(t10, x, i=1..n, [], `*`, kb, kb);

t11 := Product(exp(Sum(f(idx(x,j))+g(idx(x,k)),k=1..j)),j=1..n):
r11 := unproduct(t11, x, i=1..n, [], `*`, kb, kb);

t12 := Product(exp(Sum(f(idx(x,j),j,k)+g(idx(x,k),j,k),k=1..j)),j=1..n):
r12 := unproduct(t12, x, i=1..n, [], `*`, kb, kb);

# this is sort-of buggy as well.
t13 := Product(f(idx(x,5))*g(idx(x,7))*h(j),j=1..n):
r13 := unproduct(t13, x, i=1..n, [], `*`, kb, kb);

# Bug!
# t14 := Product(idx(x,2*j),j=1..n/2):
# r14 := unproduct(t14, x, i=1..n, [], `*`, kb, kb);

t15 := Product(f(idx(x,j),j)*g(idx(x,n+1-j),j)*h(j),j=1..n):
r15 := unproduct(t15, x, i=1..n, [], `*`, kb, kb);

t16 := Product(piecewise(a<b,f(idx(x,j),j),g(idx(x,n+1-j),j)),j=1..n):
r16 := unproduct(t16, x, i=1..n, [], `*`, kb, kb);

t17 := Product(5^idx(x,j)*f(idx(x,j))*g(idx(x,n+1-j)),j=1..n):
r17 := unproduct(t17, x, i=1..n, [], `*`, kb, kb);

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
~ quaternions are represented as a vector(4)
~      r + xi + yj + zk => [r, x, y, z]

function quaternion_rotation(theta, v)
    st = sin(theta/2)
    return [cos(theta/2), st*v[0], st*v[1], st*v[2]]
end

function quaternion_add(x, y)
    return x + y
end

function quaternion_negate(x)
    return -x
end

function quaternion_subtract(x, y)
    return x - y
end

function quaternion_multiply(x, y)
    return [ x[0]*y[0] - x[1]*y[1] - x[2]*y[2] - x[3]*y[3], \
             x[0]*y[1] + x[1]*y[0] + x[2]*y[3] - x[3]*y[2], \
	     x[0]*y[2] + x[2]*y[0] + x[3]*y[1] - x[1]*y[3], \
	     x[0]*y[3] + x[3]*y[0] + x[1]*y[2] - x[2]*y[1] ]
end

function quaternion_conjugate(x)
    return [x[0], -x[1], -x[2], -x[3]]
end

function quaternion_normalize(x)
    return x/sqrt(x*x)
end

function quaternion_inverse(x)
    return quaternion_conjugate(x)/(x*x)
end

function quaternion_divide(x, y)
    return quaternion_multiply(x, quaternion_inverse(y))
end

function quaternion_rotate(x, y)
    return quaternion_multiply(quaternion_multiply(x, y), quaternion_inverse(x))
end

function quaternion_from_matrix(m)
    if m[2,2] < 0 then
        if m[0,0] > m[1,1] then
	    t = 1 + m[0,0] - m[1,1] - m[2,2]
	    q = [m[2,1] - m[1,2], t, m[1,0] + m[0,1], m[0,2] + m[2,0]]
	else
	    t = 1 - m[0,0] + m[1,1] - m[2,2]
	    q = [m[0,2] - m[2,0], m[1,0] + m[0,1], t, m[2,1] + m[1,2]]
	end
    else if m[0,0] < -m[1,1] then
        t = 1 - m[0,0] - m[1,1] + m[2,2]
	q = [m[1,0] - m[0,1], m[0,2] + m[2,0], m[2,1] + m[1,2], t]
    else
        t = 1 + m[0,0] + m[1,1] + m[2,2]
	q = [t, m[2,1] - m[1,2], m[0,2] - m[2,0], m[1,0] - m[0,1]]
    end
    return q/(2*sqrt(t))
end

function quaternion_to_matrix(q)
    return [ 1 - 2*q[2]*q[2] - 2*q[3]*q[3], 2*q[1]*q[2] - 2*q[3]*q[0], 2*q[1]*q[3] + 2*q[2]*q[0]
             2*q[1]*q[2] + 2*q[3]*q[0], 1 - 2*q[1]*q[1] - 2*q[3]*q[3], 2*q[2]*q[3] - 2*q[1]*q[0]
	     2*q[1]*q[3] - 2*q[2]*q[0], 2*q[2]*q[3] + 2*q[1]*q[0], 1 - 2*q[1]*q[1] - 2*q[2]*q[2] ]
end

function quaternion_slerp(xa, ya, t)
    x = quaternion_normalize(xa)
    y = quaternion_normalize(ya)
    dot = x*y
    if dot < 0 then
        x = -x
	dot = -dot
    end
    if dot > 0.9995 then
        return quaternion_normalize((1-t)*x + t*y)
    end
    theta0 = acos(dot)
    theta = t*theta0
    st = sin(theta)
    st0 = sin(theta0)
    s0 = cos(theta) - dot*st/st0
    s1 = st/st0
    return s0*x + s1*y
end

max_stars_repo_name: yugyfoog/MPL

max_issues_repo_name: yugyfoog/MPL

max_forks_repo_name: yugyfoog/MPL


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

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

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

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

gdme_at := (params_a_AA + 3/5*params_a_BB)*2^(1/3)/(X_FACTOR_C*(3*Pi^2)^(2/3)):
gdme_bt := params_a_BB/(X_FACTOR_C*2^(1/3)*(3*Pi^2)^(4/3)):

gdme_f := (x, u, t) -> gdme_at + gdme_bt*((params_a_a^2 - params_a_a + 1/2)*u - 2*t):

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


text:
fun:=proc()
    uses Grid;
    printf("node %d / %d --- %a \n",MyNode(),NumNodes(),[_passed]);
    Barrier();
end proc:
Grid:-Launch(fun,1,2,3);
Grid:-Launch("rand();"):


text:
######################################################################
# This section is about the suboperad of little 1-cubes where the
# subintervals fill the whole interval.

`is_element/one_cubes_prime` := (A::set) -> proc(f)
 local X,n,i;
 global reason;

 if not(`is_element/cubes`(1)(A)(f)) then
  reason := [convert(procname,string),"f is not a 1-cube for A",f,A,reason];
  return false;
 fi;

 X := map(a -> map(op,f[a]),[op(A)]);
 X := sort(X,(u,v) -> evalb(u[1]<v[1]));
 n := nops(A);
 if X[1][1] <> 0 then
  reason := [convert(procname,string),"The first interval does not start at 0",X[1]];
  return false;
 fi;
 if X[n][2] <> 1 then
  reason := [convert(procname,string),"The last interval does not end at 1",X[n]];
  return false;
 fi;

 for i from 1 to n-1 do
  if X[i][2] <> X[i+1][1] then
   reason := [convert(procname,string),"Successive intervals do not touch",X[i],X[i+1]];
   return false;
  fi;
 od;

 return true;
end;

`is_equal/one_cubes_prime` := (A::set) -> proc(f1,f2)
 `is_equal/cubes`(1)(A)(f1,f2);
end:

`is_leq/one_cubes_prime` := NULL;

`random_element/one_cubes_prime` := (A::set) -> proc()
 local d,n,u,i,r,R,xy;

 d := 12;
 n := nops(A);
 u := table();
 u[0] := 0;
 for i from 1 to n do
  u[i] := u[i-1] + rand(1..d)();
 od;
 r := u[n];
 for i from 1 to n do u[i] := u[i]/r; od;

 R := `random_element/ord`(A)();
 xy := table();
 for i from 1 to n do
  xy[R[i]] := [[u[i-1]],[u[i]]];
 od;

 return eval(xy);
end;

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


text:
#!/usr/bin/microperl
#
# roombacmd.pl -- Roomba command-line tool,  part of RoombaComm
#
# http://roombahacking.com/
#
# Copyright (C) 2006, Tod E. Kurt, tod@todbot.com
#
#

# in case we have stty in the current directory
$ENV{'PATH'}="$ENV{PATH}:.";

sub usage() {
    printf "Usage: $0 {serialport} {init|forward|backward|spinleft|spinright|stop|poweroff}\n";
    exit(1);
}

sub roomba_init() {
    # this style stty is for linux
#    system("stty -F $PORT 115200 raw -parenb -parodd cs8 -hupcl -cstopb clocal");
#    stty ttyUSB0 115200 raw -parenb cs8 -hupcl -cstopb clocal;
    open TTYS0, ">$PORT" or die ("Cannot open $PORT");
    printf(TTYS0 "\x80"); sleep 1;
    printf(TTYS0 "\x83");
}
sub roomba_forward() {
    $vel="\x00\x40";
    $rad="\x80\x00";
    open TTYS0, ">$PORT" or die ("Cannot open $PORT");
    printf(TTYS0 "\x89$vel$rad");
}
sub roomba_backward() {
    $vel="\xff\xc0";
    $rad="\x80\x00";
    open TTYS0, ">$PORT" or die ("Cannot open $PORT");
    printf(TTYS0 "\x89$vel$rad");
}
sub roomba_spinleft() {
    $vel="\x00\x40";
    $rad="\x00\x01";
    open TTYS0, ">$PORT" or die ("Cannot open $PORT");
    printf(TTYS0 "\x89$vel$rad");
}
sub roomba_spinright() {
    $vel="\x00\x40";
    $rad="\xff\xff";
    open TTYS0, ">$PORT" or die ("Cannot open $PORT");
    printf(TTYS0 "\x89$vel$rad");
}
sub roomba_stop() {
    $vel="\x00\x00";
    $rad="\x00\x00";
    open TTYS0, ">$PORT" or die ("Cannot open $PORT");
    printf(TTYS0 "\x89$vel$rad");
}

sub roomba_poweroff() {
    open TTYS0, ">$PORT" or die ("Cannot open $PORT");
    printf(TTYS0 "\x85");
}

# If not enough arguments were passed, return
usage() if( @ARGV < 2 );

$PORT = $ARGV[0];
$CMD  = $ARGV[1];

if( $CMD eq 'init' ) {
    roomba_init();
}
elsif( $CMD eq 'forward' ) {
    roomba_forward();
}
elsif( $CMD eq 'backward' ) {
    roomba_backward();
}
elsif( $CMD eq 'spinleft' ) {
    roomba_spinleft();
}
elsif( $CMD eq 'spinright' ) {
    roomba_spinright();
}
elsif( $CMD eq 'stop' ) {
    roomba_stop();
}
elsif( $CMD eq 'poweroff' ) {
    roomba_poweroff();
}
else {
    usage();
}

max_stars_repo_name: timassman/roomba

max_issues_repo_name: timassman/roomba

max_forks_repo_name: timassman/roomba


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_lda *)
(* prefix:
  lda_c_hl_params *params;

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

$ifdef lda_c_vbh_params
params_a_r := [30, 75.0]:
params_a_c := [0.0252, 0.0127]:
$endif

xx := (k, rs) -> rs/params_a_r[k]:
hl := (k, rs) -> -params_a_c[k]*
  ((1 + xx(k, rs)^3)*log(1 + 1/xx(k, rs)) - xx(k, rs)^2 + 1/2*xx(k, rs) - 1/3):

f_hl := (rs, zeta) -> hl(1, rs) + f_zeta(zeta)*(hl(2, rs) - hl(1, rs)):
f    := (rs, zeta) -> f_hl(rs, zeta):


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

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

(* Equation (5) *)
vt84f_f0_orig := s -> 1 - params_a_mu*s^2*exp(-params_a_alpha*s^2)/(1+params_a_mu*s^2) + (1-exp(-params_a_alpha*s^4)) * (s^(-2) - 1) + 5*s^2/3:
(* Since there's a term that looks hairy at small s, do a series expansion up to s^4 *)
vt84f_f0_series := s -> eval(convert(taylor(vt84f_f0_orig(st), st = 0, 5), polynom), st=s):
(* Glue the functions together *)
vt84f_f0 := s-> my_piecewise3(s <= sqrt(DBL_EPSILON), vt84f_f0_series(s), vt84f_f0_orig(m_max(s, sqrt(DBL_EPSILON)))):

(* Convert from x to s *)
vt84f_f := x -> vt84f_f0(X2S*x):

f := (rs, z, xt, xs0, xs1) ->
  gga_kinetic(vt84f_f, rs, z, xs0, xs1):


text:
type $RR <class {@e1 i32, @e2 f32, @e3 f64}>
type $SS <class <$RR> { &method1(agg)agg,
                        &method2(void)void}>
func $foo (
  var %x <$SS>) i32 {
  var %r1 i32
  var  %r2 i32
  var  %r3 i32
  dassign %x 2 ( constval i32 32 )
  # Not correct, but for parse and dump test.
  virtualcallassigned &method2(addrof ptr %x) {
    dassign %r1 1
    dassign %r2 2
    dassign %r3 3
  }
  superclasscallassigned &method2(addrof ptr %x) {
    dassign %r1 1
    dassign %r2 2
    dassign %r3 3
  }
  interfacecallassigned &method2(addrof ptr %x) {
    dassign %r1 1
    dassign %r2 2
    dassign %r3 3
  }
  return ( dread i32 %x 2 ) }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


text:
(*
 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 *)

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

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

gx_c0 :=  0.827411:
gx_c1 := -0.643560:

gx_gx0 := a ->
  + gx_cx0/gx_cx1
  + a*(gx_c0 + gx_c1*a)/(1.0 + (gx_c0 + gx_c1 - 1)*a) * (1 - gx_cx0/gx_cx1):

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

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

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


text:
# Slider_Crank

# 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):
# Ergebnisse von Trigonometrischer Elimination lesen
read sprintf("../codeexport/fourbarprisTE/tmp/kinematic_constraints_maple_inert.m"):
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):

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


text:
program assignboolean;	{sample22}
var
x,y
:
 boolean;

begin
x := true;
y := false;
end.


text:
# P(N) is the set of real polynomials in epsilon of degree less than N

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

`is_element/P` := (N::posint) -> proc(u)
 global reason;

 if not type(u,polynom(realcons,epsilon)) then
  reason := [convert(procname,string),"u is not polynomial in epsilon",u];
  return false;
 fi;

 if degree(u,epsilon) >= N then
  reason := [convert(procname,string),"deg(u) >= N",u,N];
  return false;
 fi;

 return true;
end;

`is_leq/P` := (N::posint) -> proc(u,v)
 local w,d,c;

 w := expand(v-u);
 if (w = 0) then return true; fi;
 d := ldegree(w,epsilon);
 c := coeff(w,epsilon,d);
 return evalb(c >= 0);
end:

`random_element/P` := (N::posint) -> proc(d::posint := 4)
 local r,i;
 r := rand(-d^2..d^2);
 add(r()/d * epsilon^i,i=0..N-1);
end;

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


text:
$ifndef _TEST_H_
$define _TEST_H_
sign:=proc(x)
    if x>0 then
        1
    elif x<0 then
        -1
    else
        0
    end if;
end proc:
$endif


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

$define lda_c_pw_params
$include "lda_c_pw.mpl"

$define lda_x_params
$include "lda_x.mpl"

b98_gamma_x  := 0.11:
b98_gamma_ss := 1.6:
b98_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, u, t) -> 1 - (t - x^2/8 - u/4)/K_FACTOR_C:

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

b98_f := (rs, z, xs0, xs1, us0, us1, ts0, ts1) ->
  + lda_x_spin(rs,  z)
    * b98_g( b98_gamma_x, params_a_c_x, b98_q(xs0, us0, ts0))
  + lda_x_spin(rs, -z)
    * b98_g( b98_gamma_x, params_a_c_x, b98_q(xs1, us1, ts1))
  + lda_stoll_par(f_pw, rs,  z,  1)
    * b98_g(b98_gamma_ss, params_a_c_ss, b98_q(xs0, us0, ts0)) * Fermi_D(xs0, ts0)
  + lda_stoll_par(f_pw, rs, -z, -1)
    * b98_g(b98_gamma_ss, params_a_c_ss, b98_q(xs1, us1, ts1)) * Fermi_D(xs1, ts1)
  + lda_stoll_perp(f_pw, rs, z)
    * b98_g(b98_gamma_ab, params_a_c_ab, (b98_q(xs0, us0, ts0) + b98_q(xs1, us1, ts1))/2):

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


text:
klein_embedding := (x) ->
[(0.1*sin(3*Pi*x[1])+0.1*sin(Pi*x[1])+0.4*cos(Pi*x[1]))*sin(2*Pi*x[2])-0.5*sin(4*Pi*x[1])+sin(2*Pi*x[1]),
  0.2*sin(2*Pi*x[2])*sin(3*Pi*x[1])+2.*cos(2*Pi*x[1])+0.5,
  0.25*cos(2*Pi*x[2])*sin(2*Pi*x[1])+0.4*cos(2*Pi*x[2])];


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

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

$include "mgga_x_m08.mpl"
$include "lda_x_erf.mpl"

f_spin := (rs, z, x, t) ->
  lda_x_ax*(1 + z)^(4/3)/rs * (
    + attenuation_erf(a_cnst*rs/(1 + z)^(1/3)) * m08_f(params_a_a, params_a_b, x, t)
    + (1 - attenuation_erf(a_cnst*rs/(1 + z)^(1/3))) * m08_f(params_a_c, params_a_d, x, t)
  ):

f_m11_l := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  f_spin(rs, z, xs0, ts0) + f_spin(rs, -z, xs1, ts1):

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


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

a := -1:
b := 11.8:
c :=  0.150670:
d := 11.02e-3/RS_FACTOR:

(* I added an abs() to this equation, as otherwise the f_num is
   complex for negative values of z. Of course it is not clear at all
   what was the original intential of Wilson, or if he even considered
   this problem. *)
if evalb(Polarization = "unpol") then
    f_num := z -> a:
else
    (* without the max function we get a float exception for derivatives > 2 *)
    f_num := z -> a*sqrt(1 -  m_max(m_abs(z), 1e-10)^(5/3)):
end if:
f_den := (rs, xt) -> b + c*xt^(51/16) + d*xt^2*rs + rs:

(* Equation (25) *)
f := (rs, z, xt, xs0, xs1) ->
  f_num(z)/f_den(rs, xt):


text:
# NewtonEuler Formulation for Robot based on MDH frames

# Einleitung
# Berechnung der inversen Dynamik in Rückwärtsrekursive
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# fixb ->  Fixed Base Modell der Basis.
# dynamics -> Berechnung der Dynamik
# NewtonEuler->Newton - Euler Bewegungsgleichung
# Linkframe->Berechnung der Newton Euler im Körper-KS(KSi)
# par12 -> Parametersatz 1 (Schwerpunkt als Parameter: SX,SY,SZ) oder Parametersatz 2 (1. und 2. Moment MX,MY,MZ,...)

# Sources
# [GautierKhalil1990] Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots
# [KhalilDombre2002] Modeling, Identification and Control of Robots
# [Ortmaier2014] Vorlesungsskript Robotik I
#
# 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):
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotx":
read "../transformation/proc_troty":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../transformation/proc_trafo_mdh":
# Einstellungen für Code-Export: Optimierungsgrad (2=höchster) und Aktivierung jedes Terms.
codegen_opt := 2:
codeexport_grav := true:
codeexport_corvec := true:
codeexport_cormat := true:
codeexport_inertia := true:
codeexport_inertiaD := true:
codeexport_invdyn := true:
codegen_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 Dynamikgleichungen für %s (Herleitung im Körper-KS)\n", robot_name):
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
read sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
kin_constraints_exist := kin_constraints_exist: # nur zum Abschätzen der Komplexität
;
# Term-Vereinfachungen einstellen
if not assigned(simplify_options) or simplify_options(9)=-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
    if NJ <= 5 then # nur bei einfachen Systemen Term-Vereinfachung durchführen
      use_simplify := 1:
    else
      use_simplify := 0: # Annahme: Dauert zu lange bzw. zu rechenintensiv (nicht geprüft)
    end if:
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(9): # neunter Eintrag ist für Dynamik
end if:

read "../helper/proc_index_symmat2vector":
read "../helper/proc_symmat2vector":
# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name, base_method_name):
Trf := Trf:
Trf_c := Trf_c:
# Ergebnisse der Geschwindigkeit laden
read sprintf("../codeexport/%s/tmp/velocity_linkframe_floatb_%s_maple.m", robot_name, base_method_name):
omega_i_i:= omega_i_i:
rD_i_i:= rD_i_i:
# Zeitableitungen der MDH-Drehwinkel laden.
# Die Berechnung soll nur an einer Stelle erfolgen. Siehe robot_tree_velocity_mdh_angles.mw.
read sprintf("../codeexport/%s/tmp/velocity_mdh_angles_maple.m", robot_name):
thetaD := thetaD:
dD :=dD:
# Zeitableitungen der Geschwindigkeit laden.
# Die Berechnung soll nur an einer Stelle erfolgen. Siehe robot_tree_acceleration_mdh_angles.mw.
read sprintf("../codeexport/%s/tmp/acceleration_mdh_angles_maple.m", robot_name):
thetaDD :=thetaDD:
dDD := dDD:
# Ergebnisse der Geschwindigkeit laden. TODO: Warum doppelt?
read sprintf("../codeexport/%s/tmp/velocity_linkframe_floatb_%s_maple.m", robot_name, base_method_name):
omega_i_i:= omega_i_i:
rD_i_i:= rD_i_i:
# Ergebnisse der Beschleunigung laden
accfile := sprintf("../codeexport/%s/tmp/acceleration_linkframe_floatb_%s_maple.m", robot_name, base_method_name):
if FileTools[Exists](accfile) then
  read accfile:
else
  printf("%s. Schwerpunktsbeschleunigung wurde nicht berechnet. Abbruch der Newton-Euler-Berechnung.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  quit: # Funktioniert in GUI nicht richtig...
  robot_name := "": # ...Daher auch Löschung des Roboternamens.
end if:
omegaD_i_i:= omegaD_i_i:
rDD_i_i:= rDD_i_i:
rDD_i_Si:= rDD_i_Si:
# 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:
# Schalter zur Auswahl der unterschiedlichen Terme, die exportiert werden sollen.

# Newton-Euler formulation(mit Funktion)
OutputNewtonEuler:= NewtonEuler(f_i_i,m_i_i):
# Backward recursive computation
# Force exerted on link i by link i-1(unbekannte Kraft zwischen Teilkörper n-1 und Teilkörper n)
f_i_i:= Matrix(3,NJ+1):
# Moment exerted on link i by link i-1(unbekannte Moment zwischen Teilkörper n-1 und Teilkörper n)
m_i_i:=Matrix(3,NJ+1):

# Schnittmomente/Kräfte
tau_B := Matrix(6,1):
tau_J := Matrix(NQJ,1):
tau := Matrix(6+NQJ,1):

I_i_i_m:= Matrix(3,3,NL):
for i from 1 to NL do
  I_i_i_m(1..3,1..3,i):=Matrix(3,3):
 end do:

for i from 1 to NL do
  I_i_i_m[1,1,i]:=I_i_i[1,i]:
  I_i_i_m[1,2,i]:=I_i_i[2,i]:
  I_i_i_m[1,3,i]:=I_i_i[3,i]:
  I_i_i_m[2,1,i]:=I_i_i[2,i]:
  I_i_i_m[2,2,i]:=I_i_i[4,i]:
  I_i_i_m[2,3,i]:=I_i_i[5,i]:
  I_i_i_m[3,1,i]:=I_i_i[3,i]:
  I_i_i_m[3,2,i]:=I_i_i[5,i]:
  I_i_i_m[3,3,i]:=I_i_i[6,i]:
end do:

F_i:= Matrix(3,NL):
M_i:= Matrix(3,NL):
for i from 1 to NL do
  # External forces and moment on link i
  #Betrachte dazu nur die durch Gelenke angetriebenen Körper, nicht die Basis
  # [Khali2002](9.5.3) (9.88) & (9.89)
  F_i(1..3,i):= M(i,1)*rDD_i_i(1..3,i)+CrossProduct(omegaD_i_i(1..3,i),mr_i_i_Si(1..3,i))+ CrossProduct(omega_i_i(1..3,i), CrossProduct(omega_i_i(1..3,i), mr_i_i_Si(1..3,i))):
  M_i(1..3,i):= I_i_i_m(1..3,1..3,i).omegaD_i_i(1..3,i)+ CrossProduct(omega_i_i(1..3,i),(I_i_i_m(1..3,1..3,i).omega_i_i(1..3,i))) + CrossProduct(mr_i_i_Si(1..3,i),rDD_i_i(1..3,i)):#CrossProduct(r,F_i(1..3,i)):####
  # Terme vereinfachen
  if use_simplify>=1 then
    tmp_t1:=time():
    tmp_l11 := length(F_i(1..3,i)):
    tmp_l12 := length(M_i(1..3,i)):
    F_i(1..3,i) := simplify2(F_i(1..3,i)):
    M_i(1..3,i) := simplify2(M_i(1..3,i)):
    tmp_l21 := length(F_i(1..3,i)):
    tmp_l22 := length(M_i(1..3,i)):
    tmp_t2:=time():
    printf("%s: Terme für inneres Kraft/Moment %d vereinfacht. Länge: %d->%d / %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i, tmp_l11, tmp_l21, tmp_l12, tmp_l22, tmp_t2-tmp_t1):
  end if:
end do :

# Forces and moment exerted on link i by link i-1

for i from NL by -1 to 1 do
  f_i_i_part:= F_i(1..3,i):
  m_i_i_part:= M_i(1..3,i):# + CrossProduct(r_i_i_Si(1..3,i),(F_i(1..3,i))):

  I_nf := Matrix([ListTools[SearchAll](i-1,convert(v,list))]):

  for tmp from 1 to ColumnDimension(I_nf) do
    I_nf(tmp) := I_nf(tmp)+1:
  end do:

  if ColumnDimension(I_nf) <> 0 then
    for tmp from 1 to RowDimension(v) do
      j := I_nf(tmp):

      R_i_j := Matrix(Trf(1..3,1..3,j-1)):
      r_i_i_j := Matrix(Trf(1 .. 3, 4, j-1)):
      f_j_j := f_i_i(1..3,j):
      m_j_j := m_i_i(1..3,j):

      f_i_i_part := f_i_i_part + R_i_j.f_j_j:
      m_i_i_part := m_i_i_part + R_i_j.m_j_j + CrossProduct(r_i_i_j, R_i_j.f_j_j):

      if tmp = ColumnDimension(I_nf) then
        break: #Abbruch. Alle Nachfolger untersucht.
      end if:
    end do:
  end if:

  f_i_i(1..3,i) := f_i_i_part;
  m_i_i(1..3,i) := m_i_i_part;
  # Terme vereinfachen
  if use_simplify>=1 then
    tmp_t1:=time():
    tmp_l11 := length(f_i_i(1..3,i)):
    tmp_l12 := length(m_i_i(1..3,i)):
    f_i_i(1..3,i) := simplify2(f_i_i(1..3,i)):
    m_i_i(1..3,i) := simplify2(m_i_i(1..3,i)):
    tmp_l21 := length(f_i_i(1..3,i)):
    tmp_l22 := length(m_i_i(1..3,i)):
    tmp_t2:=time():
    printf("%s: Rekursive Terme für Kraft/Moment %d vereinfacht. Länge: %d->%d / %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i-1, tmp_l11, tmp_l21, tmp_l12, tmp_l22, tmp_t2-tmp_t1):
  end if:
end do:

for i from 2 to NL do # Schleife über Anzahl der Körper
  if sigma(i-1) <> 2 then
    if sigma(i-1) = 0 then # Drehgelenk
      tau_J(i-1,1) := <0,0,1>.m_i_i(1..3,i):
    else # Schubgelenk
      tau_J(i-1,1) := <0,0,1>.f_i_i(1..3,i):
    end if:
  end if:
end do:
# Unbekannten Kräfte und Momente zwischen den Teilkörpern bestimmen
tau_B [1..3,1]:= Transpose(Trf_c(1..3, 1..3, 1)).f_i_i(1..3,1) :
tau_B [4..6,1]:= Transpose(T_basevel).Trf_c(1..3, 1..3, 1).m_i_i(1..3,1):

tau[1..6,1]:= tau_B:
tau[7..6+NQJ,1]:= tau_J:
# Export
f_i_i := convert_t_s(f_i_i):
m_i_i := convert_t_s(m_i_i):
tau_J := convert_t_s(tau_J):
tau_B := convert_t_s(tau_B):
tau := convert_t_s(tau):
# Floating Base
f_i_i_floatb := copy(f_i_i):
m_i_i_floatb := copy(m_i_i):
tau_J_floatb := copy(tau_J):
tau_B_floatb := copy(tau_B):
tau_floatb := copy(tau):
for i from 1 to 6 do
  if i < 4 then
    var := VD_base_s[i,1] - g_world[i,1];
  else
    var := VD_base_s[i,1]:
  end if:
  f_i_i_floatb := subs({VD_base_s[i,1]=var},f_i_i_floatb):
  m_i_i_floatb := subs({VD_base_s[i,1]=var},m_i_i_floatb):
  tau_J_floatb := subs({VD_base_s[i,1]=var},tau_J_floatb):
  tau_B_floatb := subs({VD_base_s[i,1]=var},tau_B_floatb):
  tau_floatb := subs({VD_base_s[i,1]=var},tau_floatb):
end do:
# Maple Export # TODO: Wieso Bedingung "twist" gewählt?
if not(base_method_name="twist") then
  save  f_i_i_floatb,m_i_i_floatb,tau_J_floatb, tau_B_floatb, tau_floatb,sprintf("../codeexport/%s/tmp/invdyn_%s_NewtonEuler_linkframe_maple.m", robot_name, base_method_name):
  printf("%s. Maple-Ausdrücke exportiert.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
end if:
# Matlab Export
if codegen_act and codeexport_invdyn and not(base_method_name="twist") then
  MatlabExport(convert_t_s(f_i_i_floatb), sprintf("../codeexport/%s/tmp/invdyn_floatb_%s_NewtonEuler_linkframe_f_i_i_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
  MatlabExport(convert_t_s(m_i_i_floatb), sprintf("../codeexport/%s/tmp/invdyn_floatb_%s_NewtonEuler_linkframe_m_i_i_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
  MatlabExport(convert_t_s(tau_J_floatb), sprintf("../codeexport/%s/tmp/invdyn_floatb_%s_NewtonEuler_linkframe_tauJ_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
  MatlabExport(convert_t_s(tau_B_floatb), sprintf("../codeexport/%s/tmp/invdyn_floatb_%s_NewtonEuler_linkframe_tauB_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
  MatlabExport(convert_t_s(tau_floatb), sprintf("../codeexport/%s/tmp/invdyn_floatb_%s_NewtonEuler_linkframe_tauJB_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Fixed Base
for i from 1 to NQB do
  f_i_i := subs({X_base_s[i,1]=0},f_i_i):
  m_i_i := subs({X_base_s[i,1]=0},m_i_i):
  tau_J := subs({X_base_s[i,1]=0},tau_J):
  tau_B := subs({X_base_s[i,1]=0},tau_B):
  tau := subs({X_base_s[i,1]=0},tau):
end do:
for i from 1 to 6 do
  f_i_i := subs({V_base_s[i,1]=0},f_i_i):
  m_i_i := subs({V_base_s[i,1]=0},m_i_i):
  tau_J := subs({V_base_s[i,1]=0},tau_J):
  tau_B := subs({V_base_s[i,1]=0},tau_B):
  tau := subs({V_base_s[i,1]=0},tau):
  if i < 4 then
    var := -g_world[i,1];
  else
    var := 0:
  end if:
  f_i_i := subs({VD_base_s[i,1]=var},f_i_i):
  m_i_i := subs({VD_base_s[i,1]=var},m_i_i):
  tau_J := subs({VD_base_s[i,1]=var},tau_J):
  tau_B := subs({VD_base_s[i,1]=var},tau_B):
  tau := subs({VD_base_s[i,1]=var},tau):
end do:
# Maple Export
save  f_i_i,m_i_i,tau_J, tau_B, tau,sprintf("../codeexport/%s/tmp/invdyn_%s_NewtonEuler_linkframe_par%d_maple.m", robot_name, base_method_name, codegen_dynpar):
printf("%s. Maple-Ausdrücke exportiert.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
# Matlab Export
if codegen_act and codeexport_invdyn  then
  MatlabExport(f_i_i, sprintf("../codeexport/%s/tmp/invdyn_fixb_NewtonEuler_linkframe_f_i_i_par%d_matlab.m", robot_name, codegen_dynpar), codegen_opt):
  MatlabExport(m_i_i, sprintf("../codeexport/%s/tmp/invdyn_fixb_NewtonEuler_linkframe_m_i_i_par%d_matlab.m", robot_name, codegen_dynpar), codegen_opt):
  MatlabExport(tau_J, sprintf("../codeexport/%s/tmp/invdyn_fixb_NewtonEuler_linkframe_tauJ_par%d_matlab.m", robot_name, codegen_dynpar), codegen_opt):
  MatlabExport(tau_B, sprintf("../codeexport/%s/tmp/invdyn_fixb_NewtonEuler_linkframe_tauB_par%d_matlab.m", robot_name, codegen_dynpar), codegen_opt):
  MatlabExport(tau, sprintf("../codeexport/%s/tmp/invdyn_fixb_NewtonEuler_linkframe_tauJB_par%d_matlab.m", robot_name, codegen_dynpar), 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:
`gamma/permutations` := (J) -> proc(s,t)
 local k,i,m,tt;
 k := nops(J);
 i := `+`(op(J));
 tt := [];
 m := 0;
 for i from 1 to k do
  tt := [op(tt),t[i] +~ m];
  m := m + J[i];
 od:
 return [seq(op(tt[s[i]]),i=1..k)];
end:

`eta/permutations` := [1];

`circ/permutations` := (i,m,n) -> proc(s,t)
 local u,p;
 u := NULL;
 for p from 1 to m do
  if s[p] < i then
   u := u,s[p]
  elif s[p] = i then
   u := u , op(t +~ (i - 1));
  else
   u := u,s[p]+n-1;
  fi;
 od;
 return [u];
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: mgga_exc *)
(* prefix:
  mgga_x_msb_params *params;

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

$include "mgga_x_ms.mpl"

(* eq (5) in the paper *)
msb_beta := (t, x) -> ms_alpha(t, x)*K_FACTOR_C/(t + K_FACTOR_C):

(* eq (14) in the supplement *)
msb_fa := b -> (1 - (2*b)^2)^3 / (1 + (2*b)^3 + params_a_b*(2*b)^6):

msb_f := (x, u, t) -> ms_f0(X2S^2*x^2, 0) + \
  msb_fa(msb_beta(t,x))*(ms_f0(X2S^2*x^2, params_a_c) - ms_f0(X2S^2*x^2, 0)):

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


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: lda_exc *)
(* prefix:
  lda_x_sloc_params *params;

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

# https://onlinelibrary.wiley.com/doi/full/10.1002/qua.25312

f_sloc := (rs, z) ->
  -params_a_a/(2*(params_a_b + 1)) * n_total(rs)^params_a_b *
  (opz_pow_n(z,params_a_b + 1) + opz_pow_n(-z,params_a_b + 1)):

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


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

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

ml1_C := 6.187335:
ml1_b := [2.763169, 1.757515, 1.741397, 0.568985, 1.572202, 1.885389]:

ml1_alpha := z -> params_a_fc*((1 + z)^params_a_q + (1 - z)^params_a_q):
ml1_beta  := z -> (1 - z^2)^(1/3)/((1 + z)^(1/3) + (1 - z)^(1/3)):

(* From the paper: "Note that the antiparailel-spin correlation length
   diverges when the spin-polarization parameter tends to 1", which means
   that Q diverges for a ferromagnetic density *)
ml1_k := (rs, z) -> ml1_C*n_total(rs)^(1/3) * ml1_alpha(z)*ml1_beta(z):

(* Eq. 32 *)
ml1_Q := (rs, z) ->
  - ml1_b[1]/(1 + ml1_b[2]*ml1_k(rs, z))
  + ml1_b[3]*log(1 + ml1_b[4]/ml1_k(rs, z))/ml1_k(rs, z)
  + ml1_b[5]/ml1_k(rs, z)
  - ml1_b[6]/ml1_k(rs, z)^2:

(* screen for small spin densities to avoid divergences in the
  potentials.  Note that beta is zero for any polarized density and
  the whole expression for alpha*beta is symmetric in z.  Note also
  that in the expression for Q one divides by k that is zero for
  ferromagnetic densities. *)

(* there is a factor of 1/2 wrong in Eq. 31 as explained in the Erratum *)
(* With the formula below we can reproduce exactly the values of Table I.
   Note that in the Erratum the authors afirm that all the results are correct,
   and only the formulas had misspells. *)
ml1_f := (rs, z) -> n_total(rs) *
  my_piecewise3(1 - abs(z) <= p_a_zeta_threshold, 0, (1 - z^2)/4 * ml1_Q(rs, z_thr(z))):

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


text:
# We identify R^4 with the quaternions by the rule
# [x1,x2,x3,x4] |-> x1 i + x2 j + x3 k + x4.

`mu/H` := proc(a,b)
 if nargs = 0 then
  return [0,0,0,1];
 elif nargs = 1 then
  return args[1];
 elif nargs > 2 then
  return `mu/H`(a,`mu/H`(args[2..-1]));
 else
  return
  [  a[1]*b[4] + a[2]*b[3] - a[3]*b[2] + a[4]*b[1] ,
   - a[1]*b[3] + a[2]*b[4] + a[3]*b[1] + a[4]*b[2] ,
     a[1]*b[2] - a[2]*b[1] + a[3]*b[4] + a[4]*b[3] ,
   - a[1]*b[1] - a[2]*b[2] - a[3]*b[3] + a[4]*b[4]
  ];
 fi;
end:

`conj/H` := proc(a)
 [-a[1],-a[2],-a[3],a[4]];
end:

`square_norm/H` := proc(a) local i; return add(a[i]^2,i=1..4); end:

`norm/H` := proc(a) local i; return sqrt(add(a[i]^2,i=1..4)); end:

`inv/H` := (a) -> `conj/H`(a) /~ `square_norm/H`(a);

`phi/angles/H` := (t) ->
 [cos(t[1])*cos(t[2]),cos(t[1])*sin(t[2]),sin(t[1])*cos(t[3]),sin(t[1])*sin(t[3])];

# This is essentially a |-> a * i * (conjugate of a)
hopf_map := (a) ->
 [a[1]^2-a[2]^2-a[3]^2+a[4]^2,
  2*a[1]*a[2]+2*a[3]*a[4],
  2*a[1]*a[3]-2*a[2]*a[4]];

`rotation_matrix/H` := (a) -> Matrix(
 [[a[1]^2-a[2]^2-a[3]^2+a[4]^2,
   2*a[1]*a[2]-2*a[3]*a[4],
   2*a[1]*a[3]+2*a[2]*a[4]],
  [2*a[1]*a[2]+2*a[3]*a[4],
   -a[1]^2+a[2]^2-a[3]^2+a[4]^2,
   -2*a[1]*a[4]+2*a[2]*a[3]],
  [2*a[1]*a[3]-2*a[2]*a[4],
   2*a[1]*a[4]+2*a[2]*a[3],
   -a[1]^2-a[2]^2+a[3]^2+a[4]^2]]
  ) / (a[1]^2+a[2]^2+a[3]^2+a[4]^2);


text:
with(DocumentTools):
SetProperty("LogAreaSIAN", value, ""):
SetProperty("LogAreaSE", value, ""):
SetProperty("LogAreaME", value, ""):

SetProperty("GlobalParams1", expression, NULL):
SetProperty("LocalParams1", expression, NULL):
SetProperty("NoIDParams1", expression, NULL):

SetProperty("being_refined", caption, ""):

SetProperty("Bound", expression, NULL):
SetProperty("MultiFunctions", expression, NULL):
SetProperty("SingleFunctions", expression, NULL):

SetProperty("RunningTimeSIAN", value, ""):
SetProperty("RunningTimeMulti", value, ""):
SetProperty("RunningTimeSingle", value, ""):

SetProperty("Parameters", value, ""):

max_stars_repo_name: iliailmer/sian-web-app

max_issues_repo_name: iliailmer/sian-web-app

max_forks_repo_name: iliailmer/sian-web-app