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

interface(quiet=true):

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

printPolynomialIntoFile := proc(fd,s,p)
local i, hi, mi, lo:
for i from 0 to degree(p(x),x) do
	(hi,mi,lo) := hi_mi_lo(coeff(p(x),x,i)):
	if ((abs(hi) = 1.0) and (mi = 0) and (lo = 0)) then
		printf(
		"Coefficient %d of the polynomial is exactly %f and will not be stored in the table\n",i,hi):
	else
	if ((abs(hi) = 0.5) and (mi = 0) and (lo = 0)) then
		printf(
		"Coefficient %d of the polynomial is exactly %f and will not be stored in the table\n",i,hi):
	else
	if (hi <> 0) then
 		fprintf(fd,"#define %s%dh %1.50e\n",s,i,hi):
	end if:
	if (mi <> 0) then
	 	fprintf(fd,"#define %s%dm %1.50e\n",s,i,mi):
	end if:
	if (lo <> 0) then
	 	fprintf(fd,"#define %s%dl %1.50e\n",s,i,lo):
	end if:
	end if:
	end if:
od:
end proc:

# First, we compute special values

ReturnXBound := convert((ieeehexa(2^(-54)))[1],decimal,hex):

Largest := 2^(1023) * ((2^(53) - 1) / 2^(52)):

Smallest := 2^(-1023) * 1 * 2^(-51):

OverflowBound := nearest(log(Largest + 1)):

MinusOneBound := nearest(log(2^(-54))):

SimpleOverflowBound := convert(ieeehexa(OverflowBound)[1],decimal,hex):

DirectIntervalBound := convert((ieeehexa(0.25))[1],decimal,hex):

MinusOnePlusOneUlp := -1 + 2^(-53): # Attention: it's 2^(-53) because we are at a binade boundary

# Second, we have the computation of the values for the direct interval

# The function, that we approximate is

directF := unapply(exp(x) - 1,x):

# The domain is

directA := -2^(-5):
directB := 2^(-5):

# The polynomials are

quickDirectpoly := X -> X+1/2*X^2+(3360682229480701/1180591620717411303424*X^6+3660136839517697/147573952589676412928*X^5+7320130809407439/36893488147419103232*X^4+3202559734508631/2305843009213693952*X^3+4803839602572223/576460752303423488*X^2+6004799503160665/144115188075855872*X+6004799503160661/36028797018963968)*X^3:

accuDirectpoly := X -> X+1/2*X^2+(3786738884990361/4951760157141521099596496896*X^12+7100145222887513/618970019642690137449562112*X^11+6212541673969101/38685626227668133590597632*X^10+5047690109993399/2417851639229258349412352*X^9+3785767582868083/151115727451828646838272*X^8+5205430426443615/18889465931478580854784*X^7+29303968161043118891149009244865/10633823966279326983230456482242756608*X^6+65933928362347017505024866986963/2658455991569831745807614120560689152*X^5+65933928362347017505149159875899/332306998946228968225951765070086144*X^4+28846093658526820158502757550845/20769187434139310514121985316880384*X^3+21634570243895115118877068038417/2596148429267413814265248164610048*X^2+27043212804868893898596335048021/649037107316853453566312041152512*X+243583606221817153033947472119380503276473908509/1461501637330902918203684832716283019655932542976)*X^3:

# Truncate the quick phase direct interval polynomial to degree specialDegree
# for special interval |x| <= specialBound (speed-up)

specialDegree := 5:
specialBound := 2^(-12):

specialPoly := unapply(sum(coeff(quickDirectpoly(x),x,i) * x^i,i=0..specialDegree),x):

printf("Special polynomial is the direct polynomial truncated to degree %d used in |x| < 2^(%f)\n",
	specialDegree, evalf(log[2](specialBound))):

# Compute the relative errors

errDirectQuick := numapprox[infnorm](quickDirectpoly(x)/directF(x) -1,x=directA..directB):
errDirectAccu := numapprox[infnorm](accuDirectpoly(x)/directF(x) -1,x=directA..directB):

errSpecialPoly := numapprox[infnorm](specialPoly(x)/directF(x) -1,x=-specialBound..specialBound):

errDirectAccuSpecial := numapprox[infnorm](accuDirectpoly(x)/directF(x) -1,x=2^(-12)..2^(-12)):

printf("The relative approximation error of the direct interval quick polynomial is 2^(%f)\n",
	evalf(log[2](abs(errDirectQuick)))):
printf("The relative approximation error of the direct interval accurate polynomial is 2^(%f)\n",
	evalf(log[2](abs(errDirectAccu)))):
printf("The relative approximation error of the special interval special polynomial is 2^(%f)\n",
	evalf(log[2](abs(errSpecialPoly)))):
printf("The relative approximation error of the direct interval accurate polynomial in special domain is 2^(%f)\n",
	evalf(log[2](abs(errDirectAccuSpecial)))):

# Third, we have the computation of the values for the common interval

# The function, that we approximate is

commonF := unapply(exp(x),x):

# The domain is

commonA := -log(2)*2^(-12) * (1/2 + 2^(-19)):
commonB := log(2)*2^(-12) * (1/2 + 2^(-19)):

quickCommonpoly := X -> 1+X+1/2*X^2+(6004799504593679/144115188075855872*X+6004799504235425/36028797018963968)*X^3:

accuCommonpoly := X -> 1+X+1/2*X^2+(3660068549402285/18446744073709551616*X^4+6405119471061623/4611686018427387904*X^3+4803839602528529/576460752303423488*X^2+54086425609737787796676993069745/1298074214633706907132624082305024*X+54086425609737787797192670135537/324518553658426726783156020576256)*X^3:

# Compute the relative errors

errCommonQuick := numapprox[infnorm](quickCommonpoly(x)/commonF(x) -1,x=commonA..commonB):
errCommonAccu := numapprox[infnorm](accuCommonpoly(x)/commonF(x) -1,x=commonA..commonB):

printf("The relative approximation error of the common interval quick polynomial is 2^(%f)\n",
	evalf(log[2](abs(errCommonQuick)))):
printf("The relative approximation error of the common interval accurate polynomial is 2^(%f)\n",
	evalf(log[2](abs(errCommonAccu)))):

epsilonApproxRmAccurate := numapprox[infnorm]( ((1+x)/(exp(x)))-1, x=commonA*2^(-52)..commonB*2^(-52)):
epsilonApproxRlAccurate := numapprox[infnorm]( ((1+x)/(exp(x)))-1, x=commonA*2^(-105)..commonB*2^(-105)):

printf("The approximation rel error for approximating exp(rm) by 1 + rm is 2^(%2f)\n",
	log2(abs(epsilonApproxRmAccurate))):
printf("The approximation rel error for approximating exp(rl) by 1 + rl is 2^(%2f)\n",
	log2(abs(epsilonApproxRlAccurate))):

# Compute the constants for argument reduction and the tables in the common path

MsLog2Div2L := evalf(-log(2)/(2^(12))):

msLog2Div2Lh, msLog2Div2Lm, msLog2Div2Ll := hi_mi_lo(MsLog2Div2L):

epsMsLog2Div2L := evalf(abs(((msLog2Div2Lh + msLog2Div2Lm + msLog2Div2Ll) - MsLog2Div2L)/MsLog2Div2L)):
epsDDMsLog2Div2L := evalf(abs(((msLog2Div2Lh + msLog2Div2Lm) - MsLog2Div2L)/MsLog2Div2L)):

printf("The error made by storing MsLog2Div2L as a double-double is 2^(%f)\n",log[2](epsDDMsLog2Div2L)):
printf("The error made by storing MsLog2Div2L as a triple-double is 2^(%f)\n",log[2](epsMsLog2Div2L)):

gap := -floor(-log[2](abs(msLog2Div2Lm/msLog2Div2Lh))):

printf("Information: |msLog2Div2Lm| <= 2^(%f) * |msLog2Div2Lh|\n",gap):

log2InvMult2L := nearest(2^(12) / (log(2))):

shiftConst := 2^(52) + 2^(51):

indexmask1 := 2^((12)/2) - 1:
indexmask2 := indexmask1 * 2^((12)/2):

for i from 0 to 2^(12/2) - 1 do
	twoPowerIndex1hi[i], twoPowerIndex1mi[i], twoPowerIndex1lo[i] := hi_mi_lo(evalf(2^(i/(2^12)))):
	twoPowerIndex2hi[i], twoPowerIndex2mi[i], twoPowerIndex2lo[i] := hi_mi_lo(evalf(2^(i/(2^(12/2))))):
od:

# Estimate the error of the two quick phases

# ATTENTION: C EST PIFOMETRIQUE POUR L INSTANT

epsQuickDirectOverall := 2^(-62):
epsQuickCommonOverall := 2^(-62):

# Write the tables

printf("Write tables...\n"):

filename:="TEMPEXPM1/expm1.h":
fd:=fopen(filename, WRITE, TEXT):

fprintf(fd, "#include \"crlibm.h\"\n#include \"crlibm_private.h\"\n"):

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

fprintf(fd, "\#define log2InvMult2L %1.50e\n",log2InvMult2L):
fprintf(fd, "\#define msLog2Div2Lh %1.50e\n",msLog2Div2Lh):
fprintf(fd, "\#define msLog2Div2Lm %1.50e\n",msLog2Div2Lm):
fprintf(fd, "\#define msLog2Div2Ll %1.50e\n",msLog2Div2Ll):
fprintf(fd, "\#define shiftConst %1.50e\n",shiftConst):
fprintf(fd, "\#define INDEXMASK1 0x%08x\n",indexmask1):
fprintf(fd, "\#define INDEXMASK2 0x%08x\n",indexmask2):
fprintf(fd, "\#define RETURNXBOUND 0x%08x\n",ReturnXBound):
fprintf(fd, "\#define OVERFLOWBOUND %1.50e\n",OverflowBound):
fprintf(fd, "\#define LARGEST %1.50e\n",Largest):
fprintf(fd, "\#define SMALLEST %1.50e\n",Smallest):
fprintf(fd, "\#define MINUSONEBOUND %1.50e\n",MinusOneBound):
fprintf(fd, "\#define SIMPLEOVERFLOWBOUND 0x%08x\n",SimpleOverflowBound):
fprintf(fd, "\#define DIRECTINTERVALBOUND 0x%08x\n",DirectIntervalBound):
fprintf(fd, "\#define SPECIALINTERVALBOUND 0x%08x\n",convert((ieeehexa(specialBound))[1],decimal,hex)):
fprintf(fd, "\#define ROUNDCSTDIRECTRN %1.50e\n",compute_rn_constant(epsQuickDirectOverall)):
fprintf(fd, "\#define ROUNDCSTDIRECTRD %1.50e\n",epsQuickDirectOverall):
fprintf(fd, "\#define ROUNDCSTCOMMONRN %1.50e\n",compute_rn_constant(epsQuickCommonOverall)):
fprintf(fd, "\#define ROUNDCSTCOMMONRD %1.50e\n",epsQuickCommonOverall):
fprintf(fd, "\#define MINUSONEPLUSONEULP %1.50e\n",MinusOnePlusOneUlp):

fprintf(fd,"\n\n"):

printPolynomialIntoFile(fd,"quickDirectpolyC",quickDirectpoly):
fprintf(fd,"\n"):
printPolynomialIntoFile(fd,"accuDirectpolyC",accuDirectpoly):
fprintf(fd,"\n"):
printPolynomialIntoFile(fd,"quickCommonpolyC",quickCommonpoly):
fprintf(fd,"\n"):
printPolynomialIntoFile(fd,"accuCommonpolyC",accuCommonpoly):

fprintf(fd,"\n\n"):

# Print the tables
fprintf(fd, "typedef struct tPi_t_tag {double hi; double mi; double lo;} tPi_t;  \n"):
fprintf(fd, "static const tPi_t twoPowerIndex1[%d] = {\n", 2^(12/2)):
for i from 0 to 2^(12/2)-1 do
      fprintf(fd, "  { \n"):
      fprintf(fd, "    %1.50e, /* twoPowerIndex1hi[%d] */ \n", twoPowerIndex1hi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex1mi[%d] */ \n", twoPowerIndex1mi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex1lo[%d] */ \n", twoPowerIndex1lo[i], i):
      fprintf(fd, "  } "):
      if(i<2^(12/2)-1) then  fprintf(fd, ", \n"): fi
od:
fprintf(fd, "}; \n \n"):
fprintf(fd, "static const tPi_t twoPowerIndex2[%d] = {\n", 2^(12/2)):
for i from 0 to 2^(12/2)-1 do
      fprintf(fd, "  { \n"):
      fprintf(fd, "    %1.50e, /* twoPowerIndex2hi[%d] */ \n", twoPowerIndex2hi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex2mi[%d] */ \n", twoPowerIndex2mi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex2lo[%d] */ \n", twoPowerIndex2lo[i], i):
      fprintf(fd, "  } "):
      if(i<2^(12/2)-1) then  fprintf(fd, ", \n"): fi
