text:
#int fact (int n) {
#  if(n != 1)
#    return foo(n - 1);
#  else return 1;
#}

func $fact (
 var %n i32) i32 {
 if (ne i32 i32 (dread i32 %n, constval i32 1)) {
     return (constval i32 0) }
 else {
     return (constval i32 1) } }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

ext: 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) 2021 Susi Lehtola

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

(* type: mgga_exc *)

$include "mgga_x_tm.mpl"

reg_c := 3:
reg_d := 1.475:

(* text after equation 11 *)
regtm_f := (a, p) -> (1 - a)^3/(1 + (reg_d*a)^2)^(3/2) * exp(-reg_c*p):

(* equation 11 *)
regtm_zp := (a, p) -> 1 / (1 + 3/5*(a/(p+regtm_f(a,p)))):

(* equation 13 *)
regtm_w := zp -> (zp^2 + 3*zp^3) / (1 + zp^3)^2:

(* Collect all the pieces together *)
tm_w := (x, t) -> regtm_w(regtm_zp(tm_alpha(x, t), tm_p(x))):


text:
######################################################################
# Fbar(N)(A) is the Fulton-MacPherson compactification of F(N)(A),
# defined in terms of coherent families as in the LaTeX document.

`is_element/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 local TT1,B,C,P,Q,T,U;
 global reason;

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

 P := sort(map(sort,map(op,{indices(x)})));
 TT1 := select(T -> nops(T) > 1,TT);

 if P <> TT1 then
  reason := [convert(procname,string),"x is not indexed by the big sets in TT",P,TT1];
  return false;
 fi;

 C := children_map(A)(TT);

 for T in P do
  if not(`is_element/SW`(N)(T)(x[T])) then
   reason := [convert(procname,string),"x[T] is not in SW(N)(T)",eval(x[T]),N,T,reason];
   return false;
  fi;

  for U in C[T] do
   if nops(U) > 1 and not(`is_element/SWW`(N)(T,U)([x[T],x[U]])) then
    reason := [convert(procname,string),"(x[T],x[U]) is not in SWW(N)(T,U)",eval(x[T]),eval(x[U]),N,T,U,reason];
    return false;
   fi;
  od;
 od;

 return true;
end;

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

`random_element/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc()
 local TT1,children,x,C,U,m,c,a,T;

 TT1 := select(T -> nops(T) > 1,TT);
 children := children_map(A)(TT);
 TT1 := sort([op(TT1)],(T,U) -> nops(T) < nops(U));

 x := table();

 for T in TT1 do
  x[T] := `random_element/vector_functions`(N)(T)();
  C := select(U -> nops(U) > 1,children[T]);
  for U in C do
   m := rand(0..6)()/rand(1..10)();
   c := `random_element/R`(N)();
   for a in U do
    x[T][a] := m *~ x[U][a] +~ c;
   od;
  od;
 od;

 return eval(x);
end;

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

`normalise/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 local x0,T,P;

 x0 := table();
 P := select(T -> nops(T) > 1,TT);

 for T in P do x0[T] := `normalise/SW`(N)(T)(x[T]); od;
 return eval(x0);
end;

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

`bottom_normalise/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 local x0,T,P;

 x0 := table();
 P := select(T -> nops(T) > 1,TT);

 for T in P do x0[T] := `bottom_normalise/SW`(N)(T)(x[T]); od;
 return eval(x0);
end;

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

`is_equal/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x,y)
 local T,P;
 global reason;

 P := select(T -> nops(T) > 1,TT);

 for T in P do
  if not(`is_equal/SW`(N)(T)(x[T],y[T])) then
   reason := [convert(procname,string),"x[T] <> y[T]",T,x[T],y[T],reason];
   return false;
  fi;
 od;

 return true;
end;

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

`dist/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x,y)
 local T,TT1,d;
 global reason;

 TT1 := select(T -> nops(T) > 1,TT);

 d := 0;

 for T in TT1 do
  d := d + `dist/SW`(N)(T)(x[T],y[T]);
 od;

 return d;
end;

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

`is_leq/tree_Fbar` := NULL;

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

`is_interior/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 local V;

 V := map(a -> x[A][a],A);
 return evalb(nops(V) = nops(A));
end;

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

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

 x := table();
 P := select(T -> nops(T) > 1,TT);
 for T in P do
  x[T] := table();
  for t in T do
   x[T][t] := u[t];
  od;
 od;

 return eval(x);
end:

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

`res/tree_Fbar/F` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 if not(`is_interior/tree_Fbar`(N)(A)(TT)(x)) then
  return FAIL;
 fi;

 return x[A];
end;

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

`res/Fbar/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 local TT1,y,T;

 TT1 := select(T -> nops(T) > 1,TT);

 y := table():
 for T in TT1 do y[T] := eval(x[T]); od;

 return eval(y);
