7#include "ginac/parse_context.h"
16 alignas(2)
static ex SP_reader(
const exvector& ev) {
17 if(ev.size()<2)
return SP(ev[0]).subs(
SP_map);
18 else return SP(ev[0], ev[1]).subs(
SP_map);
21 alignas(2)
static ex LC_reader(
const exvector& ev) {
22 if(ev.size()==4)
return (-I)*
LC(ev[0], ev[1], ev[2], ev[3]);
26 alignas(2)
static ex SUNT_reader(
const exvector& ev) {
28 if(n==3)
return SUNT(ev[0], ev[1], ev[2]);
31 for(
int i=0; i<n-2; i++)
as.append(ev[i]);
32 return SUNT(
as,ev[n-2],ev[n-1]);
33 }
else throw Error(
"SUNT_reader: number of arguments less than 3.");
36 alignas(2)
static ex TTR_reader(
const exvector& ev) {
37 if(ev.size()==1)
return TTR(ev[0]);
39 for(
auto item : ev)
as.append(item);
43 alignas(2)
static ex DGamma_reader(
const exvector& ev) {
45 if(n==1)
return GAS(ev[0]);
48 for(
int i=1; i<n; i++) ret = ret *
GAS(ev[i],rl);
52 alignas(2)
static ex gi_reader(
const exvector& ev) {
54 if(n==0)
return GAS(1);
55 else if(n==1)
return GAS(1,
ex2int(ev[0]));
56 else throw Error(
"DGamma_reader: number of arguments more than 2.");
59 alignas(2)
static ex GMat_reader(
const exvector& ev) {
60 return GMat(ev[0],ev[1],ev[2]);
69 struct mapGamma :
public map_function {
71 ex operator()(
const ex &e) {
72 if(is_a<DGamma>(e))
return DGamma(ex_to<DGamma>(e), gline);
74 else return e.map(*
this);
76 mapGamma(
int _gline) : gline(_gline) { };
81 struct mapTR :
public map_function {
83 ex operator()(
const ex &e) {
87 trs = mapGamma(gline)(trs);
88 trs = DGamma(1, gline) * trs;
91 if(glmax>128)
throw Error(
"too large index with glmax>128.");
94 }
else if(is_a<add>(e)) {
97 unsigned glmax = gline;
100 res += (*this)(item);
101 if(glmax<gline)
glmax = gline;
105 }
else return e.map(*
this);
112 string init_script = R
"EOF(
113CFunction pow,sqrt,gamma,HF,GMat,WF;
114Tensor TTR(cyclic), f(antisymmetric), T, f4, colTp;
115Symbols reX,I2R,NF,NA,d,I,Pi;
116AutoDeclare Symbols gCF, trcN;
118AutoDeclare Index colA;
120AutoDeclare Index colF;
126 id,once,TTR(?a) = T(?a,colF1,colF1);
129 id,once,T(colA1?,colA2?,?a,colF1?,colF2?) = T(colA1,colF1,colF3)*T(colA2,?a,colF3,colF2);
131 id,once,f4(colA1?,colA2?,colA3?,colA4?) = f(colA1,colA2,colA5) * f(colA5,colA3,colA4);
137 if ( count(f,1) || match(T(colA1?,colF1?,colF2?)*T(colA1?,colF3?,colF4?)) ) redefine colXi "0";
138 id,once,f(colA1?,colA2?,colA3?) = 1/I2R/i_*T(colA1,colF1,colF2)*T(colA2,colF2,colF3)*T(colA3,colF3,colF1)-1/I2R/i_*T(colA3,colF1,colF2)*T(colA2,colF2,colF3)*T(colA1,colF3,colF1);
139 sum colF1,colF2,colF3;
140 id T(colA1?,colF1?,colF2?)*T(colA1?,colF3?,colF4?) = colTp(colF1,colF2,colF3,colF4);
142 if ( count(colTp,1) ) redefine colXj "0";
144 id,once,colTp(colF1?,colF2?,colF3?,colF4?) = I2R*(d_(colF1,colF4)*d_(colF2,colF3)-d_(colF1,colF2)*d_(colF3,colF4)/NF);
149 id T(colA1?,?a,colF1?,colF2?)*T(colA2?,?b,colF2?,colF3?) = T(colA1,?a,colA2,?b,colF1,colF3);
151id T(?a,colF1?,colF1?) = TTR(?a);
153id TTR(colA1?,colA2?) = I2R*d_(colA1,colA2);
165 ex runform(
const ex &expr_in,
int verb) {
166 static map<pid_t, Form> form_map;
168 if((form_map.find(pid)==form_map.end())) {
170 ss << init_script << endl;
171 form_map[pid].Init(
"form");
172 form_map[pid].Execute(ss.str());
174 Form &fprc = form_map[pid];
176 ex expr = expr_in.subs(
SP_map);
181 ids <<
"id " << kv.first <<
"=" << kv.second <<
";" << endl;
182 all_expr += iWF(kv.first) * iWF(kv.second);
184 string idstr = sss.str();
186 lst vec_lst, VD_lst, CF_lst, CA_lst, sym_lst;
187 for(const_preorder_iterator i = all_expr.preorder_begin(); i != all_expr.preorder_end(); ++i) {
188 if(is_a<Vector>(*i)) vec_lst.append(*i);
189 else if(is_a<Index>(*i)) {
190 if(ex_to<Index>(*i).dim==
d) VD_lst.append(*i);
191 else if(ex_to<Index>(*i).dim==
NF) CF_lst.append(*i);
192 else if(ex_to<Index>(*i).dim==
NA) CA_lst.append(*i);
193 }
else if(is_a<symbol>(*i)) sym_lst.append(*i);
194 else if(is_a<GiNaC::function>(*i)) {
195 static vector<string> fun_vec = {
196 "iWF",
"WF",
"TR",
"sin",
"cos",
"HF",
"TTR",
"GMat"
198 auto func = ex_to<GiNaC::function>(*i).get_name();
200 for(
auto fi : fun_vec) {
201 if(fi==func) { ok =
true;
break; }
204 cout << (*i) << endl;
205 throw Error(
"runform: Functions not defined in FORM: "+func);
209 vec_lst.sort(); vec_lst.unique();
210 VD_lst.sort(); VD_lst.unique();
211 CF_lst.sort(); CF_lst.unique();
212 CA_lst.sort(); CA_lst.unique();
214 sym_lst.sort(); sym_lst.unique();
218 ff <<
".store" << endl;
220 if(sym_lst.nops()>0) {
222 for(
auto sx : sym_lst) {
223 auto s = ex_to<symbol>(sx);
224 ff <<
" " << s.get_name();
225 st[s.get_name()] = sx;
229 if(vec_lst.nops()>0) {
231 for(
auto vx : vec_lst) {
232 auto v = ex_to<Vector>(vx);
234 st[v.name.get_name()] = v;
235 for(
auto vvx : vec_lst) {
236 auto vv = ex_to<Vector>(vvx);
237 st[v.name.get_name()+
"__"+vv.name.get_name()] =
SP(v,vv).subs(
SP_map);
242 if(VD_lst.nops()>0) {
244 else ff <<
"Dimension d;" << endl;
246 for(
auto ix : VD_lst) {
247 auto i = ex_to<Index>(ix);
249 st[i.name.get_name()] = i;
255 if(CF_lst.nops()>0) {
256 ff <<
"Dimension NF;" << endl;
258 for(
auto ix : CF_lst) {
259 auto i = ex_to<Index>(ix);
261 st[i.name.get_name()] = i;
267 if(CA_lst.nops()>0) {
268 ff <<
"Dimension NA;" << endl;
270 for(
auto ix : CA_lst) {
271 auto i = ex_to<Index>(ix);
273 st[i.name.get_name()] = i;
279 ff <<
".global" << endl;
283 bool islst = is_a<lst>(expr);
285 if(islst) expr_lst = ex_to<lst>(expr);
286 else expr_lst.append(expr);
287 auto total = expr_lst.nops();
294 map<ex,int,ex_is_less> e2i_map;
295 for(
auto it : expr_lst) {
296 ff <<
".store" << endl;
298 item = MapFunction([](
const ex & e, MapFunction &self)->ex{
299 if(!e.has(TR(
w)))
return e;
300 else if(e.match(TR(
w))) {
301 ex nd = e.op(0).normal().numer_denom();
302 return TR(nd.op(0))/nd.op(1);
303 }
else return e.map(self);
307 item = collect_common_factors(item);
309 auto ckey = item.op(0);
310 if(!is_a<numeric>(ckey)) {
312 if(e2i_map.find(ckey)==e2i_map.end()) {
316 }
else ckid = e2i_map[ckey];
317 string gCF =
"gCF" + to_string(ckid);
318 st[gCF] = item.op(0);
319 item = Symbol(gCF) * item.op(1);
321 item = ckey * item.op(1);
328 for(
int i=0; i<cv_lst.nops(); i++) {
329 auto it = cv_lst.op(i);
332 color_vec.push_back(vv);
333 item += cc * Symbol(
"[cl"+to_string(i)+
"]");
337 for(
int i=0; i<color_vec.size(); i++) {
338 coss <<
" [cl" << i <<
"]";
339 ff <<
"L [cl" << i <<
"]=" << color_vec[i] <<
";" << endl;
340 ff <<
".sort" << endl;
341 ff <<
"#call SUNTrace" << endl;
342 ff <<
".sort" << endl;
344 ff <<
"id NF^reX?=3^reX;" << endl;
345 ff <<
"id I2R=1/2;" << endl;
346 ff <<
"id NA^reX?=8^reX;" << endl;
347 ff <<
".sort" << endl;
355 find(item,TR(
w),trs);
363 ff <<
"L [o]=" << item <<
";" << endl;
364 ff <<
".sort" << endl;
365 ff <<
"#do i = 1,1" << endl;
366 ff <<
"id once HF(reX?)=reX;" << endl;
368 ff <<
"if(count(HF,1)>0) redefine i \"0\";" << endl;
369 ff <<
".sort" << endl;
370 ff <<
"#enddo" << endl;
371 ff <<
"contract 0;" << endl;
372 ff <<
".sort" << endl;
373 ff << idstr <<
".sort" << endl;
375 for(
int gl=1; gl<=tr.glmax; gl++) {
377 ff <<
".sort" << endl;
379 else ff <<
"Dimension d;" << endl;
380 ff <<
"Indices [g5_i1], [g5_i2], [g5_i3], [g5_i4];" << endl;
381 ff <<
"id g_(" << gl <<
",5_) = e_([g5_i1],[g5_i2],[g5_i3],[g5_i4])*g_(" << gl <<
",[g5_i1],[g5_i2],[g5_i3],[g5_i4])/24;" << endl;
382 ff <<
"sum [g5_i1],[g5_i2],[g5_i3],[g5_i4];" << endl;
383 ff <<
".sort" << endl;
386 else ff <<
"tracen " << gl <<
";" << endl;
387 ff <<
".sort" << endl;
388 ff << idstr <<
".sort" << endl;
395 tr2v[tr] = Symbol(
"[tr"+to_string(trn)+
"]");
396 auto tri = mapGamma(gid)(tr.op(0));
397 trvec.push_back(tri);
400 item = item.subs(tr2v);
401 for(
int i=0; i<trvec.size(); i++) {
402 ff <<
"L [tr" << i <<
"]=" << trvec[i] <<
";" << endl;
404 ff <<
".sort" << endl;
406 else ff <<
"Dimension d;" << endl;
407 ff <<
"Indices [g5_i1], [g5_i2], [g5_i3], [g5_i4];" << endl;
408 ff <<
"id g_(" << gid <<
",5_) = e_([g5_i1],[g5_i2],[g5_i3],[g5_i4])*g_(" << gid <<
",[g5_i1],[g5_i2],[g5_i3],[g5_i4])/24;" << endl;
409 ff <<
"sum [g5_i1],[g5_i2],[g5_i3],[g5_i4];" << endl;
412 else ff <<
"tracen " << gid <<
";" << endl;
413 ff <<
".sort" << endl;
414 ff << idstr <<
".sort" << endl;
416 ff <<
"L [o]=" << item <<
";" << endl;
417 ff <<
".sort" << endl;
418 ff <<
"#do i = 1,1" << endl;
419 ff <<
"id once HF(reX?)=reX;" << endl;
421 ff <<
"if(count(HF,1)>0) redefine i \"0\";" << endl;
422 ff <<
".sort" << endl;
423 ff <<
"#enddo" << endl;
424 ff << idstr <<
".sort" << endl;
431 tr2v[tr] = Symbol(
"trcN"+to_string(trn));
432 auto tri = mapGamma(gid)(tr.op(0));
433 trvec.push_back(tri);
436 item = item.subs(tr2v);
438 if(color_vec.size()>0) ff <<
"drop" << coss.str() <<
";" << endl;
439 ff <<
"L [o]=" << item <<
";" << endl;
440 ff <<
".sort" << endl;
441 ff <<
"#do i = 1,1" << endl;
442 ff <<
"id once HF(reX?)=reX;" << endl;
444 ff <<
"if(count(HF,1)>0) redefine i \"0\";" << endl;
445 ff <<
".sort" << endl;
446 ff <<
"#enddo" << endl;
449 for(
int i=0; i<trvec.size(); i++) {
450 ff <<
"L [tr" << i <<
"]=" << trvec[i] <<
";" << endl;
452 ff <<
".sort" << endl;
454 else ff <<
"Dimension d;" << endl;
455 ff <<
"Indices [g5_i1], [g5_i2], [g5_i3], [g5_i4];" << endl;
456 ff <<
"id g_(" << gid <<
",5_) = e_([g5_i1],[g5_i2],[g5_i3],[g5_i4])*g_(" << gid <<
",[g5_i1],[g5_i2],[g5_i3],[g5_i4])/24;" << endl;
457 ff <<
"sum [g5_i1],[g5_i2],[g5_i3],[g5_i4];" << endl;
460 else ff <<
"tracen " << gid <<
";" << endl;
461 ff <<
".sort" << endl;
462 ff << idstr <<
".sort" << endl;
463 ff <<
"drop [tr" << i <<
"];" << endl;
464 ff <<
"id trcN" << i <<
"=[tr" << i <<
"];" << endl;
465 ff <<
".sort" << endl;
466 ff << idstr <<
".sort" << endl;
470 throw Error(
"runform: unsupported form_trace_mode = " + to_string(
form_trace_mode));
473 ff <<
"contract 0;" << endl;
474 ff <<
".sort" << endl;
475 ff << idstr <<
".sort" << endl;
479 cout <<
PRE <<
"\\--Form Script @ " << c <<
" / " << total << flush;
481 cout <<
"--------------------------------------" << endl;
482 cout <<
"Form Script @ " << c <<
" / " << total << endl;
483 cout <<
"--------------------------------------" << endl;
484 cout << ss.str() << endl;
487 auto script = ss.str();
491 auto otmp = fprc.Execute(script);
493 cout <<
"--------------------------------------" << endl;
494 cout <<
"Form Output @ " << c <<
" / " << total << endl;
495 cout <<
"--------------------------------------" << endl;
496 cout << otmp << endl;
503 if(c<total) ostr +=
",";
505 }
catch(Error& err) {
510 if(verb==1) cout << endl;
518 for(
auto v : vec_lst) {
519 string pat(ex_to<Vector>(v).name.get_name());
520 string from = pat+
"(";
521 string to =
"d_("+pat+
",";
546 fp.FTable.insert({{
"SP", 2}, reader_func(SP_reader)});
547 fp.FTable.insert({{
"SP", 1}, reader_func(SP_reader)});
549 for(
int i=3; i<10; i++) fp.FTable.insert({{
"LC", i}, reader_func(LC_reader)});
550 fp.FTable.insert({{
"GMat", 3}, reader_func(GMat_reader)});
551 for(
int i=1; i<30; i++) fp.FTable.insert({{
"T", i}, reader_func(SUNT_reader)});
552 for(
int i=1; i<30; i++) fp.FTable.insert({{
"TTRX", i}, reader_func(TTR_reader)});
553 for(
int i=1; i<30; i++) fp.FTable.insert({{
"DG", i}, reader_func(DGamma_reader)});
554 fp.FTable.insert({{
"GI", 0}, reader_func(gi_reader)});
555 fp.FTable.insert({{
"GI", 1}, reader_func(gi_reader)});
556 ex ret = fp.Read(ostr);
557 if(!islst) ret = ret.op(0);
567 ex
form(
const ex & iexpr,
int verb) {
569 if(!is_a<lst>(expr)) {
570 static MapFunction mf([](
const ex & e, MapFunction & self)->ex{
571 if(!e.has(TR(w)))
return e;
572 else if(e.match(TR(w))) {
573 if(is_a<mul>(e.op(0))) {
575 for(auto const & ei : e.op(0)) {
576 if(DGamma::has(ei)) v *= ei;
581 }
else return e.map(self);
585 if(form_expand_mode==form_expand_none || is_a<lst>(expr))
return runform(expr, verb);
586 else if(form_expand_mode==form_expand_tr) {
587 auto cv_lst =
collect_lst(expr.subs(SP_map), TR(w));
589 for(
auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
590 lst out_lst = ex_to<lst>(runform(to_lst, verb));
593 for(
int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
595 return ret.subs(SP_map);
596 }
else if(form_expand_mode==form_expand_ci) {
597 auto cv_lst =
collect_lst(expr.subs(SP_map), [](
const ex & e)->bool {
598 return e.has(TR(w)) || Index::hasc(e);
601 for(
auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
602 lst out_lst = ex_to<lst>(runform(to_lst, verb));
605 for(
int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
607 return ret.subs(SP_map);
608 }
else if(form_expand_mode==form_expand_li) {
609 auto cv_lst =
collect_lst(expr.subs(SP_map), [](
const ex & e)->bool {
610 return e.has(TR(w)) || Index::hasv(e);
613 for(
auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
614 lst out_lst = ex_to<lst>(runform(to_lst, verb));
617 for(
int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
619 return ret.subs(SP_map);
620 }
else if(form_expand_mode==form_expand_all) {
621 auto cv_lst =
collect_lst(expr.subs(SP_map), [](
const ex & e)->bool {
622 return e.has(TR(w)) || SUNT::has(e) || SUNF::has(e) || SUNF4::has(e) || Index::has(e) || DGamma::has(e);
625 for(
auto cv : cv_lst) to_lst.append(cv.op(1));
626 lst out_lst = ex_to<lst>(runform(to_lst, verb));
629 for(
int i=0; i<cv_lst.nops(); i++) ret += cv_lst.op(i).op(0) * out_lst.op(i);
631 return ret.subs(SP_map);
632 }
else throw Error(
"form: unsupported form_expand_mode = " + to_string(form_expand_mode));
642 if(expr.has(GMat(w1,w2,w3)))
throw Error(
"charge_conjugate: GMat found.");
644 if(!DGamma::has(expr) && !AsGamma::has(expr))
return expr;
645 if(is_a<DGamma>(expr)) {
646 DGamma g = ex_to<DGamma>(expr);
647 if(is_a<Vector>(g.pi) || is_a<Index>(g.pi))
return -expr;
648 else if(is_zero(g.pi-1) || is_zero(g.pi-5))
return expr;
649 }
else if(is_a<AsGamma>(expr)) {
650 AsGamma ag = ex_to<AsGamma>(expr);
653 for(
int i=ag.pis.size()-1; i>=0; --i) {
655 if(is_a<Vector>(pi) || is_a<Index>(pi)) sign *= -1;
656 else if(is_zero(pi-1) || is_zero(pi-5)) sign *= 1;
657 else throw Error(
"charge_conjugate: only p/i/1/5 supported.");
660 return sign * AsGamma::from(pis_lst);
663 if(is_a<add>(expr)) {
667 }
else if(is_a<mul>(expr)) {
671 }
else if(is_a<ncmul>(expr)) {
677 cout << expr << endl;
678 throw Error(
"charge_conjugate: unexpected region.");
688 if(expr.has(GMat(w1,w2,w3)))
throw Error(
"charge_conjugate: GMat found.");
689 if(AsGamma::has(expr))
throw Error(
"gamma_transpose: Not supported for AsGamma.");
691 if(!DGamma::has(expr))
return expr;
692 if(is_a<DGamma>(expr)) {
693 if(expr.is_equal(DGamma::C))
return (-1)*DGamma::C;
694 DGamma g = ex_to<DGamma>(expr);
699 if(is_a<add>(expr)) {
703 }
else if(is_a<mul>(expr)) {
707 }
else if(is_a<ncmul>(expr)) {
713 cout <<
"gamma_transpose input:" << expr << endl;
714 throw Error(
"gamma_transpose: unexpected region.");
static bool has(const ex &e)
static bool has(const ex &e)
static bool hasc(const ex &e)
static GiNaC::exmap Dimension
static bool has(const ex &e)
lst CoPat(const ex &e, std::function< bool(const ex &)> f)
ex GAS(const ex &expr, unsigned rl)
function similar to GAD/GSD in FeynClac
const int form_trace_each_all
ex gamma_transpose(const ex &expr)
make the transpose operaton M --> M^T
ex charge_conjugate(const ex &expr)
make the charge conjugate operaton, M -> C^{-1} . M^T . C w.r.t. a GMat object
lst collect_lst(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica, reture the lst { {c1,v1}, {c2,v2}, ... }
void string_replace_all(string &str, const string &from, const string &to)
const int form_trace_each_each
const int form_trace_auto
ex LC(ex pi1, ex pi2, ex pi3, ex pi4)
function similar to LCD in FeynCalc
ex form(const ex &iexpr, int verb)
evalulate expr in form program, see also the form_trace_mode and form_expand_mode
ex SP(const ex &a, bool use_map=false)
int ex2int(ex num)
ex to integer