od:
fprintf(fd, "}; \n \n"):

fprintf(fd, "\n\n"):

fclose(fd):

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

max_stars_repo_name: squarePenguin/parvsl

max_issues_repo_name: squarePenguin/parvsl

max_forks_repo_name: squarePenguin/parvsl


text:
program assignchar;
var x, y : char;
begin
    x := 'X';
    y := 'Y'
end.

ext: mpl

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing

max_forks_repo_name: Hiroya-W/lang-processing


text:
#int fact (int n) {
#  if(n != 1)
#    return foo(n - 1);
#  else return 1;
#}
func $foo varargs ( var %i i32, ... ) i32 {
  return (dread i32 %i) }
func $foo1 varargs ( ... ) i32 {
  return (constval i32 99) }

func $fact (
 var %n i32) i32 {
 if (ne i32 i32 (dread i32 %n, constval i32 1)) {
    call &foo(
      sub i32 (dread i32 %n, constval i32 1))
    return (regread i32 %%retval)}
   else {
     return (constval i32 1) } }
 # 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:
######################################################################

`is_element/equiv` := (A::set) -> proc(R)
 global reason;

 if not `is_element/autorel`(A)(R) then
  reason := [convert(procname,string),"R is not a relation on A",R,A,reason];
  return false;
 fi;

 if not(`is_reflexive/autorel`(A)(R)) then
  reason := [convert(procname,string),"R is not reflexive",R];
  return false;
 fi;

 if not(`is_symmetric/autorel`(A)(R)) then
  reason := [convert(procname,string),"R is not symmetric",R];
  return false;
 fi;

 if not(`is_transitive/autorel`(A)(R)) then
  reason := [convert(procname,string),"R is not transitive",R];
  return false;
 fi;

 return true;
end;

`is_equal/equiv`     := eval(`is_equal/autorel`);
`is_leq/equiv`       := eval(`is_leq/autorel`);
`is_separated/equiv` := eval(`is_separated/autorel`);

`is_total/equiv`     := (A::set) -> proc(R)
 local U,a,b;

 U := {seq(seq([a,b],b in A),a in A)} minus R;

 return evalb(U = {});
end;

`block_partition/equiv` := (A::set) -> proc(R)
 local pi,X,x,B;

 pi := NULL;
 X := A;

 while X <> {} do
  x := X[1];
  B := map(ab -> ab[2],select(ab -> ab[1] = x,R));
  pi := pi,B;
  X := X minus B;
 od:

 pi := {pi};

 return pi;
end:

`block_equiv/partition` := (A::set) -> proc(pi)
 local B,x,y;
 {seq(seq(seq([x,y],y in B),x in B),B in pi)};
end;

`equiv_is_equiv` := (A::set) -> (R) -> (a,b) ->
  member([a,b],R);

`random_element/equiv` := (A::set) -> proc()
 `block_equiv/partition`(A)(`random_element/partitions`(A)());
end;

`list_elements/equiv` := proc(A::set)
 map(`block_equiv/partition`(A),
     `list_elements/partitions`(A));
end:

`count_elements/equiv` := (A::set) -> combinat[bell](nops(A));


text:
A:B;
C:D;
|A B C D|
|^A _A ^^^B ^^_A|
[
  |A' B, A'' ^B,',|
  [
    |A B C D|
  ]
]

max_stars_repo_name: dwysocki/music-programming-language

max_issues_repo_name: dwysocki/music-programming-language

max_forks_repo_name: dwysocki/music-programming-language


text:
program assigninteger;	{sample21}
var n:integer;begin n:=1 end.

max_issues_repo_name: shouth/LanguageProcessing


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

`is_element/standard_stasheff_star` := (n::posint) -> proc(t)
 local JJ,TT,J,T,u,i,j,k;
 global reason;

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

 JJ := {seq(seq({seq(k,k=i..j)},j=i..n),i=1..n)};

 if map(op,{indices(t)}) <> JJ then
  reason := [convert(procname,string),"t is not indexed by the set of intervals",t,JJ];
  return false;
 fi;

 TT := NULL;

 for J in JJ do
  u := t[J];

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

  if (nops(J) = 1 or nops(J) = n) and u <> 1 then
   return false;
  fi;

  if u > 0 then TT := TT,J; fi;
 od;
 TT := {TT};

 if not(`is_element/trees`({seq(i,i=1..n)})(TT)) then
  reason := [convert(procname,string),"TT is not a tree",TT,reason];
  return false;
 fi;

 return true;
end;

`is_equal/standard_stasheff_star` := (n::posint) -> proc(t,u)
 local JJ,J,i,j,k;

 JJ := {seq(seq({seq(k,k=i..j)},j=i..n),i=1..n)};

 for J in JJ do
  if t[J] <> u[J] then
   return false;
  fi;
 od;
 return true;
end:

`is_leq/standard_stasheff_star` := NULL;

`random_element/standard_stasheff_star` := (n::posint) -> proc()
 local JJ,TT,T,t,i,j,k,d;

 d := 12;

 JJ := {seq(seq({seq(k,k=i..j)},j=i..n),i=1..n)};
 TT := `random_element/standard_stasheff_trees`(n)();
 t := table();

 for T in JJ do t[T] := 0; od;
 t[{seq(i,i=1..n)}] := 1;
 for i from 1 to n do t[{i}] := 1; od;

 for T in TT do
  if nops(T) > 1 and nops(T) < n then
   t[T] := rand(0..d)()/d;
  fi;
 od;

 return eval(t);
end;

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


text:
with(DifferentialThomas):
with(Tools):
Ranking([t], [[x15, x13, x7, x9, x2, x8, x14, x5, x4, x12, x6, x1, x10, x11, x3], [y4, y6, y2, y1, y3, y5] , [u]]):
sys := [
diff(x15(t), t) - (-1/2500*x15(t) + 1/2000000*x7(t)),
diff(x13(t), t) - (e2a*x14(t) + 1/2*x10(t)*x6(t) - x2(t)*x13(t) - 1/50000*x13(t)),
diff(x7(t), t) - (5*i1*x6(t) - 1/2*x11(t)*x7(t)),
diff(x9(t), t) - (-1/2500*x9(t) - 1/2000000*x7(t)),
diff(x2(t), t) - (-1/5*x10(t)*x2(t) - x2(t)*x13(t) - x2(t)*k_deg - x2(t)*k2*u(t)*x8(t) - x2(t)*k3 + x1(t)*u(t)*k1 + x4(t)*t1 + x5(t)*t2),
diff(x8(t), t) - (1/2*x9(t) - c5*x8(t)),
diff(x14(t), t) - (-5*e2a*x14(t) + 1/2*x11(t)*x7(t)),
diff(x5(t), t) - (x2(t)*x13(t) - x5(t)*t2),
diff(x4(t), t) - (1/5*x10(t)*x2(t) - x4(t)*t1),
diff(x12(t), t) - (1/2000000*x7(t) - c3a*x12(t)),
diff(x6(t), t) - (-1/2*x10(t)*x6(t) + 1/50000*x13(t) - i1*x6(t) - x5(t)*t2),
diff(x1(t), t) - (-x1(t)*k_deg - x1(t)*u(t)*k1 + k_prod),
diff(x10(t), t) - (-1/5*x10(t)*x2(t) - x10(t)*i1a - 1/2*x10(t)*x6(t) - 1/10000*x10(t) + 1/2000*x11(t) + x12(t)*c4a),
diff(x11(t), t) - (5*x10(t)*i1a - 1/2*x11(t)*x7(t) - 1/400*x11(t)),
diff(x3(t), t) - (x2(t)*k2*u(t)*x8(t) + x2(t)*k3 - x3(t)*k_deg),
y4(t) - (x2(t) + x3(t) + x1(t)),
y6(t) - (x12(t)),
y2(t) - (x10(t) + x13(t)),
y1(t) - (x7(t)),
y3(t) - (x9(t)),
y5(t) - (x2(t))
];
res := CodeTools[CPUTime](ThomasDecomposition(sys));


text:
repeat process ACV_antiBridges {
// Wake up one idle anti bridge after beeing teleported
/execute @e[type=area_effect_cloud,name=ACV_AntiBridge_Stopping] ~ ~ ~ entitydata @e[type=area_effect_cloud,name=ACV_AntiBridge_Idle,dy=0,c=1] {CustomName:"ACV_AntiBridge"}
/scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge] add ACV_TpProof

/testfor @e[type=area_effect_cloud,name=ACV_AntiBridge]
invert: stop

/kill @e[type=area_effect_cloud,name=ACV_AntiBridge_Stopping]

/scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=ACV_InCarpet] remove ACV_InCarpet
/execute @e[type=area_effect_cloud,name=ACV_AntiBridge] ~ ~ ~ detect ~ ~ ~ carpet 3 scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge,r=0,c=1] add ACV_InCarpet

/entitydata @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=!ACV_InCarpet] {CustomName:"ACV_AntiBridge_Stopping"}
/scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge_Stopping] remove ACV_TpProof

// To avoid race conditions, if there are multiple antibridges in the same block one needs to be processed (setblock) after the others, so that it removes the bridge correctly
/scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge] add ACV_Multiple
/execute @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=ACV_Multiple] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge,dy=0,c=-1] remove ACV_Multiple

/execute @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=ACV_Multiple] ~ ~ ~ kill @e[type=area_effect_cloud,name=ACV_BridgeCross,dy=0,c=1]

/scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=ACV_KeepCarpet] remove ACV_KeepCarpet
/scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=ACV_KeepPistonExtension] remove ACV_KeepPistonExtension
/execute @e[type=area_effect_cloud,name=ACV_BridgeCross] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge,dy=0] add ACV_KeepCarpet
/execute @e[type=area_effect_cloud,name=ACV_BridgeBelow] ~ ~-1 ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_AntiBridge,dy=0] add ACV_KeepPistonExtension

/execute @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=!ACV_Multiple] ~ ~ ~ kill @e[type=area_effect_cloud,name=ACV_BridgeCross,dy=0,c=1]

/execute @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=!ACV_KeepCarpet] ~ ~ ~ kill @e[type=area_effect_cloud,name=ACV_BridgeBelow,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=!ACV_KeepCarpet] ~ ~1 ~ kill @e[type=area_effect_cloud,name=ACV_BridgeBelow,dy=0]

///setblock ~ ~ ~ piston_extension Not: /fill ~ ~ ~ ~ ~ ~ piston_extension 0 replace carpet 3 due to https://bugs.mojang.com/browse/MC-61864
/execute @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=!ACV_KeepCarpet] ~ ~ ~ execute @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=ACV_KeepPistonExtension,r=0,c=1] ~ ~ ~ setblock ~ ~ ~ piston_extension
/execute @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=!ACV_KeepCarpet] ~ ~ ~ fill ~ ~ ~ ~ ~ ~ air 0 replace carpet 3
/execute @e[type=area_effect_cloud,name=ACV_AntiBridge,tag=!ACV_KeepCarpet] ~ ~ ~ fill ~ ~-1 ~ ~ ~-1 ~ air 0 replace piston_extension

/tp @e[type=area_effect_cloud,name=ACV_AntiBridge,rym=45,ry=135] ~-1 ~ ~
/tp @e[type=area_effect_cloud,name=ACV_AntiBridge,rym=135,ry=-135] ~ ~ ~-1
/tp @e[type=area_effect_cloud,name=ACV_AntiBridge,rym=-135,ry=-45] ~1 ~ ~
/tp @e[type=area_effect_cloud,name=ACV_AntiBridge,rym=-45,ry=45] ~ ~ ~1

/execute @e[type=area_effect_cloud,name=ACV_AntiBridge] ~ ~ ~ kill @e[type=area_effect_cloud,name=ACV_LightBridgeEnd,dy=0]
}


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

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

n := rs -> 1/(2*rs):

f := (rs, zeta) -> \
 (params_a_a1 + params_a_a2*n(rs) + params_a_a3*n(rs)^2) * n(rs)^params_a_alpha:

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
$include <binomsums.mpl>

$ifdef BUILD_MLA
savelib(BinomSums);
$endif

max_stars_repo_name: lairez/binomsum

max_issues_repo_name: lairez/binomsum

max_forks_repo_name: lairez/binomsum


text:
with(DifferentialAlgebra):
ring_diff := DifferentialRing(blocks = [[N, S, M, P, E], [y1, y3, y2] ], derivations = [t]):
sys := [
diff(N(t), t) - (-delta_NE*N(t)*P(t) - mu_N*N(t)),
diff(S(t), t) - (-mu_LE*S(t)*E(t) + delta_EL*S(t) - S(t)^2*mu_LL - S(t)*delta_LM),
diff(M(t), t) - (S(t)*delta_LM - M(t)*mu_M),
diff(P(t), t) - (rho_P*P(t)^2 - S(t)*P(t)*mu_PL - E(t)*mu_PE*P(t) - P(t)*mu_P),
diff(E(t), t) - (-mu_EE*E(t)^2 + delta_NE*N(t)*P(t) - delta_EL*E(t) + E(t)*P(t)*rho_E),
y1(t) - (N(t)),
y3(t) - (M(t)),
y2(t) - (S(t) + E(t))
];
res := CodeTools[CPUTime](RosenfeldGroebner(sys, ring_diff, singsol=none));