end:

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

`ext/tree_Fbar/Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 local TT1,UU,TU,T,U,m,v,ok,a,y;

 TT1 := select(T -> nops(T) > 1,TT);
 UU := `list_elements/big_subsets`(A);

 for U in UU do
  TU := select(T -> U minus T = {},TT1);
  m := min(op(map(nops,TU)));
  TU := select(U -> nops(U) = m,TU);
  T := TU[1];
  y[U] := table();
  v := x[T][U[1]];
  ok := false;
  for a in U do
   y[U][a] := x[T][a];
   if not(`is_equal/R`(N)(x[T][a],v)) then
    ok := true;
   fi;
  od:

  if not(ok) then return FAIL; fi;
 od;

 return eval(y);
end:

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

`perturb/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x,epsilon)
 local UU,n,y,i,a,f;

 UU := sort([op(TT)],(T,U) -> nops(T) > nops(U));
 UU := select(U -> nops(U) > 1,UU);
 n := nops(UU);
 y := table();
 for i from 1 to n do
  y[i] := `normalise/SW`(N)(UU[i])(x[UU[i]]);
  for a in A minus UU[i] do
   y[i][a] := [0$N];
  od;
 od:
 f := table():
 for a in A do
  f[a] := [0$N];
  for i from 1 to n do
   f[a] := combine(simplify(expand(f[a] +~ epsilon^(i-1) *~ y[i][a])));
  od:
 od:

 return `res/Fbar/tree_Fbar`(N)(A)(TT)(`inc/F/Fbar`(N)(A)(f));
end:

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

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


text:
# D/Fbar(A,B) is
# D(A,NWB) = {(u,v) in Map(A,[0,infinity)) x NWB : |u|^2 + |v|^2 = 1}

`is_element/D/Fbar` := (N::posint) -> (A::set,B::set) -> proc(uv)
 global reason;
 local u,v,nu,nv;

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

 u,v := op(uv);

 if not (`is_element/real_functions`(A)(u) and
         `is_nonnegative/real_functions`(A)(u)) then
  reason := [convert(procname,string),"u is not a nonnegative real function on A",u,A];
  return false;
 fi;

 nu := `norm/real_functions`(A)(u)^2;

 if not `is_element/prime_W`(N)(B)(v) then
  reason := [convert(procname,string),"v is not in NWB",v,B];
  return false;
 fi;

 nv := `norm/prime_W`(N)(B)(v)^2;

 if nu + nv <> 1 then
  reason := [convert(procname,string),"(u,v) is not in normalised",u,v];
  return false;
 fi;

 return true;
end:

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

`is_equal/D/Fbar` := (N::posint) -> (A::set,B::set) -> proc(uv1,uv2)
 local u1,v1,u2,v2;
 u1,v1 := op(uv1);
 u2,v2 := op(uv2);

 return `is_equal/real_functions`(A)(u1,u2) and
        `is_equal/vector_functions`(N)(B)(v1,v2);
end:

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

`random_element/D/Fbar` := (N::posint) -> (A::set,B::set) -> proc()
 local n,r,u,v,a,b;

 n := nops(A);
 r := 0;

 while r = 0 do
  u := `random_element/real_functions`(A)();
  v := `random_element/prime_W`(N)(B)();
  r := `norm/real_functions`(A)(u)^2 +
       `norm/vector_functions`(N)(B)(v)^2;
 od;

 r := sqrt(r);
 for a in A do u[a] := simplify(abs(u[a] / r)); od:
 for b in B do v[b] := simplify(v[b] /~ r); od:

 return [eval(u),eval(v)];
end:

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

`is_leq/D/Fbar` := NULL:
`list_elements/D/Fbar` := NULL:
`count_elements/D/Fbar` := NULL:


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

mbeef_k  := 6.5124: (* PBEsol transformation *)
mbeef_xi := p -> 2*p/(mbeef_k + p) - 1:

(* expansion is better than 1e-13 for large a *)
mbeef_xj0 := a -> - (1 - a^2)^3/(1 + a^3*(1 + a^3)):
mbeef_xj := a -> enforce_smooth_lr(mbeef_xj0, a, 1e4, 4):

with(orthopoly):
mbeef_expansion := (x, t) -> add(add(
  + mbeef_coefs[i][j]
  * P(j-1, mbeef_xi(X2S^2*x^2))
  * P(i-1, mbeef_xj((t - x^2/8)/K_FACTOR_C)),
i=1..mbeef_n), j=1..mbeef_n):