HepLib
Loading...
Searching...
No Matches
Form.cpp
Go to the documentation of this file.
1
6#include "HEP.h"
7#include "ginac/parse_context.h"
8
9namespace HepLib {
10
11 //-----------------------------------------------------------
12 // Extend Parser for form, copied from ginac/parser of GiNaC
13 //-----------------------------------------------------------
14 namespace {
15
16 alignas(2) static ex SP_reader(const exvector& ev) {
17 if(ev.size()<2) return SP(ev[0]);
18 else return SP(ev[0], ev[1]);
19 }
20
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]);
23 else return Eps(ev);
24 }
25
26 alignas(2) static ex SUNT_reader(const exvector& ev) {
27 int n = ev.size();
28 if(n==3) return SUNT(ev[0], ev[1], ev[2]);
29 else if(n>3) {
30 lst as;
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.");
34 }
35
36 alignas(2) static ex TTR_reader(const exvector& ev) {
37 if(ev.size()==1) return TTR(ev[0]);
38 lst as;
39 for(auto item : ev) as.append(item);
40 return TTR(as);
41 }
42
43 alignas(2) static ex DGamma_reader(const exvector& ev) {
44 int n = ev.size();
45 if(n==1) return GAS(ev[0]);
46 ex ret = 1;
47 int rl = ex2int(ev[0]);
48 for(int i=1; i<n; i++) ret = ret * GAS(ev[i],rl);
49 return ret;
50 }
51
52 alignas(2) static ex gi_reader(const exvector& ev) {
53 int n = ev.size();
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.");
57 }
58
59 alignas(2) static ex GMat_reader(const exvector& ev) {
60 return GMat(ev[0],ev[1],ev[2]);
61 }
62 }
63
64 //-----------------------------------------------------------
65 // Run Form Program
66 //-----------------------------------------------------------
67 namespace {
68
69 struct mapGamma : public map_function {
70 public:
71 ex operator()(const ex &e) {
72 if(is_a<DGamma>(e)) return DGamma(ex_to<DGamma>(e), gline);
73 else if(!DGamma::has(e)) return e;
74 else return e.map(*this);
75 }
76 mapGamma(int _gline) : gline(_gline) { };
77 private:
78 unsigned gline = 0;
79 };
80
81 struct mapTR : public map_function {
82 public:
83 ex operator()(const ex &e) {
84 if (e.match(TR(w))) {
85 ex trs = e.op(0);
86 gline++;
87 trs = mapGamma(gline)(trs);
88 trs = DGamma(1, gline) * trs;
89 if(glmax<gline) {
90 glmax = gline;
91 if(glmax>128) throw Error("too large index with glmax>128.");
92 }
93 return TR(trs);
94 } else if(is_a<add>(e)) {
95 ex res = 0;
96 unsigned gl = gline;
97 unsigned glmax = gline;
98 for(auto item : e) {
99 gline = gl;
100 res += (*this)(item);
101 if(glmax<gline) glmax = gline;
102 }
103 gline = glmax;
104 return res;
105 } else if(is_a<power>(e)) {
106 auto ee = e.op(0);
107 auto nn = e.op(1);
108 if(!ee.has(TR(w))) return e;
109 if(!nn.info(info_flags::posint)) return e; // no need to handle negative powers
110 ex res = 1;
111 int pn = ex_to<numeric>(nn).to_int();
112 for(int i=0; i<pn; i++) {
113 res *= (*this)(ee);
114 }
115 return res;
116 } else return e.map(*this);
117 }
118 unsigned glmax = 0;
119 private:
120 unsigned gline = 0;
121 };
122
123 string init_script = R"EOF(
124CFunction pow,sqrt,gamma,HF,GMat,WF;
125Tensor TTR(cyclic), f(antisymmetric), T, f4, colTp;
126Symbols reX,I2R,NF,NA,d;
127AutoDeclare Symbols gCF, trcN;
128Dimension NA;
129AutoDeclare Index colA;
130Dimension NF;
131AutoDeclare Index colF;
132
133#procedure SUNTrace
134Dimension NF;
135
136repeat;
137 id,once,TTR(?a) = T(?a,colF1,colF1);
138 sum colF1;
139 repeat;
140 id,once,T(colA1?,colA2?,?a,colF1?,colF2?) = T(colA1,colF1,colF3)*T(colA2,?a,colF3,colF2);
141 sum colF3;
142 id,once,f4(colA1?,colA2?,colA3?,colA4?) = f(colA1,colA2,colA5) * f(colA5,colA3,colA4);
143 sum colA5;
144 endrepeat;
145endrepeat;
146
147#do colXi = 1,1
148 if ( count(f,1) || match(T(colA1?,colF1?,colF2?)*T(colA1?,colF3?,colF4?)) ) redefine colXi "0";
149 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);
150 sum colF1,colF2,colF3;
151 id T(colA1?,colF1?,colF2?)*T(colA1?,colF3?,colF4?) = colTp(colF1,colF2,colF3,colF4);
152 #do colXj = 1,1
153 if ( count(colTp,1) ) redefine colXj "0";
154 .sort
155 id,once,colTp(colF1?,colF2?,colF3?,colF4?) = I2R*(d_(colF1,colF4)*d_(colF2,colF3)-d_(colF1,colF2)*d_(colF3,colF4)/NF);
156 #enddo
157#enddo
158
159repeat;
160 id T(colA1?,?a,colF1?,colF2?)*T(colA2?,?b,colF2?,colF3?) = T(colA1,?a,colA2,?b,colF1,colF3);
161endrepeat;
162id T(?a,colF1?,colF1?) = TTR(?a);
163id TTR(colA1?) = 0;
164id TTR(colA1?,colA2?) = I2R*d_(colA1,colA2);
165.sort
166
167Local I = i_;
168Local Pi = pi_;
169.sort
170
171#endprocedure
172.global
173 )EOF";
174 }
175
176 //-----------------------------------------------------------
177 // form
178 //-----------------------------------------------------------
179 namespace {
180 ex runform(const ex &expr_in, int verb) {
181 static map<pid_t, Form> form_map;
182 auto pid = getpid();
183 if((form_map.find(pid)==form_map.end())) { // init section
184 ostringstream ss;
185 ss << init_script << endl;
186 form_map[pid].Init("form");
187 form_map[pid].Execute(ss.str());
188 }
189 Form &fprc = form_map[pid];
190
191 ex expr = expr_in.subs(SP_map);
192 ex all_expr = expr;
193 stringstream sss;
194 FormFormat ids(sss);
195 for(auto kv : SP_map) {
196 ids << "id " << kv.first << "=" << kv.second << ";" << endl;
197 all_expr += iWF(kv.first) * iWF(kv.second);
198 }
199 string idstr = sss.str();
200
201 lst vec_lst, VD_lst, CF_lst, CA_lst, sym_lst;
202 for(const_preorder_iterator i = all_expr.preorder_begin(); i != all_expr.preorder_end(); ++i) {
203 if(is_a<Vector>(*i)) vec_lst.append(*i);
204 else if(is_a<Index>(*i)) {
205 if(ex_to<Index>(*i).dim==d) VD_lst.append(*i);
206 else if(ex_to<Index>(*i).dim==NF) CF_lst.append(*i);
207 else if(ex_to<Index>(*i).dim==NA) CA_lst.append(*i);
208 } else if(is_a<symbol>(*i)) sym_lst.append(*i);
209 else if(is_a<GiNaC::function>(*i)) {
210 static vector<string> fun_vec = {
211 "iWF", "WF", "TR", "sin", "cos", "HF", "TTR", "GMat"
212 };
213 auto func = ex_to<GiNaC::function>(*i).get_name();
214 bool ok = false;
215 for(auto fi : fun_vec) {
216 if(fi==func) { ok = true; break; }
217 }
218 if(!ok) {
219 cout << (*i) << endl;
220 throw Error("runform: Functions not defined in FORM: "+func);
221 }
222 }
223 }
224 vec_lst.sort(); vec_lst.unique();
225 VD_lst.sort(); VD_lst.unique();
226 CF_lst.sort(); CF_lst.unique();
227 CA_lst.sort(); CA_lst.unique();
228 sym_lst.append(d);
229 sym_lst.sort(); sym_lst.unique();
230
231 stringstream ss;
232 FormFormat ff(ss);
233 ff << ".store" << endl;
234 symtab st;
235 if(sym_lst.nops()>0) {
236 ff << "Symbols";
237 for(auto sx : sym_lst) {
238 auto s = ex_to<symbol>(sx);
239 ff << " " << s.get_name();
240 st[s.get_name()] = sx;
241 }
242 ff << ";" << endl;
243 }
244 if(vec_lst.nops()>0) {
245 ff << "Vectors";
246 for(auto vx : vec_lst) {
247 auto v = ex_to<Vector>(vx);
248 ff << " " << v;
249 st[v.name.get_name()] = v;
250 for(auto vvx : vec_lst) {
251 auto vv = ex_to<Vector>(vvx);
252 st[v.name.get_name()+"__"+vv.name.get_name()] = SP(v,vv);
253 }
254 }
255 ff << ";" << endl;
256 }
257 if(VD_lst.nops()>0) {
258 if(form_using_dim4) ff << "Dimension 4;" << endl;
259 else ff << "Dimension d;" << endl;
260 ff << "Indices";
261 for(auto ix : VD_lst) {
262 auto i = ex_to<Index>(ix);
263 ff << " " << i;
264 st[i.name.get_name()] = i;
265 auto itr = Index::Dimension.find(ix);
266 if(itr!=Index::Dimension.end()) ff << "=" << itr->second;
267 }
268 ff << ";" << endl;
269 }
270 if(CF_lst.nops()>0) {
271 ff << "Dimension NF;" << endl;
272 ff << "Indices";
273 for(auto ix : CF_lst) {
274 auto i = ex_to<Index>(ix);
275 ff << " " << i;
276 st[i.name.get_name()] = i;
277 auto itr = Index::Dimension.find(ix);
278 if(itr!=Index::Dimension.end()) ff << "=" << itr->second;
279 }
280 ff << ";" << endl;
281 }
282 if(CA_lst.nops()>0) {
283 ff << "Dimension NA;" << endl;
284 ff << "Indices";
285 for(auto ix : CA_lst) {
286 auto i = ex_to<Index>(ix);
287 ff << " " << i;
288 st[i.name.get_name()] = i;
289 auto itr = Index::Dimension.find(ix);
290 if(itr!=Index::Dimension.end()) ff << "=" << itr->second;
291 }
292 ff << ";" << endl;
293 }
294 ff << ".global" << endl;
295 ff << endl;
296
297 // trace and contract
298 bool islst = is_a<lst>(expr);
299 lst expr_lst;
300 if(islst) expr_lst = ex_to<lst>(expr);
301 else expr_lst.append(expr);
302 auto total = expr_lst.nops();
303
304 string ostr;
305 int gid = 1;
306 ostr = "{";
307 int c = 1;
308 int kid = 0;
309 map<ex,int,ex_is_less> e2i_map;
310 for(auto it : expr_lst) {
311 ff << ".store" << endl;
312 ex item = it;
313 item = MapFunction([](const ex & e, MapFunction &self)->ex{
314 if(!e.has(TR(w))) return e;
315 else if(e.match(TR(w))) {
316 ex nd = e.op(0).normal().numer_denom();
317 return TR(nd.op(0))/nd.op(1);
318 } else return e.map(self);
319 })(item);
320 item = normal(item); // TODO: normal
321 // pull out global common factor
322 item = collect_common_factors(item);
323 item = CoPat(item,[](const ex &e)->bool{return Index::has(e) || DGamma::has(e) || Eps::has(e);});
324 auto ckey = item.op(0);
325 if(!is_a<numeric>(ckey)) {
326 int ckid;
327 if(e2i_map.find(ckey)==e2i_map.end()) {
328 kid++;
329 e2i_map[ckey] = kid;
330 ckid = kid;
331 } else ckid = e2i_map[ckey];
332 string gCF = "gCF" + to_string(ckid);
333 st[gCF] = item.op(0);
334 item = Symbol(gCF) * item.op(1);
335 } else {
336 item = ckey * item.op(1);
337 }
338
339 // pull out color factor
340 auto cv_lst = collect_lst(item, [](const ex &e)->bool{return Index::hasc(e);});
341 item=0;
342 exvector color_vec;
343 for(int i=0; i<cv_lst.nops(); i++) {
344 auto it = cv_lst.op(i);
345 auto cc = it.op(0);
346 auto vv = it.op(1);
347 color_vec.push_back(vv);
348 item += cc * Symbol("[cl"+to_string(i)+"]");
349 }
350
351 ostringstream coss;
352 for(int i=0; i<color_vec.size(); i++) {
353 coss << " [cl" << i << "]";
354 ff << "L [cl" << i << "]=" << color_vec[i] << ";" << endl;
355 ff << ".sort" << endl;
356 ff << "#call SUNTrace" << endl;
357 ff << ".sort" << endl;
358 if(form_using_su3) {
359 ff << "id NF^reX?=3^reX;" << endl;
360 ff << "id I2R=1/2;" << endl;
361 ff << "id NA^reX?=8^reX;" << endl;
362 ff << ".sort" << endl;
363 }
364 }
365 ff << endl;
366
367 // methods to handle TR objects
368 auto tr_mode = form_trace_mode;
369 exset trs;
370 find(item,TR(w),trs);
371 if(form_trace_mode==form_trace_auto) {
372 if(trs.size()<2) tr_mode = form_trace_all;
373 else tr_mode = form_trace_each_each;
374 }
375 if(tr_mode==form_trace_all) {
376 mapTR tr;
377 item = tr(item);
378 ff << "L [o]=" << item << ";" << endl;
379 ff << ".sort" << endl;
380 ff << "#do i = 1,1" << endl;
381 ff << "id once HF(reX?)=reX;" << endl;
382 ff << idstr;
383 ff << "if(count(HF,1)>0) redefine i \"0\";" << endl;
384 ff << ".sort" << endl;
385 ff << "#enddo" << endl;
386 ff << "contract 0;" << endl;
387 ff << ".sort" << endl;
388 ff << idstr << ".sort" << endl;
389 ff << endl;
390 for(int gl=1; gl<=tr.glmax; gl++) {
391 if(form_using_gamma5) {
392 ff << ".sort" << endl;
393 if(form_using_dim4) ff << "Dimension 4;" << endl;
394 else ff << "Dimension d;" << endl;
395 ff << "Indices [g5_i1], [g5_i2], [g5_i3], [g5_i4];" << endl;
396 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;
397 ff << "sum [g5_i1],[g5_i2],[g5_i3],[g5_i4];" << endl;
398 ff << ".sort" << endl;
399 }
400 if(form_using_dim4) ff << "trace4 " << gl << ";" << endl;
401 else ff << "tracen " << gl << ";" << endl;
402 ff << ".sort" << endl;
403 ff << idstr << ".sort" << endl;
404 }
405 } else if(tr_mode==form_trace_each_all) {
406 exmap tr2v;
407 int trn=0;
408 exvector trvec;
409 for(auto tr : trs) {
410 tr2v[tr] = Symbol("[tr"+to_string(trn)+"]");
411 auto tri = mapGamma(gid)(tr.op(0));
412 trvec.push_back(tri);
413 trn++;
414 }
415 item = item.subs(tr2v);
416 for(int i=0; i<trvec.size(); i++) {
417 ff << "L [tr" << i << "]=" << trvec[i] << ";" << endl;
418 if(form_using_gamma5) {
419 ff << ".sort" << endl;
420 if(form_using_dim4) ff << "Dimension 4;" << endl;
421 else ff << "Dimension d;" << endl;
422 ff << "Indices [g5_i1], [g5_i2], [g5_i3], [g5_i4];" << endl;
423 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;
424 ff << "sum [g5_i1],[g5_i2],[g5_i3],[g5_i4];" << endl;
425 }
426 if(form_using_dim4) ff << "trace4 " << gid << ";" << endl;
427 else ff << "tracen " << gid << ";" << endl;
428 ff << ".sort" << endl;
429 ff << idstr << ".sort" << endl;
430 }
431 ff << "L [o]=" << item << ";" << endl;
432 ff << ".sort" << endl;
433 ff << "#do i = 1,1" << endl;
434 ff << "id once HF(reX?)=reX;" << endl;
435 ff << idstr;
436 ff << "if(count(HF,1)>0) redefine i \"0\";" << endl;
437 ff << ".sort" << endl;
438 ff << "#enddo" << endl;
439 ff << idstr << ".sort" << endl;
440 ff << endl;
441 } else if(tr_mode==form_trace_each_each) {
442 exmap tr2v;
443 int trn=0;
444 exvector trvec;
445 for(auto tr : trs) {
446 tr2v[tr] = Symbol("trcN"+to_string(trn));
447 auto tri = mapGamma(gid)(tr.op(0));
448 trvec.push_back(tri);
449 trn++;
450 }
451 item = item.subs(tr2v);
452
453 if(color_vec.size()>0) ff << "drop" << coss.str() << ";" << endl;
454 ff << "L [o]=" << item << ";" << endl;
455 ff << ".sort" << endl;
456 ff << "#do i = 1,1" << endl;
457 ff << "id once HF(reX?)=reX;" << endl;
458 ff << idstr;
459 ff << "if(count(HF,1)>0) redefine i \"0\";" << endl;
460 ff << ".sort" << endl;
461 ff << "#enddo" << endl;
462 ff << endl;
463
464 for(int i=0; i<trvec.size(); i++) {
465 ff << "L [tr" << i << "]=" << trvec[i] << ";" << endl;
466 if(form_using_gamma5) {
467 ff << ".sort" << endl;
468 if(form_using_dim4) ff << "Dimension 4;" << endl;
469 else ff << "Dimension d;" << endl;
470 ff << "Indices [g5_i1], [g5_i2], [g5_i3], [g5_i4];" << endl;
471 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;
472 ff << "sum [g5_i1],[g5_i2],[g5_i3],[g5_i4];" << endl;
473 }
474 if(form_using_dim4) ff << "trace4 " << gid << ";" << endl;
475 else ff << "tracen " << gid << ";" << endl;
476 ff << ".sort" << endl;
477 ff << idstr << ".sort" << endl;
478 ff << "drop [tr" << i << "];" << endl;
479 ff << "id trcN" << i << "=[tr" << i << "];" << endl;
480 ff << ".sort" << endl;
481 ff << idstr << ".sort" << endl;
482 ff << endl;
483 }
484 } else {
485 throw Error("runform: unsupported form_trace_mode = " + to_string(form_trace_mode));
486 }
487
488 ff << "contract 0;" << endl;
489 ff << ".sort" << endl;
490 ff << idstr << ".sort" << endl;
491
492 if(verb==1) {
493 cout << "\r \r";
494 cout << PRE << "\\--Form Script @ " << c << " / " << total << flush;
495 } else if(verb>1) {
496 cout << "--------------------------------------" << endl;
497 cout << "Form Script @ " << c << " / " << total << endl;
498 cout << "--------------------------------------" << endl;
499 cout << ss.str() << endl;
500 }
501
502 auto script = ss.str();
503 string_replace_all(script, "sin(", "sin_(");
504 string_replace_all(script, "cos(", "cos_(");
505 try {
506 auto otmp = fprc.Execute(script);
507 if(verb>2) {
508 cout << "--------------------------------------" << endl;
509 cout << "Form Output @ " << c << " / " << total << endl;
510 cout << "--------------------------------------" << endl;
511 cout << otmp << endl;
512 }
513
514 ostr += otmp;
515 ss.clear();
516 ss.str("");
517
518 if(c<total) ostr += ",";
519 c++;
520 } catch(Error& err) {
521 form_map.erase(pid);
522 throw;
523 }
524 }
525 if(verb==1) cout << endl;
526 ostr += "}";
527
528 string_replace_all(ostr, "[", "(");
529 string_replace_all(ostr, "]", ")");
530 string_replace_all(ostr, "\\\n", "");
531 string_replace_all(ostr, "\n", "");
532 string_replace_all(ostr, " ","");
533 for(auto v : vec_lst) {
534 string pat(ex_to<Vector>(v).name.get_name());
535 string from = pat+"(";
536 string to = "d_("+pat+",";
537 string_replace_all(ostr, from, to);
538 from = pat+".";
539 to = pat+"__";
540 string_replace_all(ostr, from, to);
541 }
542
543 string_replace_all(ostr, "d_(", "SP(");
544 string_replace_all(ostr, "e_(", "LC(");
545 string_replace_all(ostr, "sin_(", "sin(");
546 string_replace_all(ostr, "cos_(", "cos(");
547 string_replace_all(ostr, "g_(", "DG(");
548 string_replace_all(ostr, ",5_", ",5");
549 string_replace_all(ostr, ",6_", ",6");
550 string_replace_all(ostr, ",7_", ",7");
551 string_replace_all(ostr, "gi_(", "GI(");
552 string_replace_all(ostr, "TTR(", "TTRX(");
553
554 st["I2R"] = TF;
555 st["NA"] = NA;
556 st["NF"] = NF;
557 st["I"] = I;
558 st["i_"] = I;
559 st["pi_"] = Pi;
560
561 Parser fp(st);
562 fp.FTable.insert({{"SP", 2}, reader_func(SP_reader)});
563 fp.FTable.insert({{"SP", 1}, reader_func(SP_reader)});
564 //fp.FTable.insert({{"LC", 4}, reader_func(LC_reader)});
565 for(int i=3; i<10; i++) fp.FTable.insert({{"LC", i}, reader_func(LC_reader)});
566 fp.FTable.insert({{"GMat", 3}, reader_func(GMat_reader)});
567 for(int i=1; i<30; i++) fp.FTable.insert({{"T", i}, reader_func(SUNT_reader)});
568 for(int i=1; i<30; i++) fp.FTable.insert({{"TTRX", i}, reader_func(TTR_reader)});
569 for(int i=1; i<30; i++) fp.FTable.insert({{"DG", i}, reader_func(DGamma_reader)});
570 fp.FTable.insert({{"GI", 0}, reader_func(gi_reader)});
571 fp.FTable.insert({{"GI", 1}, reader_func(gi_reader)});
572 ex ret = fp.Read(ostr);
573 if(!islst) ret = ret.op(0);
574 return ret;
575 }}
576
583 ex form(const ex & iexpr, int verb) {
584 ex expr = iexpr;
585 if(!is_a<lst>(expr)) {
586 static MapFunction mf([](const ex & e, MapFunction & self)->ex{
587 if(!e.has(TR(w))) return e;
588 else if(e.match(TR(w))) {
589 if(is_a<mul>(e.op(0))) {
590 ex c = 1, v = 1;
591 for(auto const & ei : e.op(0)) {
592 if(DGamma::has(ei)) v *= ei;
593 else c *= ei;
594 }
595 return c*TR(v);
596 } else return e;
597 } else return e.map(self);
598 });
599 expr = mf(expr);
600 }
601 if(form_expand_mode==form_expand_none || is_a<lst>(expr)) return runform(expr, verb);
602 else if(form_expand_mode==form_expand_tr) {
603 auto cv_lst = collect_lst(expr.subs(SP_map), TR(w));
604 lst to_lst;
605 for(auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
606 lst out_lst = ex_to<lst>(runform(to_lst, verb));
607
608 ex ret = 0;
609 for(int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
610
611 return ret.subs(SP_map);
612 } else if(form_expand_mode==form_expand_ci) {
613 auto cv_lst = collect_lst(expr.subs(SP_map), [](const ex & e)->bool {
614 return e.has(TR(w)) || Index::hasc(e);
615 });
616 lst to_lst;
617 for(auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
618 lst out_lst = ex_to<lst>(runform(to_lst, verb));
619
620 ex ret = 0;
621 for(int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
622
623 return ret.subs(SP_map);
624 } else if(form_expand_mode==form_expand_li) {
625 auto cv_lst = collect_lst(expr.subs(SP_map), [](const ex & e)->bool {
626 return e.has(TR(w)) || Index::hasv(e);
627 });
628 lst to_lst;
629 for(auto cv : cv_lst) to_lst.append(cv.op(0)*cv.op(1));
630 lst out_lst = ex_to<lst>(runform(to_lst, verb));
631
632 ex ret = 0;
633 for(int i=0; i<cv_lst.nops(); i++) ret += out_lst.op(i);
634
635 return ret.subs(SP_map);
636 } else if(form_expand_mode==form_expand_all) {
637 auto cv_lst = collect_lst(expr.subs(SP_map), [](const ex & e)->bool {
638 return e.has(TR(w)) || SUNT::has(e) || SUNF::has(e) || SUNF4::has(e) || Index::has(e) || DGamma::has(e);
639 });
640 lst to_lst;
641 for(auto cv : cv_lst) to_lst.append(cv.op(1));
642 lst out_lst = ex_to<lst>(runform(to_lst, verb));
643
644 ex ret = 0;
645 for(int i=0; i<cv_lst.nops(); i++) ret += cv_lst.op(i).op(0) * out_lst.op(i);
646
647 return ret.subs(SP_map);
648 } else throw Error("form: unsupported form_expand_mode = " + to_string(form_expand_mode));
649 return 0;
650 }
651
657 ex charge_conjugate(const ex & expr) {
658 if(expr.has(GMat(w1,w2,w3))) throw Error("charge_conjugate: GMat found.");
659
660 if(!DGamma::has(expr) && !AsGamma::has(expr)) return expr;
661 if(is_a<DGamma>(expr)) {
662 DGamma g = ex_to<DGamma>(expr);
663 if(is_a<Vector>(g.pi) || is_a<Index>(g.pi)) return -expr;
664 else if(is_zero(g.pi-1) || is_zero(g.pi-5)) return expr;
665 } else if(is_a<AsGamma>(expr)) {
666 AsGamma ag = ex_to<AsGamma>(expr);
667 ex sign = 1;
668 lst pis_lst;
669 for(int i=ag.pis.size()-1; i>=0; --i) {
670 ex pi = ag.pis[i];
671 if(is_a<Vector>(pi) || is_a<Index>(pi)) sign *= -1;
672 else if(is_zero(pi-1) || is_zero(pi-5)) sign *= 1;
673 else throw Error("charge_conjugate: only p/i/1/5 supported.");
674 pis_lst.append(pi);
675 }
676 return sign * AsGamma::from(pis_lst);
677 }
678
679 if(is_a<add>(expr)) {
680 ex ret = 0;
681 for(auto item : expr) ret += charge_conjugate(item);
682 return ret;
683 } else if(is_a<mul>(expr)) {
684 ex ret = 1;
685 for(auto item : expr) ret *= charge_conjugate(item);
686 return ret;
687 } else if(is_a<ncmul>(expr)) {
688 ex ret = 1;
689 int n = expr.nops();
690 for(int i=n-1; i>-1; i--) ret *= charge_conjugate(expr.op(i));
691 return ret;
692 }
693 cout << expr << endl;
694 throw Error("charge_conjugate: unexpected region.");
695 return 0;
696 }
697
703 ex gamma_transpose(const ex & expr) {
704 if(expr.has(GMat(w1,w2,w3))) throw Error("charge_conjugate: GMat found.");
705 if(AsGamma::has(expr)) throw Error("gamma_transpose: Not supported for AsGamma.");
706
707 if(!DGamma::has(expr)) return expr;
708 if(is_a<DGamma>(expr)) {
709 if(expr.is_equal(DGamma::C)) return (-1)*DGamma::C;
710 DGamma g = ex_to<DGamma>(expr);
711 g.isTr = !g.isTr;
712 return g;
713 }
714
715 if(is_a<add>(expr)) {
716 ex ret = 0;
717 for(auto item : expr) ret += gamma_transpose(item);
718 return ret;
719 } else if(is_a<mul>(expr)) {
720 ex ret = 1;
721 for(auto item : expr) ret *= gamma_transpose(item);
722 return ret;
723 } else if(is_a<ncmul>(expr)) {
724 ex ret = 1;
725 int n = expr.nops();
726 for(int i=n-1; i>-1; i--) ret *= gamma_transpose(expr.op(i));
727 return ret;
728 }
729 cout << "gamma_transpose input:" << expr << endl;
730 throw Error("gamma_transpose: unexpected region.");
731 return 0;
732 }
733
734}
735
glmax
unsigned glmax
Definition Form.cpp:118
HEP.h
HEP header file.
pn
int pn
Definition Others.cpp:25
HepLib::DGamma::has
static bool has(const ex &e)
HepLib::Eps::has
static bool has(const ex &e)
HepLib::Index::hasc
static bool hasc(const ex &e)
Definition Basic.cpp:234
HepLib::Index::Dimension
static GiNaC::exmap Dimension
Definition HEP.h:109
HepLib::Index::has
static bool has(const ex &e)
HepLib
HepLib namespace.
Definition BASIC.cpp:17
HepLib::form_using_dim4
bool form_using_dim4
Definition Init.cpp:218
HepLib::NA
const Symbol NA
HepLib::CoPat
lst CoPat(const ex &e, std::function< bool(const ex &)> f)
Definition BASIC.h:772
HepLib::NF
const Symbol NF
HepLib::SP_map
exmap SP_map
Definition Init.cpp:188
HepLib::form_trace_all
const int form_trace_all
Definition Init.cpp:205
HepLib::GAS
ex GAS(const ex &expr, unsigned rl)
function similar to GAD/GSD in FeynClac
Definition DGamma.cpp:280
HepLib::form_trace_mode
int form_trace_mode
Definition Init.cpp:208
HepLib::form_using_su3
bool form_using_su3
Definition Init.cpp:217
HepLib::w
ex w
Definition Init.cpp:93
HepLib::form_trace_each_all
const int form_trace_each_all
Definition Init.cpp:206
HepLib::d
const Symbol d
HepLib::TF
const Symbol TF
HepLib::gamma_transpose
ex gamma_transpose(const ex &expr)
make the transpose operaton M --> M^T
Definition Form.cpp:701
HepLib::as
const Symbol as
HepLib::charge_conjugate
ex charge_conjugate(const ex &expr)
make the charge conjugate operaton, M -> C^{-1} . M^T . C w.r.t. a GMat object
Definition Form.cpp:655
HepLib::collect_lst
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}, ... }
Definition BASIC.cpp:1223
HepLib::form_using_gamma5
bool form_using_gamma5
Definition Init.cpp:219
HepLib::string_replace_all
void string_replace_all(string &str, const string &from, const string &to)
Definition Functions.cpp:148
HepLib::form_trace_each_each
const int form_trace_each_each
Definition Init.cpp:207
HepLib::form_trace_auto
const int form_trace_auto
Definition Init.cpp:204
HepLib::SP
ex SP(const ex &a, bool use_map=true)
Definition Pair.cpp:166
HepLib::LC
ex LC(ex pi1, ex pi2, ex pi3, ex pi4)
function similar to LCD in FeynCalc
Definition Eps.cpp:179
HepLib::form
ex form(const ex &iexpr, int verb)
evalulate expr in form program, see also the form_trace_mode and form_expand_mode
Definition Form.cpp:581
HepLib::PRE
string PRE
Definition Init.cpp:146
HepLib::ex2int
int ex2int(ex num)
ex to integer
Definition BASIC.cpp:894
Generated by doxygen 1.9.8