text:
# E(N) is the set of polynomials of the form 0 or +/- epsilon^d
# with 0 <= d < N.

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

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

 if not(`is_element/P`(N)(u)) then
  reason := [convert(procname,string)];
  return false;
 fi;

 if (u = 0) then return true; fi;

 d := degree(u,epsilon);
 if u <> epsilon^d and u <> -epsilon^d then
  reason := [convert(procname,string),"u <> +/- epsilon^d",u,d];
  return false;
 fi;

 return true;
end;

######################################################################
# `is_leq/E`(N)(u,v) is true if v-u is zero or is c epsilon^d + higher
# with c > 0.

`is_leq/E` :=  (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:

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

`is_preceq/E` := (N::posint) -> proc(u,v)
 if (u = v) then return true; fi;
 if (u = 0) then return true; fi;
 if (v = 0) then return false; fi;
 return evalb(ldegree(u,epsilon) > ldegree(v,epsilon));
end:

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

`list_elements/E` := proc(N::posint)
 local i;
 return [seq(-epsilon^i,i=0..N-1),0,seq(epsilon^(N-1-i),i=0..N-1)];
end:

`random_element/E` := (N::posint) -> proc()
 local L;
 L := `list_elements/E`(N);
 L[rand(1..nops(L))()];
end:

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

`count_elements/E` := (N::posint) -> 2*N+1;

`abs/E` := (N::posint) -> proc(u)
 if (u = 0) then return 0; fi;
 return epsilon^degree(u,epsilon);
end:

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

`lambda/P/E` := (N::posint) -> proc(u)
 local c,d;

 if u = 0 then
  return 0;
 else
  d := ldegree(u,epsilon);
  c := coeff(u,epsilon,d);
  return signum(c) * epsilon^d;
 fi;
end:

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

`delta/P` := (N::posint) -> (u,v) -> `lambda/P/E`(N)(v-u);


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

cases := [
    [
        [[
            diff(x1(t), t) - (a * x1(t) - b * x1(t) * x2(t)),
            diff(x2(t), t) - (-c * x2(t) + d * x1(t) * x2(t) + u(t)),
            y(t) = x1(t)
        ],
        [x1, x2], [y], [u], [a, b, c, d]],
        [y(t) * diff(diff(y(t), t), t) - diff(y(t), t)^2 - y(t)^2 * diff(y(t), t) * d + y(t) * diff(y(t), t) * c + a * d * y(t)^3 + u(t) * y(t)^2 * b - a * c * y(t)^2]
    ],
    [
        [[
            diff(xa(t), t) + k1 * xa(t),
            diff(xb(t), t) - k1 * xa(t) + k2 * xb(t),
            diff(xc(t), t) - k2 * xb(t),
            diff(eA(t), t),
            diff(eC(t), t),
            y1(t) - xc(t),
            y2(t) - eA(t) * xa(t) - eB * xb(t) - eC(t) * xc(t),
            y3(t) - eA(t),
            y4(t) - eC(t)
        ],
        [xa, xb, xc, eA, eC], [y2, y1, y3, y4], [], [k1, k2, eB]],
        [
            k1 * k2 * (y2(t) - y1(t) * y4(t)) - eB * k1 * diff(y1(t), t) - k2 * y3(t) * diff(y1(t), t) - y3(t) * diff(diff(y1(t), t), t),
            diff(diff(diff(y1(t), t), t), t) + (k1 + k2) * diff(diff(y1(t), t), t) + k1 * k2 * diff(y1(t), t),
            diff(y3(t), t),
            diff(y4(t), t)
        ]
    ]
]:

num_passed := 0:
num_failed := 0:

for case in cases do
    input := case[1]:
    print(input);
    correct := case[2]:
    result := GetIOEquations(op(input), 0):
    passed := true:
    if nops(result) <> nops(correct) then
        passed := false:
    else
        for i from 1 to nops(result) do
            if simplify(result[i] - correct[i]) <> 0 then
                passed := false:
            end if:
        end do:
    end if:
    if passed = true then
        printf("PASSED");
        num_passed := num_passed + 1:
    else
        printf("FAILED");
        print(result);
        print(correct);
        num_failed := num_failed + 1:
    end if:
end do:

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

max_stars_repo_name: iliailmer/AllIdentifiableFunctions

max_issues_repo_name: iliailmer/AllIdentifiableFunctions

max_forks_repo_name: iliailmer/AllIdentifiableFunctions


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

$define lda_c_vbh_params
$include "lda_c_hl.mpl"

malpha := (4*Pi/3)^(1/6):

(* 4.28e-3/2, where the 2 comes from the covertion from Ryd. to Hartree *)
a1     := Pi/(16*(3*Pi^2)^(4/3)):
a2     := -0.262:
a3     := -7/(9*2^(5/3)):

t1 := (rs, z, xt) ->
  + xt^2*exp(a2*xt/(malpha*sqrt(rs)))
  * sqrt(2)/sqrt((1 + z)^(5/3) + (1 - z)^(5/3)):
t2 := (z, xs0, xs1) ->
  + a3*(xs0^2*(1 + z)^(4/3) + xs1^2*(1 - z)^(4/3)):

f := (rs, z, xt, xs0, xs1) ->
  + f_hl(rs, z)
  + a1/(malpha^2*rs)*(t1(rs, z, xt) + t2(z, xs0, xs1)):


text:
######################################################################
#
# General operad framework.

`check_semi_gamma_args` := proc(J,U,V,is_el_CC)
 local i,k,jj;
 global reason;

 if not type(J,list(nonnegint)) then
  reason := [convert(procname,string),"J is not a list of natural numbers",J];
  return false;
 fi;

 k := nops(J);
 jj := `+`(op(JJ));

 if not is_el_CC(k)(U) then
  reason := [convert(procname,string),"U is not in C(k)",eval(U),k];
  return false;
 fi;

 if not (type(V,list) and nops(V) = k) then
  reason := [convert(procname,string),"V is not a list of length k",eval(V),k];
  return false;
 fi;

 for i from 1 to k do
  if not is_el_CC(J[i])(V[i]) then
   reason := [convert(procname,string),"V[i] is not in C(j[i])",i,j[i],eval(V[i])];
   return false;
  fi;
 od;

 return true;
end;

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

`check_semi_gamma_axiom` := proc(L,U,V,W,is_el_CC,is_equal_CC,gamma_CC)
 local J,L_,k,jj,ii,H,k0,j0,W0,W_,UV,VW0,VW,UVW0,UVW1;
 global reason;

 if not type(L,list(list(nonnegint))) then
  reason := [convert(procname,string),"L is not a list of lists of natural numbers",L];
  return false;
 fi;

 J := map(nops,L);
 L_ := map(op,L);
 k := nops(J);
 jj := `+`(op(J));
 ii := `+`(op(L_));
 H := [seq(`+`(op(L[k0])),k0=1..k)];

 if not is_el_CC(k)(U) then
  reason := [convert(procname,string),"U is not in O(k)",eval(U),k];
  return false;
 fi;

 if not (type(V,list) and nops(V) = k) then
  reason := [convert(procname,string),"V is not a list of length k",eval(V),k];
  return false;
 fi;

 for k0 from 1 to k do
  if not is_el_CC(J[k0])(V[k0]) then
   reason := [convert(procname,string),"V[k0] is not in O(J[k0])",k0,J[k0],eval(V[k0])];
   return false;
  fi;
 od;

 if not (type(W,list) and nops(W) = k) then
  reason := [convert(procname,string),"W is not a list of length k",eval(W),k];
  return false;
 fi;

 for k0 from 1 to k do
  W0 := W[k0];
  if not (type(W0,list) and nops(W0) = J[k0]) then
   reason := [convert(procname,string),"W[k0] is not a list of length J[k0]",k0,J[k0],W0];
   return false;
  fi;

  for j0 from 1 to J[k0] do
   if not is_el_CC(L[k0][j0])(W[k0][j0]) then
    reason := [convert(procname,string),"W[k0][j0] is not in O(L[k0][j0])",k0,j0,L[k0][j0],eval(W[k0][j0])];
    return false;
   fi;
  od:
 od:

 UV := gamma_CC(J)(U,V);
 if not is_el_CC(jj)(UV) then
  reason := [convert(procname,string),"UV is not in O(jj)",jj,J,eval(U),eval(V),eval(UV)];
  return false;
 fi;

 W_ := map(op,W);

 if not `check_semi_gamma_args`(L_,UV,W_,is_el_CC) then
  reason := [convert(procname,string),"invalid arguments to gamma 1",reason];
  return false;
 fi;

 UVW0 := gamma_CC(L_)(UV,W_);
 if not is_el_CC(ii)(UVW0) then
  reason := [convert(procname,string),"UVW0 is not in O(ii)",ii,eval(UVW0)];
  return false;
 fi;

 VW := [];
 for k0 from 1 to k do
  if not `check_semi_gamma_args`(L[k0],V[k0],W[k0],is_el_CC) then
   reason := [convert(procname,string),"invalid arguments to gamma 2",reason];
   return false;
  fi;
  VW0 := gamma_CC(L[k0])(V[k0],W[k0]);
  if not is_el_CC(H[k0])(VW0) then
   reason := [convert(procname,string),"VW[k0] is not in O(H[k0])",k0,H[k0],eval(VW0)];
   return false;
  fi;
  VW := [op(VW),VW0];
 od;

 if not `check_semi_gamma_args`(H,U,VW,is_el_CC) then
  reason := [convert(procname,string),"invalid arguments to gamma 3",reason];
  return false;
 fi;

 UVW1 := gamma_CC(H)(U,VW);
 if not is_el_CC(ii)(UVW1) then
  reason := [convert(procname,string),"UVW1 is not in O(ii)",ii,eval(UVW1)];
  return false;
 fi;

 if not is_equal_CC(ii)(UVW0,UVW1) then
  reason := [convert(procname,string),"UVW0 != UVW1",eval(UVW0),eval(UVW1)];
  return false;
 fi;

 return true;
end;


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

`eta/trees` := (A::set) -> `if`(nops(A)=1,{A},FAIL);

`gamma/trees` := (A::set,B::set) -> (p) -> proc(UU,TT)
 local F,U,V,VV,b;

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

 VV := NULL;
 for U in UU do
  V := `union`(seq(F[b],b in U));
  if V <> {} then VV := VV,V; fi;
 od;

 for b in B do
  VV := VV,op(TT[b]);
 od;

 VV := sort(map(sort,{VV}));

 return VV;
end;


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

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

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

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

(* standard PBE *)
$ifdef gga_x_pbe_params
params_a_kappa := 0.8040:
params_a_mu    := 0.2195149727645171:
$endif

(* PBE_SOL *)
$ifdef gga_x_pbe_sol_params
params_a_kappa := 0.8040:
params_a_mu    := MU_GE:
$endif

$ifdef gga_x_pbe_tca_params
params_a_kappa := 1.227:
params_a_mu    := 0.2195149727645171:
$endif

pbe_f0 := s -> 1 + params_a_kappa*(1 - params_a_kappa/(params_a_kappa + params_a_mu*s^2)):
pbe_f  := x -> pbe_f0(X2S*x):

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


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

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

(* type: lda_exc *)

$define lda_c_pw_params
$define lda_c_pw_modified_params
$include "lda_c_pw.mpl"

alpha := (4/(9*Pi))^(1/3):

# Eq. (46) of Phys. Rev. A 73, 032506 (2006)
# I am not sure if this is the on-top pair correlation functional
# intended by the authors
g0_C    :=  0.08193:
g0_D    := -0.01277:
g0_E    :=  0.001859:
g0_dd   :=  0.7524:
g0_aHD  := -alpha*(Pi^2 + 6*log(2) - 3)/(5 * Pi):
g0_B    := -2*g0_aHD - g0_dd:
pmgb_g0 := rs ->
        (1 - g0_B*rs + g0_C*rs^2 + g0_D*rs^3 + g0_E*rs^4)*exp(-g0_dd*rs)/2:

pmgb_gc0 := rs -> pmgb_g0(rs) - 1/2:

# Eq. (14)
pmgb_phi := (n, z) ->
         (opz_pow_n(z,n/3) + opz_pow_n(-z,n/3))/2:

# Eq. (22)
pmgb_Q_a := 5.84605:
pmgb_Q_c := 3.91744:
pmgb_Q_d := 3.44851:
pmgb_Q_b := pmgb_Q_d - 3*Pi*alpha/(4*log(2) - 4):

pmgb_Q := x ->
       ((2*log(2) - 2)/Pi^2) * log(
        (1 + pmgb_Q_a*x + pmgb_Q_b*x^2 + pmgb_Q_c*x^3)/
        (1 + pmgb_Q_a*x + pmgb_Q_d*x^2)
       ):

# Eq. (31)
pmgb_b0 := (rs) -> 0.784949*rs:

# Eq. (32)
pmgb_gpp1 := (rs) -> (2^(5/3)/(5*alpha^2*rs^2)) *
          (1 - 0.02267*rs)/(1 + 0.4319*rs + 0.04*rs^2):

# Eq. (33)
pmgb_D2 := (rs) -> (-0.388*rs + 0.676*rs^2)*exp(-0.547*rs)/rs^2:

# Eq. (34)
pmgb_D3 := (rs) -> (-4.95*rs + rs^2)*exp(-0.31*rs)/rs^3:

# Eq. (28)
pmgb_cc4 := (rs, z) ->
       + (opz_pow_n( z,2)/4) * pmgb_gpp1(rs*(2/(1 + z))^(1/3))
       + (opz_pow_n(-z,2)/4) * pmgb_gpp1(rs*(2/(1 - z))^(1/3))
       + (1 - z^2) * pmgb_D2(rs)
       - pmgb_phi(8, z)/(5*alpha^2*rs^2):

# Eq. (29)
pmgb_cc5 := (rs, z) ->
       + (opz_pow_n( z,2)/4) * pmgb_gpp1(rs*(2/(1 + z))^(1/3))
       + (opz_pow_n(-z,2)/4) * pmgb_gpp1(rs*(2/(1 - z))^(1/3))
       + (1 - z^2) * pmgb_D3(rs):

# Eq. (30)
pmgb_C2 := (rs, z) ->
        -3*(1 - z^2)*pmgb_gc0(rs)/(8*rs^3):

pmgb_C3 := (rs, z) ->
        -(1 - z^2)*pmgb_g0(rs)/(sqrt(2*Pi)*rs^3):

pmgb_C4 := (rs, z) ->
        -9*pmgb_cc4(rs, z)/(64*rs^3):

pmgb_C5 := (rs, z) ->
        -9*pmgb_cc5(rs, z)/(40*sqrt(2*Pi)*rs^3):

pmgb_a1 := (rs, z) ->
        4*pmgb_b0(rs)^6*pmgb_C3(rs, z) + pmgb_b0(rs)^8*pmgb_C5(rs, z):

pmgb_a2 := (rs, z) ->
        4*pmgb_b0(rs)^6*pmgb_C2(rs, z) + pmgb_b0(rs)^8*pmgb_C4(rs, z)
        + 6*pmgb_b0(rs)^4*f_pw(rs, z):

pmgb_a3 := (rs, z) ->
        pmgb_b0(rs)^8*pmgb_C3(rs, z):

pmgb_a4 := (rs, z) ->
        pmgb_b0(rs)^8*pmgb_C2(rs, z) + 4*pmgb_b0(rs)^6*f_pw(rs, z):

pmgb_a5 := (rs, z) ->
        pmgb_b0(rs)^8*f_pw(rs, z):

# Eq. (26)
pmgb_ec_LR := (rs, z) -> (
           + pmgb_phi(2, z)^3*pmgb_Q(p_a_hyb_omega_0_*sqrt(rs)/pmgb_phi(2, z))
           + pmgb_a1(rs, z)*p_a_hyb_omega_0_^3
           + pmgb_a2(rs, z)*p_a_hyb_omega_0_^4
           + pmgb_a3(rs, z)*p_a_hyb_omega_0_^5
           + pmgb_a4(rs, z)*p_a_hyb_omega_0_^6
           + pmgb_a5(rs, z)*p_a_hyb_omega_0_^8
           )/(1 + pmgb_b0(rs)^2*p_a_hyb_omega_0_^2)^4:

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


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

cases := [
    [
        [a * x + a * y, [x, y], [a]], 1
    ],
    [
        [(a + b) * (x + y) - (b + c) * (y + z) + 3 * (a + c) * (x + z), [x, y, z], [a, b, c]], 3
    ]
];

num_passed := 0:
num_failed := 0:

for case in cases do
    input := case[1]:
    correct := case[2]:
    passed := true:
    result := DecomposePolynomial(op(input), 0):
    if nops(result[1]) <> correct then
        passed := false:
    end if:
    orig_poly := 0:
    for i from 1 to nops(result[1]) do
        orig_poly := orig_poly + result[1][i] * result[2][i]:
    end do:
    if simplify(orig_poly - input[1]) <> 0 then
        passed := false:
    end if:
    if passed = true then
        printf("PASSED");
        num_passed := num_passed + 1:
    else
        printf("FAILED");
        num_failed := num_failed + 1:
        print(result);
    end if:
end do:

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


text:
flavor 1
srclang 3
id 65535
numfuncs 2
import "HelloWorld.mplt"
import "/home/bravewtz/Desktop/openarkcompiler/libjava-core/libjava-core.mplt"
entryfunc &LHelloWorld_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V
fileinfo {
  @INFO_filename "HelloWorld.jar"}
srcfileinfo {
  1 "HelloWorld.java"}
type $__class_meta__ <struct {
  @shadow <* void> final,
  @monitor i32 final,
  @classloader u16 final,
  @objsize u16 final,
  @itab <* void> final,
  @vtab <* void> final,
  @gctib <* void> final,
  @classinforo <* void> final,
  @clinitbridge <* void> final}>
type $__class_meta_ro__ <struct {
  @classname <* void> final,
  @ifields <* void> final,
  @methods <* void> final,
  @superclass_or_componentclass <* void> final,
  @numoffields u16 final,
  @numofmethods u16 final,
  @flag u16 final,
  @numofsuperclasses u16 final,
  @padding u32 final,
  @mod i32 final,
  @annotation i32 final,
  @clinitAddr i32 final}>
type $__method_info__ <struct {
  @method_in_vtab_index i64 final,
  @declaringclass i64 final,
  @addr i64 final,
  @mod i32 final,
  @methodname i32 final,
  @signaturename i32 final,
  @annotationvalue i32 final,
  @flag u16 final,
  @argsize u16 final,
  @padding u32 final}>
type $__method_info_compact__ <struct {
  @method_in_vtab_index i32 final,
  @addr i32 final,
  @lebPadding0 u8 final}>
type $__field_info__ <struct {
  @offset u64 final,
  @mod u32 final,
  @flag u16 final,
  @index u16 final,
  @typeName i64 final,
  @fieldname u32 final,
  @annotation u32 final,
  @declaringclass <* void> final}>
type $__field_info_compact__ <struct {
  @offset u32 final,
  @lebPadding0 u8 final}>
type $__superclass_meta__ <struct {
  @superclassinfo <* void> final}>
type $MUIDFuncDefTabEntry <struct {
  @funcUnifiedAddr u64}>
type $MUIDFuncInfTabEntry <struct {
  @funcSize u32,
  @funcName u32}>
type $MUIDFuncDefMuidTabEntry <struct {
  @muidLow u64,
  @muidHigh u64}>
type $MUIDDataDefTabEntry <struct {
  @dataUnifiedAddr u64}>
type $MUIDDataDefMuidTabEntry <struct {
  @muidLow u64,
  @muidHigh u64}>
type $MUIDUnifiedUndefTabEntry <struct {
  @globalAddress u64}>
type $MUIDUnifiedUndefMuidTabEntry <struct {
  @muidLow u64,
  @muidHigh u64}>
type $MUIDRangeTabEntry <struct {
  @tabBegin <* void>,
  @tabEnd <* void>}>
javaclass $LHelloWorld_3B <$LHelloWorld_3B> public
var $_C_STR_d9bd4a0c2c2117158ed933ab7468a461 fstatic <[4] u64> readonly = [0, 0x1900000000, 0x6c6c6548c63cb61d, 0x21646c726f57206f]
var $Ljava_2Flang_2FSystem_3B_7Cout extern <* <$Ljava_2Fio_2FPrintStream_3B>> final public static
var $__cinf_Ljava_2Flang_2FString_3B <$__class_meta__>
func &MCC_GetOrInsertLiteral () <* <$Ljava_2Flang_2FString_3B>>
var $__vtb_LHelloWorld_3B fstatic <[11] <* void>> = [16, 24, 32, 36, 8, 4, 48, 12, 20, 44, 28]
var $__cinf_LHelloWorld_3B <$__class_meta__> public
var $__methods_info__LHelloWorld_3B fstatic <[2] <$__method_info__>> public = [[1= 0xfff6, 2= addrof ptr $__cinf_LHelloWorld_3B, 3= addroffunc ptr &LHelloWorld_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V, 4= 9, 5= 56, 6= 76, 7= 168, 8= 0xa201, 9= 1, 10= 0], [1= 0xfff6, 2= addrof ptr $__cinf_LHelloWorld_3B, 3= addroffunc ptr &LHelloWorld_3B_7C_3Cinit_3E_7C_28_29V, 4= 0x10001, 5= 184, 6= 212, 7= 168, 8= 0xad81, 9= 1, 10= 0]]
var $__cinf_Ljava_2Flang_2FObject_3B extern <$__class_meta__> public
var $__superclasses__LHelloWorld_3B fstatic <[1] <$__superclass_meta__>> public = [[1= 0x4000000000000000]]
var $__classinforo__LHelloWorld_3B fstatic <$__class_meta_ro__> public = [1= 4, 2= 0, 3= addrof ptr $__methods_info__LHelloWorld_3B, 4= addrof ptr $__superclasses__LHelloWorld_3B, 5= 0, 6= 2, 7= 0, 8= 1, 9= 0, 10= 33, 11= 168, 12= 0]
var $MCC_GCTIB__LHelloWorld_3B fstatic <* void> public
var $classStateInitialized u64
var $__cinf_LHelloWorld_3B <$__class_meta__> public = [1= 0x26ea31dfa0409ef, 2= 0, 3= 0xffff, 4= 0, 5= 0, 6= addrof ptr $__vtb_LHelloWorld_3B, 7= addrof ptr $MCC_GCTIB__LHelloWorld_3B, 8= addrof ptr $__classinforo__LHelloWorld_3B, 9= addrof ptr $classStateInitialized]
var $__muid_classmetadata_bucket$$HelloWorld_jar <[1] <* void>> public = [addrof ptr $__cinf_LHelloWorld_3B]
func &MCC_Reflect_ThrowCastException nosideeffect () void
func &MCC_Reflect_Check_Casting_NoArray nosideeffect () void
func &MCC_Reflect_Check_Casting_Array nosideeffect () void
func &MCC_CheckThrowPendingException nosideeffect () void
func &MCC_PreNativeCall (var %caller ref) <* void>
func &MCC_PostNativeCall (var %env <* void>) void
func &MCC_DecodeReference nosideeffect (var %obj ref) ref
func &MCC_CallFastNative (var %func <* void>) <* void>
func &MCC_CallSlowNative0 (var %func <* void>) <* void>
func &MCC_CallSlowNative1 (var %func <* void>) <* void>
func &MCC_CallSlowNative2 (var %func <* void>) <* void>
func &MCC_CallSlowNative3 (var %func <* void>) <* void>
func &MCC_CallSlowNative4 (var %func <* void>) <* void>
func &MCC_CallSlowNative5 (var %func <* void>) <* void>
func &MCC_CallSlowNative6 (var %func <* void>) <* void>
func &MCC_CallSlowNative7 (var %func <* void>) <* void>
func &MCC_CallSlowNative8 (var %func <* void>) <* void>
func &MCC_CallFastNativeExt (var %func <* void>) <* void>
func &MCC_CallSlowNativeExt (var %func <* void>) <* void>
func &MCC_SetReliableUnwindContext nosideeffect () void
var $__reg_jni_func_tab$$HelloWorld_jar <[0] <* void>>
var $__cinf_Ljava_2Flang_2FSystem_3B extern <$__class_meta__>
func &MCC_getFuncPtrFromItabSecondHash64 nosideeffect () ptr
var $__muid_func_def_tab$$HelloWorld_jar fstatic <[2] <$MUIDFuncDefTabEntry>> = [[1= addroffunc ptr &LHelloWorld_3B_7C_3Cinit_3E_7C_28_29V], [1= addroffunc ptr &LHelloWorld_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V]]
var $__muid_func_inf_tab$$HelloWorld_jar fstatic <[2] <$MUIDFuncInfTabEntry>> = [[1= addroffunc ptr &LHelloWorld_3B_7C_3Cinit_3E_7C_28_29V, 2= addroffunc ptr &LHelloWorld_3B_7C_3Cinit_3E_7C_28_29V], [1= addroffunc ptr &LHelloWorld_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V, 2= addroffunc ptr &LHelloWorld_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V]]
var $__muid_func_def_muid_tab$$HelloWorld_jar fstatic <[2] <$MUIDFuncDefMuidTabEntry>> = [[1= 0x960a4f24512c54e, 2= -1012656116509287985], [1= 0x406c510221f3746, 2= -1731984239878388817]]
var $__muid_func_muid_idx_tab$$HelloWorld_jar fstatic <[2] u32> = [1, 0]
var $__muid_data_def_tab$$HelloWorld_jar fstatic <[1] <$MUIDDataDefTabEntry>> = [[1= addrof ptr $__cinf_LHelloWorld_3B]]
var $__muid_data_def_muid_tab$$HelloWorld_jar fstatic <[1] <$MUIDDataDefMuidTabEntry>> = [[1= 0x52921ad85a4ce56, 2= -3592158421343910012]]
var $__muid_func_undef_tab$$HelloWorld_jar fstatic <[12] <$MUIDUnifiedUndefTabEntry>> = [[1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0]]
var $__muid_func_undef_muid_tab$$HelloWorld_jar fstatic <[12] <$MUIDUnifiedUndefMuidTabEntry>> = [[1= 0x3e32352aee789835, 2= -3887705395317205813], [1= -4187412136968710015, 2= -3800091941095621250], [1= 0x6742c234127e0a27, 2= -3762262047879347071], [1= 0x7230554331c55d92, 2= -3676689525926909155], [1= -715372855679083712, 2= -2647497990906227723], [1= -7464356948810446352, 2= -2259485500590180091], [1= 0x783627f2afd1cbde, 2= -2046851302095768916], [1= -2701934576591406938, 2= -1693831364093527548], [1= 0x126777a7fe39e1fb, 2= -1314856249532362766], [1= 0xbf40578f3343f7a, 2= -1198421541845410999], [1= 0x7ca2bdf69a6c7c94, 2= -801329978528900548], [1= 0x477aafa4d7dd102b, 2= -442561182569419835]]
var $__muid_data_undef_tab$$HelloWorld_jar fstatic <[3] <$MUIDUnifiedUndefTabEntry>> = [[1= addrof ptr $__cinf_Ljava_2Flang_2FObject_3B], [1= addrof ptr $Ljava_2Flang_2FSystem_3B_7Cout], [1= addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B]]
var $__muid_data_undef_muid_tab$$HelloWorld_jar fstatic <[3] <$MUIDUnifiedUndefMuidTabEntry>> = [[1= -567417612161374449, 2= -3298852447504547670], [1= 0x191283ac418c4bb9, 2= -1676204161023949463], [1= -5921653145571052587, 2= -171150348656858163]]
var $__muid_range_tab$$HelloWorld_jar fstatic <[29] <$MUIDRangeTabEntry>> = [[1= 0x7ad679d34e1fa4f0, 2= -434896966956065287], [1= 0x2b9bf9e7d26c703b, 2= -3810538229946804202], [1= 2, 2= 2], [1= 3, 2= 3], [1= 4, 2= 4], [1= 5, 2= 5], [1= 6, 2= 6], [1= 7, 2= 7], [1= 8, 2= 8], [1= 9, 2= 9], [1= 0, 2= 0], [1= 11, 2= 11], [1= 12, 2= 12], [1= 13, 2= 13], [1= 14, 2= 14], [1= 15, 2= 15], [1= addrof ptr $__muid_func_def_tab$$HelloWorld_jar, 2= addrof ptr $__muid_func_def_tab$$HelloWorld_jar], [1= 0, 2= 0], [1= addrof ptr $__muid_func_inf_tab$$HelloWorld_jar, 2= addrof ptr $__muid_func_inf_tab$$HelloWorld_jar], [1= addrof ptr $__muid_func_undef_tab$$HelloWorld_jar, 2= addrof ptr $__muid_func_undef_tab$$HelloWorld_jar], [1= addrof ptr $__muid_data_def_tab$$HelloWorld_jar, 2= addrof ptr $__muid_data_def_tab$$HelloWorld_jar], [1= 0, 2= 0], [1= addrof ptr $__muid_data_undef_tab$$HelloWorld_jar, 2= addrof ptr $__muid_data_undef_tab$$HelloWorld_jar], [1= addrof ptr $__muid_func_def_muid_tab$$HelloWorld_jar, 2= addrof ptr $__muid_func_def_muid_tab$$HelloWorld_jar], [1= addrof ptr $__muid_func_undef_muid_tab$$HelloWorld_jar, 2= addrof ptr $__muid_func_undef_muid_tab$$HelloWorld_jar], [1= addrof ptr $__muid_data_def_muid_tab$$HelloWorld_jar, 2= addrof ptr $__muid_data_def_muid_tab$$HelloWorld_jar], [1= addrof ptr $__muid_data_undef_muid_tab$$HelloWorld_jar, 2= addrof ptr $__muid_data_undef_muid_tab$$HelloWorld_jar], [1= addrof ptr $__muid_func_muid_idx_tab$$HelloWorld_jar, 2= addrof ptr $__muid_func_muid_idx_tab$$HelloWorld_jar], [1= 0, 2= 0]]
var $__reflection_strtab$$HelloWorld_jar fstatic <[57] u8> = [0, 76, 72, 101, 108, 108, 111, 87, 111, 114, 108, 100, 59, 0, 109, 97, 105, 110, 0, 40, 91, 76, 106, 97, 118, 97, 47, 108, 97, 110, 103, 47, 83, 116, 114, 105, 110, 103, 59, 41, 86, 0, 48, 33, 48, 0, 60, 105, 110, 105, 116, 62, 0, 40, 41, 86, 0]
var $__compilerVersionNum$$HelloWorld_jar <[0] <* void>> = [1, 0]
func &LHelloWorld_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$LHelloWorld_3B>>) void {
  funcid 48153
  var %Reg1_R43694 <* <$LHelloWorld_3B>> localrefvar
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>> localrefvar
  var %__muid_symptr <* void>

  intrinsiccall MCCIncRef (dread ref %_this)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43694)
  dassign %Reg1_R43694 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  regassign ptr %1 (dread ref %Reg1_R57)
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43694))
  intrinsiccall MCCIncRef (dread ref %Reg1_R57)
  intrinsiccall MCCDecRef (regread ptr %1)
  #Call function:Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V
  dassign %__muid_symptr 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[12] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_func_undef_tab$$HelloWorld_jar, constval i64 9)))
  icallassigned (dread ptr %__muid_symptr, dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  intrinsiccall MPL_CLEANUP_LOCALREFVARS (dread ref %Reg1_R43694, dread ref %Reg1_R57)
  return ()
}
func &LHelloWorld_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg2_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void {
  funcid 48154
  var %Reg0_R562 <* <$Ljava_2Fio_2FPrintStream_3B>> localrefvar
  var %Reg1_R43 <* <$Ljava_2Flang_2FString_3B>> localrefvar
  var %L_STR_161334 <* <$Ljava_2Flang_2FString_3B>>

  #INSTIDX : 0||0000:  getstatic
  intrinsiccall MPL_CLINIT_CHECK (addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B)
  regassign ptr %1 (dread ref %Reg0_R562)
  #Read from: Ljava_2Flang_2FSystem_3B_7Cout
  dassign %Reg0_R562 0 (iread ref <* <* <$Ljava_2Fio_2FPrintStream_3B>>> 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[3] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_data_undef_tab$$HelloWorld_jar, constval i64 1))))
  intrinsiccall MCCIncRef (dread ref %Reg0_R562)
  intrinsiccall MCCDecRef (regread ptr %1)
  #INSTIDX : 3||0003:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_d9bd4a0c2c2117158ed933ab7468a461) { dassign %L_STR_161334 0 }
  intrinsiccall MCCIncRef (dread ptr %L_STR_161334)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43)
  dassign %Reg1_R43 0 (dread ptr %L_STR_161334)
  #INSTIDX : 5||0005:  invokevirtual
  icallassigned (
    iread u64 <* u64> 0 (add ptr (
      iread ptr <* <$__class_meta__>> 6 (iread ref <* <$Ljava_2Flang_2FObject_3B>> 1 (dread ref %Reg0_R562)),
      constval u32 312)),
    dread ref %Reg0_R562,
    dread ref %Reg1_R43) {}
  #INSTIDX : 8||0008:  return
  intrinsiccall MPL_CLEANUP_LOCALREFVARS (dread ref %Reg0_R562, dread ref %Reg1_R43)
  return ()
}

max_stars_repo_name: Bravewtz/OpenArkCompiler

max_issues_repo_name: Bravewtz/OpenArkCompiler

max_forks_repo_name: Bravewtz/OpenArkCompiler


text:
(*
    变量可以划分为：
    + 自由变量：在满足非零的前提下可以全部取0。
    + 确定性变量：取值已经确定，或者由其它变量确定。
    + 有约束变量：其约束来自于解的定义域，由解的性质可知，其值不由其它变量所确定。

    处理时主要处理约束条件，自由变量可以加单取0，确定性变量的值已知。在完全处理约束的前提下，无需验证是否满足约束。

    至于如何处理约束：
    + 和c有关的约束：为了能够始终满足约束，特解应当和c有关，可以采取变成等式约束进行求解的方式进行处理。(注意对于单独的c的约束，直接删除)
    + 和c无关的约束：可以简单求解不等式，求解后的结果，都是单边的不等式约束，可以看情况取解。

    对于求解后的结果，可以细分为：
    + 单变量约束，可以按照简单的方法取定特解。
    + 多变量约束，可以在取定其它变量的值之后，转化单变量约束进行处理。
*)

$ifndef _FETCH_
$define _FETCH_

$include "Condition.mpl"
$include "Utils.mpl"

# 根据方程和约束取特解
fetchSpecSol:=proc(sol::list,con::set,{nonzero::boolean:=false})
    local vars,vf,eq,ca,cc,res,ccv,sc;
    vars:=lhs~(sol);
    vf,eq,ca,cc:=filtCon({sol[],classifySolve(con)[]});
    if cc<>{} then
        # 对于和c有关的不等式约束，转化成等式约束后再进行求解
        cc:=transIeq~(cc);
        ccv:=indets(cc,name);
        ccv:=select(type,ccv,specindex(a));
        sc:=RealDomain:-solve(cc,[ccv[]],explicit);
        res:=map(rebuildAndFetch,sc,vars,vf,eq,ca,nonzero);
        res:=sortByComplexity(res);
        return res[1];
    else
        return fetchByCons(vars,vf,eq,ca,nonzero);
    end if;
end proc:

# 转化关于c的不等式约束之后，再对约束进行分类
rebuildAndFetch:=proc(s,vars,_vf,_eq,_ca,nonzero)
    local vf,eq,ca,cc;
    vf,eq,ca,cc:=filtCon({s[]} union findSolutionDomain(s));
    vf:=_vf minus indets(lhs~(eq),name);
    eq:=eq union _eq;
    ca:=ca union _ca;
    return fetchByCons(vars,vf,eq,ca,nonzero);
end proc:

# 求解关于a的约束之后，转化为单边约束，再逐解取特解
fetchByCons:=proc(vars,vf,eq,ca,nonzero)
    local sca,res;
    sca:=[RealDomain:-solve(ca)];
    res:=map[4](fetchBySingleCons,vars,vf,eq,sca,nonzero);
    if nonzero then
        # 在满足非零要求时，目的是给出所有可能的基
        # 而求解a的约束产生的每个解是不相交的
        # 所以需要合并各个解的基
        return `union`(res[]);
    else
        return sortByComplexity(res)[1];
    end if;
end proc:

# 对于和a有关的约束，采用逐步替换的方法取特解
fetchBySingleCons:=proc(vars,_vf,_eq,_ca,nonzero)
    local vf,eq,ca,sca,res,req,tmp,i,n,ra;
    vf:=_vf;# 自由变量
    eq:=_eq;# 不等式方程的解等式约束
    ca:=_ca;# 求解跟a有关的约束得到的解
    # 对于单变量约束，采用最邻近整数方式取特解
    # 逐步替换多变量约束中的变量，最终都取得特解
    while true do
        sca,ca:=selectremove(isSingleCon,ca);
        if sca={} then
            break;
        end if;
        sca:=fetchIeq~(sca);
        eq:=eq union sca;
        ca:=subs(sca[],ca);
    end do;
    # 由于对不等约束的求解可能存在等式约束，因此最后剩下的是等式约束
    # BUG  出错的问题在于上面得到的ca可能会有等式约束
    ca,ra:=selectremove(type,ca,equation);
    vf:={seq(x=0,x in vf)};
    eq:=eq union ca union vf;
    req:=vars;
    # 带入两次是因为带入一次可能不完全
    res:=eval(subs(eq[],req));
    res:=eval(subs(eq[],res));
    # 假设还有变量存在，就重新求解
    if indets(res,specindex(a))<>{} then
        res:=fetchSpecSol([seq(a[i]=res[i],i=1..numelems(vars))],{});
    end if;
    if nonzero then
        if andmap(type,res,0) then
            n:=numelems(vf);
            res:={};
            for i from 1 to n do
                tmp:=subsop([i,2]=1,vf);
                res:=res union {eval(subs(tmp[],eq[],req))};
            end do;
        else
            res:={res};
        end if;
    end if;
    return res;
end proc:

# 是否是单变量约束
isSingleCon:=proc(x)
    return not type(x,`=`) and numelems(indets(x,name))=1;
end proc:

# 筛选：自由变量，确定性约束，和c无关的不等式约束，和c有关的不等式约束
# 会删除之和c有关的约束
filtCon:=proc(con::set)
    local ef,eq,ca,cc,rst,vf;
    ef,rst:=selectremove(x->lhs(x)=rhs(x),con);
    eq,rst:=selectremove(type,rst,`=`);
    vf:=indets(ef,name) minus lhs~(eq) minus indets(rst,name);
    cc,ca:=selectremove(has,rst,c);
    cc:=select(has,cc,a);
    return vf,eq,ca,cc;
end proc:

# 将不等式约束转化为等式约束
transIeq:=proc(x)
    if op(0,x)=`<>` or op(0,x)=`<` then
        return lhs(x)=rhs(x)-1;
    elif op(0,x)=`>` then
        return lhs(x)=rhs(x)+1;
    else
        return lhs(x)=rhs(x);
    end if;
end proc:

# 对于 <> < <= 的约束条件取特解
# 取满足条件的最邻近整数
fetchIeq:=proc(x)
    local r;
    if type(x,`<>`) then
        r:=lhs(x)=floor(rhs(x))+1;
    elif type(lhs(x),numeric) then
        if type(x,`<`) then
            r:=rhs(x)=floor(lhs(x))+1;
        else
            r:=rhs(x)=ceil(lhs(x));
        end if;
    else
        if type(x,`<`) then
            r:=lhs(x)=ceil(rhs(x))-1;
        else
            r:=lhs(x)=floor(rhs(x));
        end if;
    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: gga_exc *)

gamma_ss := 0.2:
cc_ss    := [0.0136823, 0.268920, -0.550769,  1.03947, 0.0]:

gamma_ab := 0.006:
cc_ab    := [0.836897,  1.72051,  -2.78498,  -4.57504, 0.0]:

$include "lda_c_vwn.mpl"
$include "b97.mpl"

f := (rs, z, xt, xs0, xs1) ->
  b97_f(f_vwn, gamma_ss, cc_ss, gamma_ab, cc_ab,
        rs, z, xs0, xs1):


text:
# Velocity Calculation for the Robot based on MDH frames
# Introduction
# Berechnung der Geschwindigkeit von Koordinatensystemen und Schwerpunkten
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# floatb -> floating base wird durch base twist (Geschwindigkeit der Basis) oder vollständige Orientierung (Euler-Winkel) berücksichtigt
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# velocity -> Berechnung der Geschwindigkeit aller Segmente
# worldframe -> Berechnung der Geschwindigkeit im Welt-KS (KSW)
# par1 -> Parametersatz 1 (Schwerpunkt als Parameter: SX,SY,SZ)

# Authors
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-03
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
#
# Sources
# [GautierKhalil1990] Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots
# [KhalilDombre2002] Modeling, Identification and Control of Robots
# [Ortmaier2014] Vorlesungsskript Robotik I (WS 2014/15)
# [Ott2008] Cartesian Impedance Control of Redundant and Flexible-Joint Robots
# 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):
codegen_act := true:
codegen_opt := 2:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
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":
read "../robot_codegen_definitions/robot_env":
printf("%s. Generiere Geschwindigkeit für %s (Herleitung im Welt-KS)\n", FormatTime("%Y-%m-%d %H:%M:%S"), 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(4)=-1 then # Standard-Einstellungen:
  if not kin_constraints_exist then # normale serielle Ketten und Baumstrukturen
    use_simplify := 0: # Standardmäßig aus
  else # mit kinematischen Zwangsbedingungen
    use_simplify := 1: # standardmäßig simplify-Befehle anwenden
  end if:
else
  use_simplify := simplify_options(4): # vierter Eintrag ist für Geschwindigkeit
end if:

# 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:
read sprintf("../codeexport/%s/tmp/kinematics_com_worldframe_floatb_%s_par1_maple.m", robot_name, base_method_name):
r_W_W_Si := r_W_W_Si:
r_W_i_Si := r_W_i_Si:
# 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:

# Calculate Velocities
# First assume fixed base model with base velocity and acceleration set to zero
# Anfangsgeschwindigkeit definieren für floating base model
# Velocities of Frames
# Anfangsgeschwindigkeit der Basis
# Analog zu [Ortmaier2014], Gl. (7.6)
# Geschwindigkeit der Basis der Koordinatensysteme der Körper
rD_W_i := Matrix(3, NL):
rD_W_i(1..3,1) := V_base_t(1..3,1):
# Winkelgeschwindigkeit von Körper i ausgedrückt im KS W
omega_W_i := Matrix(3, NL):
if base_method_name = "twist" then:
  omega_W_i(1..3,1) := V_base_t(4..6,1):
end:
if base_method_name = "eulxyz" then:
  omega_W_i(1..3,1) := T_basevel . V_base_t(4..6,1):
end:
printf("%s. Nutze die Methode %s für die Basis-Drehgeschwindigkeit\n", FormatTime("%Y-%m-%d %H:%M:%S"), base_method_name):
# Erhöhe die Geschwindigkeit jedes Körpers
# Betrachte dazu nur die durch Gelenke angetriebenen Körper, nicht die Basis
for i from 1 to NJ do # Gelenke durchgehen
  # Körper der von Gelenkwinkel i bewegt wird: Index i+1
  # Vorgängerkörper bestimmen
  j := v(i) + 1:
  if sigma(i) = 0 then # Drehgelenk
    # Geschwindigkeit des Vorgängers; Trf_c(...,i+1) enthält z-Achse des Körperkoordinatensystems, das von qi bewegt wird.
    # [Ortmaier2014] (7.7) (S.115) (dort falsche Indizes für MDH), [KhalilDombre2002] (9.14)
    omega_W_i(1 .. 3, i+1) := omega_W_i(1 .. 3, j) + thetaD(i,1)*Trf_c(1 .. 3, 3, i+1):
  else: # Schubgelenk
    # [Ortmaier2014] (7.12) (S.116)
    omega_W_i(1 .. 3, i+1) := omega_W_i(1 .. 3, j):
  end if:
  # Terme vereinfachen (Teil 1)
  if use_simplify=1 then
    tmp_t11:=time():
    tmp_l11 := length(omega_W_i(1 .. 3, i+1)):
    omega_W_i(1 .. 3, i+1) := simplify2(omega_W_i(1 .. 3, i+1)):
    tmp_l21 := length(omega_W_i(1 .. 3, i+1)):
    tmp_t21:=time():
  end if:
  # Vektor vom Ursprung des vorherigen Koordinatensystems zu diesem KS
  r_W_im1_i := -Trf_c(1 .. 3, 4, j) + Trf_c(1 .. 3, 4, i+1):
  if sigma(i) = 0 then # Drehgelenk
    # [Ortmaier2014] (7.10) (S.115), [KhalilDombre2002] (9.15)
    rD_W_i(1 .. 3, i+1) := rD_W_i(1 .. 3, j) + CrossProduct(omega_W_i(1 .. 3, j), r_W_im1_i):
  else: # Schubgelenk
    # [Ortmaier2014] (7.15) (S.116), [KhalilDombre2002] (9.15)
    rD_W_i(1 .. 3, i+1) := rD_W_i(1 .. 3, j) + CrossProduct(omega_W_i(1 .. 3, j), r_W_im1_i) + dD(i,1)*Trf_c(1 .. 3, 3, i+1):
  end if:
  printf("%s. Geschwindigkeit für Körperkoordinatensystem %d aufgestellt (Herleitung im Welt-KS).\n",\
    FormatTime("%Y-%m-%d %H:%M:%S"), i-1):#0=Basis
  # Terme vereinfachen (Teil 2)
  if use_simplify=1 then
    tmp_t12:=time():
    tmp_l12 := length(rD_W_i(1 .. 3, i+1)):
    rD_W_i(1 .. 3, i+1) := simplify2(rD_W_i(1 .. 3, i+1)):
    tmp_l22 := length(rD_W_i(1 .. 3, i+1)):
    tmp_t22:=time():
    printf("%s: Terme für Geschwindigkeiten vereinfacht. Länge: %d->%d / %d->%d. Rechenzeit %1.1fs und %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l11, tmp_l21, tmp_l12, tmp_l22, tmp_t21-tmp_t11, tmp_t22-tmp_t12):
  end if:
end do:
# Velocities of Center of Mass
#
rD_W_Si := Matrix(3, NL):
for i to NL do
  # [Ortmaier2014] (7.17) (S.118), [KhalilDombre2002] (9.15)
  rD_W_Si(1 .. 3, i) := Matrix(rD_W_i(1 .. 3, i)) + Matrix(CrossProduct(omega_W_i(1 .. 3, i), r_W_i_Si(1 .. 3, i))):
  printf("%s. Geschwindigkeit für Körperschwerpunkt %d aufgestellt.\n", FormatTime("%Y-%m-%d %H:%M:%S"), i-1):#0=Basis
  if use_simplify=1 then
    tmp_t1:=time():
    tmp_l1 := length(rD_W_Si(1 .. 3, i)):
    rD_W_Si(1 .. 3, i) := simplify2(rD_W_Si(1 .. 3, i)):
    tmp_l2 := length(rD_W_Si(1 .. 3, i)):
    tmp_t2:=time():
    printf("%s: Terme für Schwerpunkts-Geschwindigkeiten vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end if:
end do:

# Exportieren
# Basis-Position aus Termen entfernen. Die Basis-Position kann keinen Einfluss auf die Geschwindigkeiten der einzelnen Körper haben.
# Bei manchen komplexen Systemen scheint der Term aber trotzdem nicht zu verschwinden, da Maple die Vereinfachung nicht erkennt
for i from 1 to 3 do
  omega_W_i := subs({X_base_s[i,1]=0},omega_W_i):
  rD_W_i := subs({X_base_s[i,1]=0},rD_W_i):
  rD_W_Si := subs({X_base_s[i,1]=0},rD_W_Si):
end do:
# Maple Export
save omega_W_i, rD_W_i, rD_W_Si, sprintf("../codeexport/%s/tmp/velocity_worldframe_floatbase_%s_par1_maple.m", robot_name, base_method_name):
printf("%s. Maple-Ausdrücke exportiert.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
# Matlab Export
if codegen_act then
  MatlabExport(convert_t_s(omega_W_i), sprintf("../codeexport/%s/tmp/velocity_omegai0_floatb_%s_worldframe_matlab.m", robot_name, base_method_name), codegen_opt):
  MatlabExport(convert_t_s(rD_W_i), sprintf("../codeexport/%s/tmp/velocity_rDi0_floatb_%s_worldframe_matlab.m", robot_name, base_method_name), codegen_opt):
  MatlabExport(convert_t_s(rD_W_Si), sprintf("../codeexport/%s/tmp/velocity_rDSi0_floatb_%s_worldframe_par1_matlab.m", robot_name, base_method_name), codegen_opt):
  printf("Geschwindigkeiten in Matlab exportiert. %s\n", FormatTime("%Y-%m-%d %H:%M:%S")):
end if:

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


text:
process ACV_handleSuccessfulPortals {
// TODO: playsound
}


text:
# Acceleration Calculation for the Robot based on MDH frames
# Introduction
# Berechnung der Beschleunigung von Koordinatensystemen und Schwerpunkten
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# acceleration_mdh_angles  -> Berechnung der Beschleunigung der MDH-Koordinaten (Drehung und Verschiebung in z-Richtung. Diese können zeitabhängig sein.)
# Sources
# [KhalilDombre2002] Modeling, Identification and Control of Robots
# [Ortmaier2014] Vorlesungsskript Robotik I (WS 2014/15)
# 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):
codegen_act := true:
codegen_opt := 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 Beschleunigung für %s\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(5)=-1 then # Standard-Einstellungen:
  if not kin_constraints_exist then # normale serielle Ketten und Baumstrukturen
    use_simplify := 0: # Standardmäßig aus
  else # mit kinematischen Zwangsbedingungen
    use_simplify := 1: # standardmäßig simplify-Befehle anwenden
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(5): # fünfter Eintrag ist für Beschleunigung
end if:

# Ergebnisse der Geschwindigkeit laden
read sprintf("../codeexport/%s/tmp/velocity_mdh_angles_maple.m", robot_name):
thetaD:= thetaD:
dD:=dD:
# Lade Marker für Existenz kinematischer Zwangsbedingungen (Verknüpfung von MDH-Gelenkwinkeln durch verallgemeinerte Koordinaten)
read sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
kin_constraints_exist := kin_constraints_exist:

# Zeitableitung der Drehwinkelgeschwindigkeit berechnen
# Im Falle kinematischer Zwangsbedingungen wurden diese schon im Arbeitsblatt für die Geschwindigkeit eingesetzt
# und müssen hier nicht mehr betrachtet werden.
# Zeitableitung der Geschwindigkeit in Abhängigkeit der verallgemeinerten Koordinaten
thetaD_qt := thetaD:
dD_qt := dD:
thetaDD := Matrix(NJ, 1):
dDD:= Matrix(NJ, 1):

for i from 1 to NJ do:
  thetaDD(i,1):= diff(thetaD_qt(i,1),t):
  dDD(i,1)    := diff(dD_qt(i,1),t):
end do:
if kin_constraints_exist then # ist nur rechenaufwändig, wenn ZB vorliegen
  printf("%s. Zweite Zeitableitung der MDH-Winkel gebildet.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
end if:
# Terme vereinfachen

if use_simplify>=1 and kin_constraints_exist then # ist nur sinnvoll, wenn ZB vorliegen
  tmp_t0 := time():
  tmp_l0 := length(thetaDD)+length(dDD):
  printf("%s: Vereinfache MDH-Beschleunigungen. Länge: %d.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l0):
  for i from 1 to NJ do
    tmp_t1 := time():
    tmp_l1 := length(thetaDD(i,1)) + length(dDD(i,1)): # es kann sowieso nur einer der beiden Informationen enthalten
    printf("%s: Vereinfache MDH-Beschl. %d. Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i, tmp_l1):
    thetaDD(i,1) := simplify2(thetaDD(i,1)):
    dDD(i,1)     := simplify2(dDD(i,1)):
    tmp_t2 := time():
    tmp_l2 := length(thetaDD(i,1)) + length(dDD(i,1)):
    printf("%s: MDH-Beschl. %d vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end do:
  tmp_l3 := length(thetaDD)+length(dDD):
  printf("%s: MDH-Beschleunigungen vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l0, tmp_l3, tmp_t2-tmp_t0):
end if:
# Ausdruck für Zeitableitungen der Beschleunigungen exportieren
if codegen_act then
  MatlabExport(convert_t_s(thetaDD), sprintf("../codeexport/%s/tmp/acceleration_mdh_angles_matlab.m", robot_name), codegen_opt):
  MatlabExport(convert_t_s(dDD), sprintf("../codeexport/%s/tmp/acceleration_mdh_deltaz_matlab.m", robot_name), codegen_opt):
end if:
# Ausdruck für Maple speichern
save thetaDD, dDD, sprintf("../codeexport/%s/tmp/acceleration_mdh_angles_maple.m", robot_name):
save thetaDD, dDD, sprintf("../codeexport/%s/tmp/acceleration_mdh_angles_maple", robot_name):


text:
`is_element/sphere` := (N::nonnegint) -> proc(x)
 local err,i;

 if not `is_element/R`(N+1)(x) then
  return false;
 fi;

 err := simplify(add(x[i]^2,i=1..N+1) - 1);
 return evalb(err = 0);
end:

`is_equal/sphere` := (N::nonnegint) -> (x,y) -> evalb(simplify(x -~ y) = [0$(N+1)]);

`is_leq/sphere` := NULL;
`list_elements/sphere` := NULL;
`count_elements/sphere` := NULL;

`random_element/sphere` := (N::nonnegint) -> proc(d::posint := 5)
 local r,x,n,k,e,i;

 r := rand(-d..d);
 x := [seq(r(),i=1..N)];
 n := add(x[i]^2,i=1..N);
 x := x *~ (2/(1+n));
 k := rand(0..N)();
 e := (-1) ^ rand(0..1)(); # only really needed for N = 0
 x := e *~ [op(1..k,x),(1-n)/(1+n),op(k+1..N,x)];
 return x;
end:


text:
IRATH:=module()
    option  package;
    export  findTanhSolutions;
    local   nomalizeEqn,
            getODE,
            simplifyEqn,
            findDiffOrder,
            findElemOrder,
            findItemOrder,
            findOrders,
            findInflexion,
            solveAllOrder,
            solveForOrder,
            selectSolutions,
            printSolutions,
            solveForPS,
            transSolution,
            subsIeq,
            subsSol,
            checkSolution,
            isEquationOf;

    (*
     * 求解函数
     *)
    findTanhSolutions:=proc(eqi)
        local eq,ms,params,oeq,st;
        st:=time();

        printf("输入的方程为");
        eq:=nomalizeEqn(eqi);
        oeq:=eq;
        print(eq);

        # printf("方程中的参数为");
        params:=indets(eq,name) minus {t,x};
        # print(params);

        # printf("行波变换后的方程为");
        eq:=getODE(eq);
        # print(eq);

        # printf("方程各项阶数为");
        ms:=findOrders(eq);
        # print(ms);

        # printf("拐点为");
        ms:=findInflexion(ms);
        # print(ms);

        if evalb(ms={0}) then
            printf("无解\n");
            printf("时间已过 %f 秒\n",time()-st);
            return;
        end if;

        # printf("开始求解\n");
        solveAllOrder(eq,ms,params,oeq);

        printf("时间已过 %f 秒\n",time()-st);
    end proc:

    (*
     * 标准化方程
     * 方程变量是u(x,t)
    *)
    nomalizeEqn:=proc(ee)
        local e;
        e:=expand(numer(ee));
        return e;
    end proc:

    (*
     * 进行行波变换，将偏微分方程转化为常微分方程
     * 要求函数为u(x,t),替换为u(xi),
     * xi=k*(x-ct)
    *)
    getODE:=proc(eqi)
        local eq;
        eq:=subs(u(x,t)=u(x,t,xi),eqi);
        eq:=PDETools[dsubs](diff(u(x,t,xi),t)=-c*k*diff(u(x,t,xi),xi),
        diff(u(x,t,xi),x)=k*diff(u(x,t,xi),xi),
        eq);
        eq:=eval(subs(u(x,t,xi)=u(xi),eq));
        return simplifyEqn(eq);
    end proc:

    (*
     * 化简方程
    *)
    simplifyEqn:=proc(ee)
        local e,v,p;
        e:=numer(ee);
        e:=factor(e);
        if type(e,`*`) then
            p:=1;
            for v in e do
                if type(v,`+`) then
                    p:=p*v;
                 elif type(v,`^`) and type(expand(op(1,v)),`+`) then
                    p:=p*(simplifyEqn(op(1,v))^op(2,v));
                end if;
            end do;
            p:=expand(p);
        else
            p:=expand(e);
        end if;
        return p;
    end proc:

    (*
     * 计算导数项阶数
    *)
    findDiffOrder:=proc(ee)
        local e;
        e:=convert(ee,'D');
        e:=op(0,op(0,e));
        if evalb(op(0,e)=`@@`) then
            op(2,e);
        else
            1;
        end if;
    end proc:

    (*
     * 计算元素的阶数
    *)
    findElemOrder:=proc(ee)
        local e;
        if evalb(op(0,ee)='u') then
            e:=[1,0];
        elif evalb(op(0,ee)='diff') then
            e:=[1,findDiffOrder(ee)];
        elif type(ee,`^`) then
            e:=findElemOrder(op(1,ee))*op(2,ee);
        else
            e:=[0,0];
        end if;
        return e;
    end proc:

    (*
     * 计算方程某一项的阶数
    *)
    findItemOrder:=proc(ee)
        local e,p;
        if type(ee,`*`) then
            p:=[0,0];
            for e in ee do
                p:=p+findElemOrder(e);
            end do;
        else
            p:=findElemOrder(ee);
        end if;
        return p;
    end proc:

    (*
     * 计算方程每一项的阶数
     * 假设至少有两项
     * 假设指数上不含参数
    *)
    findOrders:=proc(ee)
        findItemOrder~([op(ee)]);
    end proc:

    (*
     * 计算拐点
    *)
    findInflexion:=proc(mss::list(list(integer)))
        local ms,i,j,n,mm,m,fun;
        fun:=(x,m)->x[1]*m+x[2];
        mm:={};
        ms:={op(mss)};
        n:=numelems(ms);
        for i from 1 to n-1 do
            for j from i+1 to n do
                if evalb(ms[i][1]<>ms[j][1]) then
                    m:=-(ms[i][2]-ms[j][2])/(ms[i][1]-ms[j][1]);
                    if evalb(fun(ms[i],m)=max(map(fun,ms,m))) then
                        mm:=mm union {m};
                    end if;
                end if;
            end do;
        end do;
        return mm;
    end proc:

    (*
     * 对所有阶数进行求解
    *)
    solveAllOrder:=proc(eqi,mms::set(rational),params::set(name),oeq)
        local trans,tr,eq,ms,m;
        # 计算全部变换
        trans:=map(x->sign(x)/denom(x),mms);
        for tr in trans do
            eq:=eval(subs(u(xi)=u(xi)^tr,eqi));
            eq:=simplifyEqn(eq);
            ms:=findOrders(eq);
            m:=select(type,findInflexion(ms),posint);
            m:=max(m);
            printf("m=%a,方程为:",tr*m);
            print(eq);
            if not type(ms,list(list(integer))) then
                printf("方程不是多项式，无解\n");
                next;
            end if;
            solveForOrder(eq,m,params,tr,oeq);
        end do;
    end proc:

    (*
     * 对特定阶数进行求解
    *)
    solveForOrder:=proc(eqi,m::posint,params::set(name),tr::rational,oeq)
        local eq,f,PS,vars,sols;
        f:=add(seq(a[i]*tanh(xi)^i,i=0..m));
        vars:=[seq(a[i],i=0..m)];
        eq:=eval(subs(u(xi)=f,eqi));
        eq:=subs(tanh(xi)=T,eq);
        eq:=collect(eq,T);
        PS:=PolynomialTools[CoefficientList](eq,T,termorder=reverse);
        sols:=solveForPS(PS,m,params,vars);
        printSolutions(sols,m,tr);
    end proc:

    (*
     * 筛选解
     * 删去k,c=0的解
     * 删除参数为0的解
     * 删除a[1]..a[m]都为0的解
    *)
    selectSolutions:=proc(solss::{set,list},
                         params::set(name),m::posint)
        local zs,ps,sols,as;
        as:={seq(a[i],i=1..m)};
        sols:=convert(solss,set);
        sols:=select(type,sols,equation);
        ps:=params union {k,c};
        zs:=lhs~(select(_x->type(_x,`=`(name,0)),sols));
        if evalb( (zs intersect ps) <> {} ) then
            return false;
        else
            return not evalb(as subset zs);
        end if;
    end proc:

    (*
     * 输出解
    *)
    printSolutions:=proc(sols::set({list,set}),
                        m::posint,tr::rational)
        local sol;
        if evalb(sols={}) then
            printf("无解\n");
            return;
        end if;
        printf("共有%d个解",numelems(sols));
        print(u(xi)=add(seq(a[i]*tanh(xi)^i,i=0..m))^tr,xi=k*(x-ct)+xi[0]);
        printf("其中");
        for sol in sols do
            print(expand(simplify(sol)));
        end do;
    end proc:

    (*
     * 求解PS
     * 基于csolve的版本
     * 只有剩余方程不能用csolve求解时，才使用solve求解
    *)
    solveForPS:=proc(PS,m::posint,params::set(name),vars::list(name))
        local pa,pp,sola,sol,vs,res,np,nsol,ssol,nres,sssol;
        res:={};# 所有解
        pa:=PS[1..(m+1)];
        pp:=PS[(m+1)..-1];
        sola:=[csolve(pa,{vars[],k,c})];# 求解前m+1个方程
        sola:=select(selectSolutions,sola,params,m);# 除去平凡的解
        for sol in sola do
            # 代入化简
            np:=simplify(subs(op(sol),pp));
            np:=simplifyEqn~(np);# 化简方程，去除非零项
            np:=remove(type,np,0);

            # 如果剩余的方程全0，则不需要继续求解
            if evalb(np=[]) then
                res:=res union {sol};
                next;
            end if;

            # 求解剩余方程
            try
                nsol:=[csolve(np,indets(np,name))];
            catch:
                nsol:=[RealDomain[solve](np,indets(np,name))];
            end try;

            nsol:=select(selectSolutions,nsol,params,m);

            for ssol in nsol do
                # 合并解，因为前面的解，在后面已经被带入
                # 所以，后面的解不会包含前面已有的符号
                # 因此，直接合并即可
                try # 有可能出现分母为0的情况
                    sssol:=simplify(subs(op(ssol),sol));
                    res:=res union {sssol union ssol};
                catch:
                    next;
                end try;
            end do;
        end do;

        res:=select(selectSolutions,res,params,m);

        # printf("未转换的解");
        # print~(res);

        nres:={};
        map(transSolution,res,'nres',params,vars);
        res:=select(selectSolutions,nres,params,m);

        return res;
    end proc:

    (*
     * 转化csolve的结果
    *)
    transSolution:=proc(res::set(equation),nres::evaln(set),
                       params::set(name),vars::list(name))
        local r,rv,rp,np,rrv,rrp,sp,newSol;
        # 先求解a[0]..a[m],k,c
        sp:=select(isEquationOf,res,{vars[],k,c});
        sp:=map(_x->simplifyEqn(lhs(_x)-rhs(_x)),sp);# 化简方程，去除非零项
        rv:=RealDomain[solve](sp,[vars[],k,c]);

        # 如果求解失败则直接原样返回
        if evalb(rv=[]) then
            nres:=eval(nres) union {[res[]]};
            WARNING("solve求解失败");
            return;
        end if;

        # 再求解参数
        for rrv in rv do
            try
                np:=simplify(subs(op(rrv),res));
            catch:
                next;
            end try;
            np:=map(_x->simplifyEqn(lhs(_x)-rhs(_x)),np);# 化简方程，去除非零项
            np:=select(isEquationOf,np,params);
            rp:=RealDomain[solve](np,[params[]]);
            for rrp in rp do
                rrp:=subsIeq~(rrp);
                try
                    newSol:=[op(subsSol(rrv,rrp)),rrp[]];
                    newSol:=simplify(rationalize(newSol));
                    nres:=eval(nres) union {newSol};
                catch:
                    next;
                end try;
            end do;
        end do;
        return;
    end proc:

    (*
     * 消去不等式范围约束
    *)
    subsIeq:=proc(ee)
        local v;
        if type(ee,equation) then
            return ee;
        else
            v:=indets(ee);
            if evalb(nops(v)=1) then
                v:=v[1];
                return (v=v);
            else
                return ee;
            end if;
        end if;
    end proc:

    (*
     * 解集带入
     * sol2代入sol1
    *)
    subsSol:=proc(sol1,sol2)
        local sol;
        sol:=select(type,sol2,equation);
        expand(subs(op(sol),sol1));
    end proc:

    (*
     * 代回原方程进行验证
    *)
    checkSolution:=proc(ssol,oeq,m,tr)
        local uxi,sol,r;
        sol:=select(type,ssol,equation);
        r:=subs(u(x,t)=add(seq(a[i]*tanh(xi)^i,i=0..m))^tr,oeq);
        r:=subs(xi=k*(x-c*t),r);
        r:=subs(op(sol),r);
        r:=simplify(r);
        evalb(r=0);
    end proc:

    (*
     * 选择关于给定参数的方程
    *)
    isEquationOf:=proc(eq,vars::set)
        local vs;
        vs:=indets(eq,name);
        evalb( (vs intersect vars) <> {} );
    end proc:

end module:


text:
######################################################################
#
# General operad framework.

`check_gamma_args` := proc(A,B,p,U,V,is_el_CC)
 local F,b;
 global reason;

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

 if not `is_element/maps`(A,B)(p) then
  reason := [convert(procname,string),"not a map from A to B",A,B,reason];
  return false;
 fi;

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

 if not is_el_CC(B)(U) then
  reason := [convert(procname,string),"U is not in C(B)",eval(U),B];
  return false;
 fi;

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

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

 for b in B do
  if not is_el_CC(F[b])(V[b]) then
   reason := [convert(procname,string),"V[b] is not in C(F[b])",b,eval(V[b]),eval(F[b])];
   return false;
  fi;
 od;

 return true;
end;

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

`check_gamma_axiom` := proc(A,B,C,p,q,U,V,W,is_el_CC,is_equal_CC,gamma_CC)
# local a,b,c,qp,Fp,Fq,Fqp,UV,VW,UVW0,UVW1,W0,p0;
# global reason,check_gamma_axiom_args;
 local a,b,c,qp,Fp,Fq,Fqp,UV,VW,W0,p0;
 global UVW0,UVW1,reason,check_gamma_axiom_args;

 check_gamma_axiom_args := [args];

 if not type(A,set) then
  reason := [convert(procname,string),"A is not a set",eval(A)];
  return false;
 fi;
 if not type(B,set) then
  reason := [convert(procname,string),"B is not a set",eval(B)];
  return false;
 fi;
 if not type(C,set) then
  reason := [convert(procname,string),"C is not a set",eval(C)];
  return false;
 fi;

 if not `is_element/maps`(A,B)(p) then
  reason := [convert(procname,string),"p is not a map from A to B",eval(p),A,B,reason];
  return false;
 fi;

 if not `is_element/maps`(B,C)(q) then
  reason := [convert(procname,string),"q is not a map from B to C",eval(q),B,C,reason];
  return false;
 fi;

 qp := `compose/maps`(A,B,C)(p,q);

 Fp := fibres(A,B)(p);
 Fq := fibres(B,C)(q);
 Fqp := fibres(A,C)(qp);

 if not is_el_CC(C)(U) then
  reason := [convert(procname,string),"U is not in O(C)",eval(U),C];
  return false;
 fi;

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

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

 for c in C do
  if not is_el_CC(Fq[c])(V[c]) then
   reason := [convert(procname,string),"V[c] is not in O(Fq[c])",c,eval(V[c]),eval(Fq[c])];
   return false;
  fi;
 od;

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

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

 for b in B do
  if not is_el_CC(Fp[b])(W[b]) then
   reason := [convert(procname,string),"W[b] is not in O(Fp[b])",b,eval(W[b]),eval(Fp[b])];
   return false;
  fi;
 od;

 UV := gamma_CC(B,C)(q)(U,V);
 if not is_el_CC(B)(UV) then
  reason := [convert(procname,string),"UV is not in O(B)",eval(UV),B,C,eval(q),eval(U),eval(V)];
  return false;
 fi;

 if not `check_gamma_args`(A,B,p,UV,W,is_el_CC) then
  reason := [convert(procname,string),"invalid arguments to gamma 1",reason];
  return false;
 fi;

 UVW0 := gamma_CC(A,B)(p)(UV,W);
 if not is_el_CC(A)(UVW0) then
  reason := [convert(procname,string),"UVW0 is not in C(A)",eval(UVW0),A];
  return false;
 fi;

 VW := table();
 for c in C do
  W0 := table();
  p0 := table();
  for b in Fq[c] do W0[b] := W[b]; od;
  for a in Fqp[c] do p0[a] := p[a]; od;
  if not `check_gamma_args`(Fqp[c],Fq[c],p0,V[c],W0,is_el_CC) then
   reason := [convert(procname,string),"invalid arguments to gamma 2",reason];
   return false;
  fi;
  VW[c] := gamma_CC(Fqp[c],Fq[c])(p0)(V[c],W0);
 od;

 if not `check_gamma_args`(A,C,qp,U,VW,is_el_CC) then
  reason := [convert(procname,string),"invalid arguments to gamma 3",reason];
  return false;
 fi;

 UVW1 := gamma_CC(A,C)(qp)(U,VW);
 if not is_el_CC(A)(UVW1) then
  reason := [convert(procname,string),"UVW1 is not in C(A)",eval(UVW1),A];
  return false;
 fi;

 if not is_equal_CC(A)(UVW0,UVW1) then
  reason := [convert(procname,string),"UVW0 != UVW1",eval(UVW0),eval(UVW1)];
  return false;
 fi;

 return true;
end;

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

check_operad_morphism_axiom :=
 proc(A,B,p,U,V,
      is_el_CC,is_equal_CC,gamma_CC,
      is_el_DD,is_equal_DD,gamma_DD,
      phi)

 local b,Fp,UV,U1,V1,UV1,UV2;

 global reason,check_operad_morphism_axiom_args;

 check_operad_morphism_axiom_args := [args];

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

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

 if not `is_element/maps`(A,B)(p) then
  reason := [convert(procname,string),"p is not a map from A to B",eval(p),A,B,reason];
  return false;
 fi;

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

 if not is_el_CC(B)(U) then
  reason := [convert(procname,string),"U is not in C(B)",eval(U),B];
  return false;
 fi;

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

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

 for b in B do
  if not is_el_CC(Fp[b])(V[b]) then
   reason := [convert(procname,string),"V[b] is not in C(Fp[b])",b,eval(V[b]),eval(Fp[b])];
   return false;
  fi;
 od;

 UV := gamma_CC(A,B)(p)(U,V);
 if not is_el_CC(A)(UV) then
  reason := [convert(procname,string),"UV is not in C(A)",eval(UV),A,B,eval(p),eval(U),eval(V)];
  return false;
 fi;

 U1 := phi(B)(U);
 if not is_el_DD(B)(U1) then
  reason := [convert(procname,string),"phi(U) is not in D(B)",eval(U1),B];
  return false;
 fi;

 V1 := table();
 for b in B do
  V1[b] := phi(Fp[b])(V[b]);
  if not is_el_DD(Fp[b])(V1[b]) then
   reason := [convert(procname,string),"phi(V[b]) is not in D(Fp[b])",eval(V1),b,eval(Fp[b])];
   return false;
  fi;
 od;

 UV1 := gamma_DD(A,B)(p)(U1,V1);
 if not is_el_DD(A)(UV1) then
  reason := [convert(procname,string),"UV1 is not in D(A)",eval(UV1),A];
  return false;
 fi;

 UV2 := phi(A)(UV);
 if not is_el_DD(A)(UV1) then
  reason := [convert(procname,string),"phi(UV) is not in D(A)",eval(UV2),A];
  return false;
 fi;

 if not(is_equal_DD(A)(UV1,UV2)) then
  reason := [convert(procname,string),"UV1 != UV2",eval(UV1),eval(UV2)];
  return false;
 fi;

 return true;
end:

check_operad_morphism_axiom_globalise := proc()
 global check_operad_morphism_axiom_args,
         A,B,p,U,V,
         is_el_CC,is_equal_CC,gamma_CC,
         is_el_DD,is_equal_DD,gamma_DD,
         phi;
 local ARGS;
 ARGS := check_operad_morphism_axiom_args;

 A               := eval(ARGS[ 1]);
 B		 := eval(ARGS[ 2]);
 p		 := eval(ARGS[ 3]);
 U		 := eval(ARGS[ 4]);
 V		 := eval(ARGS[ 5]);
 is_el_CC	 := eval(ARGS[ 6]);
 is_equal_CC	 := eval(ARGS[ 7]);
 gamma_CC	 := eval(ARGS[ 8]);
 is_el_DD	 := eval(ARGS[ 9]);
 is_equal_DD	 := eval(ARGS[10]);
 gamma_DD	 := eval(ARGS[11]);
 phi 		 := eval(ARGS[12]);
 NULL;
end:

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

# An operad E consists of sets E(A) for each finite set A, functorial
# for bijections of A.  We have formulated the definitions in such a
# way that this functoriality is determined by the rest of the
# structure.  The code here implements this.

`action/operad` := (gen_name,eta,gamma) -> (A::set,B::set) -> (f) -> proc(x)
 local f_inv,e,a;

 f_inv := table();
 e := table();
 for a in A do
  f_inv[f[a]] := a;
  e[a] := eta({f[a]});
 od;
 return gamma(B,A)(f_inv)(x,e);
end;

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

extend_gamma_linear := (g0) -> (A::set,B::set) -> (p) -> proc(y,x)
 local c,y1,F,BL,M,b,t,z,m,mt,u,i;

 if y = 0 then
  return 0;
 elif type(y,`+`) then
  return map(t -> extend_gamma_linear(g0)(A,B)(p)(t,x),y);
 fi;

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

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

 BL := sort([op(B)]);
 M := [[1]];
 for b in BL do
  t := map(coeff_split,sum_terms(x[b]));
  M := [seq(seq([m[1] * u[1],seq(m[i],i=2..nops(m)),u[2]],u in t),m in M)];
 od:

 z := 0;

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

 return z;

end;


text:
est_pi := proc(num_steps::integer)
   local step, s, x, mypi, i, start_time, stop_time:
   printf("  Estimating Pi using: %1d slices\n", num_steps):
   start_time := time[real]():
   step := 1.0/num_steps:
   s:=0:
   for i from 1 to num_steps do
      x:=(i-0.5)*step:
      s:=s+(4.0/(1.0 + x*x)):
   end do:
   mypi:=s*step:
   stop_time := time[real]():
   printf("Estimated value of Pi: %.16g\n", mypi):
   printf("      Elapsed time(s): %.16g\n", stop_time - start_time):
   printf("  Observed difference: %.16g\n", abs(Pi - mypi)):
   mypi;
end proc:

native_est_pi := proc(num_steps::integer)
   local step, mypi, k, start_time, stop_time:
   printf("  Estimating Pi using: %1d slices\n", num_steps):
   start_time := time[real]():
   mypi:=Re(sum((4/(1+(((k-0.5)/num_steps)^2))),k=1..num_steps)/num_steps):
   stop_time := time[real]():
   printf("Estimated value of Pi: %.16g\n", mypi):
   printf("      Elapsed time(s): %.16g\n", stop_time - start_time):
   printf("  Observed difference: %.16g\n", abs(Pi - mypi)):
   mypi;
end proc: