HepLib
Loading...
Searching...
No Matches
Basic.cpp
Go to the documentation of this file.
1
6#include "HEP.h"
7#include "flint/ulong_extras.h"
8
9namespace HepLib {
10
11 DEFAULT_CTOR(Index)
12 IMPLEMENT_HAS(Index)
13 IMPLEMENT_ALL(Index)
14
15 DEFAULT_CTOR(Vector)
16 IMPLEMENT_HAS(Vector)
17 IMPLEMENT_ALL(Vector)
18
19 DEFAULT_CTOR(SUNT)
20 IMPLEMENT_HAS(SUNT)
21 IMPLEMENT_ALL(SUNT)
22
23 DEFAULT_CTOR(SUNF)
24 IMPLEMENT_HAS(SUNF)
25 IMPLEMENT_ALL(SUNF)
26
27 DEFAULT_CTOR(SUNF4)
28 IMPLEMENT_HAS(SUNF4)
29 IMPLEMENT_ALL(SUNF4)
30
31 //-----------------------------------------------------------
32 // FormFormat Output
33 //-----------------------------------------------------------
34 FormFormat::FormFormat(ostream &os, unsigned opt) : print_dflt(os, opt) {}
35 FormFormat::FormFormat() : print_dflt(std::cout) {}
36 GINAC_IMPLEMENT_PRINT_CONTEXT(FormFormat, print_dflt)
37
38 const FormFormat & FormFormat::operator << (const basic & v) const {
39 v.print(*this);
40 return *this;
41 }
42 const FormFormat & FormFormat::operator << (const ex & v) const {
43 v.print(*this);
44 return *this;
45 }
46 const FormFormat & FormFormat::operator << (const lst & v) const {
47 v.print(*this);
48 return *this;
49 }
50 const FormFormat & FormFormat::operator<<(std::ostream& (*v)(std::ostream&)) const {
51 s << v;
52 return *this;
53 }
54
55 void FormFormat::power_print(const power & p, const FormFormat & c, unsigned level) {
56 if(p.op(1)==2 && !DGamma::has(p)) {
57 c << "((" << p.op(0) << ")*(" << p.op(0) << "))";
58 } else {
59 c << "((" << p.op(0) << ")^(" << p.op(1) << "))";
60 }
61 }
62
63 //-----------------------------------------------------------
64 // FCFormat Output
65 //-----------------------------------------------------------
66 FCFormat::FCFormat(ostream &os, unsigned opt) : print_dflt(os, opt) {}
67 FCFormat::FCFormat() : print_dflt(std::cout) {}
68 GINAC_IMPLEMENT_PRINT_CONTEXT(FCFormat, print_dflt)
69
70 const FCFormat & FCFormat::operator << (const basic & v) const {
71 v.print(*this);
72 return *this;
73 }
74 const FCFormat & FCFormat::operator << (const ex & v) const {
75 v.print(*this);
76 return *this;
77 }
78 const FCFormat & FCFormat::operator << (const lst & v) const {
79 v.print(*this);
80 return *this;
81 }
82 const FCFormat & FCFormat::operator<<(std::ostream& (*v)(std::ostream&)) const {
83 s << v;
84 return *this;
85 }
86
87 const FCFormat & FCFormat::operator << (const matrix & mat) const {
88 s << "{";
89 int nr = mat.rows();
90 int nc = mat.cols();
91 for(int r=0; r<nr; r++) {
92 s << "{";
93 for(int c=0; c<nc; c++) {
94 mat(r,c).print(*this);
95 if(c+1!=nc) s << ",";
96 }
97 s << "}";
98 if(r+1!=nr) s << ",";
99 }
100 s << "}";
101 return *this;
102 }
103
104 const FCFormat & FCFormat::operator << (const exvector & e) const {
105 auto i = e.begin();
106 auto vend = e.end();
107 if (i==vend) { s << "{}"; return *this; }
108 s << "{";
109 while (true) {
110 i->print(*this);
111 ++i;
112 if(i==vend) break;
113 s << ",";
114 }
115 s << "}";
116 return *this;
117 }
118
119 const FCFormat & FCFormat::operator << (const exset & e) const {
120 auto i = e.begin();
121 auto send = e.end();
122 if (i==send) { s << "{}"; return *this; }
123 s << "{";
124 while (true) {
125 i->print(*this);
126 ++i;
127 if(i==send) break;
128 s << ",";
129 }
130 s << "}";
131 return *this;
132 }
133
134 const FCFormat & FCFormat::operator << (const exmap & e) const {
135 auto i = e.begin();
136 auto mend = e.end();
137 if (i==mend) { s << "{}"; return *this; }
138 s << "{";
139 while (true) {
140 i->first.print(*this);
141 s << "->";
142 i->second.print(*this);
143 ++i;
144 if(i==mend) break;
145 s << ",";
146 }
147 s << "}";
148 return *this;
149 }
150
151 namespace {
152 class ncmul_hack : public ncmul { // due to printseq is protected
153 public:
154 ncmul_hack(ncmul _nm) : ncmul(_nm){ }
155 void print(const FCFormat & c, unsigned level) {
156 printseq(c, '(', '.', ')', precedence(), level);
157 }
158 };
159 ex mat_conj(const ex & e1, const ex & e2, const ex & e3) {
160 return GMat(e1.conjugate(), e3, e2);
161 }
162 }
163 void FCFormat::ncmul_print(const ncmul & nm, const FCFormat & c, unsigned level) {
164 ncmul_hack(nm).print(c, level);
165 }
166
167 //-----------------------------------------------------------
168 // Index Class
169 //-----------------------------------------------------------
170 //GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Index, basic,print_func<print_context>(&Index::print))
171 GiNaC::registered_class_info & Index::get_class_info_static() { return reg_info; }
173 Index * Index::duplicate() const { Index * bp = new Index(*this); bp->setflag(GiNaC::status_flags::dynallocated); return bp; }
174 void Index::accept(GiNaC::visitor & v) const { if (visitor *p = dynamic_cast<visitor *>(&v)) p->visit(*this); else inherited::accept(v); }
175 const GiNaC::registered_class_info &Index::get_class_info() const { return get_class_info_static(); }
176 GiNaC::registered_class_info &Index::get_class_info() { return get_class_info_static(); }
177 const char *Index::class_name() const { return get_class_info_static().options.get_name(); }
178 //GINAC_IMPLEMENT_REGISTERED_CLASS END
179
180 Index::Index(const string &s, const ex & dimension) : name(s), dim(dimension) { }
181 int Index::compare_same_type(const basic &other) const {
182 if(!is_a<Index>(other)) throw Error("Index::compare_same_type");
183 const Index &o = static_cast<const Index &>(other);
184 auto ret = name.get_name().compare(o.name.get_name());
185 if(ret<0) return -1;
186 else if(ret>0) return 1;
187 else return dim.compare(o.dim);
188 }
189
190 bool Index::is_equal_same_type(const basic & other) const {
191 if(!is_a<Index>(other)) throw Error("Index::is_equal_same_type");
192 const Index &o = static_cast<const Index &>(other);
193 auto ret = name.get_name() == o.name.get_name();
194 if(!ret) return false;
195 return dim.is_equal(o.dim);
196 }
197
198 void Index::print(const print_context &c, unsigned level) const {
199 c.s << name;
200 }
201
203 return Pair(*this, i);
204 }
205
207 return Pair(p, *this);
208 }
209
210 void Index::archive(archive_node & n) const {
211 inherited::archive(n);
212 n.add_string("name", name.get_name());
213 n.add_ex("dim", dim);
214 }
215
216 void Index::read_archive(const archive_node& n) {
217 inherited::read_archive(n);
218 string nstr;
219 n.find_string("name", nstr);
220 name = Symbol(nstr);
221 n.find_ex("dim", dim);
222 }
223
224 ex Index::derivative(const symbol & s) const {
225 return 0;
226 }
227
228 bool Index::has(const ex & e, const ex & DIM) {
229 for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i)
230 if(is_a<Index>(*i) && ex_to<Index>(*i).dim==DIM) return true;
231 return false;
232 }
233
234 bool Index::hasc(const ex & e) {
235 for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i)
236 if(is_a<Index>(*i) && (ex_to<Index>(*i).dim==NA || ex_to<Index>(*i).dim==NF)) return true;
237 return false;
238 }
239
240 bool Index::hasv(const ex & e) {
241 for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i)
242 if(is_a<Index>(*i) && ex_to<Index>(*i).dim==d) return true;
243 return false;
244 }
245
246 //-----------------------------------------------------------
247 // Vector Class
248 //-----------------------------------------------------------
249 //GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Vector, basic,print_func<print_context>(&Vector::print))
250 GiNaC::registered_class_info & Vector::get_class_info_static() { return reg_info; }
252 Vector * Vector::duplicate() const { Vector * bp = new Vector(*this); bp->setflag(GiNaC::status_flags::dynallocated); return bp; }
253 void Vector::accept(GiNaC::visitor & v) const { if (visitor *p = dynamic_cast<visitor *>(&v)) p->visit(*this); else inherited::accept(v); }
254 const GiNaC::registered_class_info &Vector::get_class_info() const { return get_class_info_static(); }
255 GiNaC::registered_class_info &Vector::get_class_info() { return get_class_info_static(); }
256 const char *Vector::class_name() const { return get_class_info_static().options.get_name(); }
257 //GINAC_IMPLEMENT_REGISTERED_CLASS END
258
259 Vector::Vector(const string &s) : name(s) { }
260 int Vector::compare_same_type(const basic &other) const {
261 if(!is_a<Vector>(other)) throw Error("Vector::compare_same_type");
262 const Vector &o = static_cast<const Vector &>(other);
263 auto ret = name.get_name().compare(o.name.get_name());
264 if(ret==0) return 0;
265 else if(ret<0) return -1;
266 else return 1;
267 }
268
269 bool Vector::is_equal_same_type(const basic & other) const {
270 if(!is_a<Vector>(other)) throw Error("Vector::is_equal_same_type");
271 const Vector &o = static_cast<const Vector &>(other);
272 return (name.get_name() == o.name.get_name());
273 }
274
275 void Vector::print(const print_context &c, unsigned level) const {
276 c.s << name;
277 }
278
280 return Pair(*this, p);
281 }
282
284 return Pair(*this, i);
285 }
286
287 void Vector::archive(archive_node & n) const {
288 inherited::archive(n);
289 n.add_string("name", name.get_name());
290 }
291
292 void Vector::read_archive(const archive_node& n) {
293 inherited::read_archive(n);
294 string nstr;
295 unsigned t;
296 n.find_string("name", nstr);
297 name = Symbol(nstr);
298 }
299
300 ex Vector::derivative(const symbol & s) const {
301 return 0;
302 }
303
304 //GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(SUNT, basic,print_func<print_dflt>(&SUNT::print).print_func<FormFormat>(&SUNT::form_print).print_func<FCFormat>(&SUNT::fc_print))
305 GiNaC::registered_class_info & SUNT::get_class_info_static() { return reg_info; }
307 SUNT * SUNT::duplicate() const { SUNT * bp = new SUNT(*this); bp->setflag(GiNaC::status_flags::dynallocated); return bp; }
308 void SUNT::accept(GiNaC::visitor & v) const { if (visitor *p = dynamic_cast<visitor *>(&v)) p->visit(*this); else inherited::accept(v); }
309 const GiNaC::registered_class_info &SUNT::get_class_info() const { return get_class_info_static(); }
310 GiNaC::registered_class_info &SUNT::get_class_info() { return get_class_info_static(); }
311 const char *SUNT::class_name() const { return get_class_info_static().options.get_name(); }
312 //GINAC_IMPLEMENT_REGISTERED_CLASS END
313
314 SUNT::SUNT(ex a, ex i, ex j) : aij{a,i,j} { }
315 int SUNT::compare_same_type(const basic &other) const {
316 if(!is_a<SUNT>(other)) throw Error("SUNT::compare_same_type");
317 const SUNT &o = static_cast<const SUNT &>(other);
318 for(int i=0; i<3; i++) {
319 auto c = aij[i].compare(o.aij[i]);
320 if(c!=0) return c;
321 }
322 return 0;
323 }
324
325 bool SUNT::is_equal_same_type(const basic & other) const {
326 if(!is_a<SUNT>(other)) throw Error("SUNT::is_equal_same_type");
327 const SUNT &o = static_cast<const SUNT &>(other);
328 for(int i=0; i<3; i++) {
329 if(!aij[i].is_equal(o.aij[i])) return false;
330 }
331 return true;
332 }
333
334 void SUNT::form_print(const FormFormat &c, unsigned level) const {
335 if(is_a<lst>(aij[0])) {
336 bool first = true;
337 for(auto item : aij[0]) {
338 if(first) { first=false; c << "T(" << item; }
339 else c << "," << item;
340 }
341 } else c << "T(" << aij[0];
342 c << "," << aij[1] << "," << aij[2] << ")";
343 }
344
345 void SUNT::fc_print(const FCFormat &c, unsigned level) const {
346 c << "SUNTF[" << aij[0] << "," << aij[1] << "," << aij[2] << "]";
347 }
348
349 void SUNT::print(const print_dflt &c, unsigned level) const {
350 c.s << "T(" << aij[0] << "," << aij[1] << "," << aij[2] << ")";
351 }
352
353 size_t SUNT::nops() const { return 3; }
354 ex SUNT::op(size_t i) const {
355 return aij[i];
356 }
357 ex & SUNT::let_op(size_t i) {
358 ensure_if_modifiable();
359 return aij[i];
360 }
361
362 void SUNT::archive(archive_node & n) const {
363 inherited::archive(n);
364 n.add_ex("a", aij[0]);
365 n.add_ex("i", aij[1]);
366 n.add_ex("j", aij[2]);
367 }
368
369 void SUNT::read_archive(const archive_node& n) {
370 inherited::read_archive(n);
371 ex o;
372 n.find_ex("a", o);
373 aij[0] = ex_to<Index>(o);
374 n.find_ex("i", o);
375 aij[1] = ex_to<Index>(o);
376 n.find_ex("j", o);
377 aij[2] = ex_to<Index>(o);
378 }
379
380 ex SUNT::derivative(const symbol & s) const {
381 return 0;
382 }
383
384 ex SUNT::conjugate() const {
385 if(!is_a<lst>(aij[0])) return SUNT(aij[0], aij[2], aij[1]);
386 lst argv = ex_to<lst>(aij[0]);
387 lst as;
388 for(auto it=argv.rbegin(); it!=argv.rend(); ++it) as.append(*it);
389 return SUNT(as, aij[2], aij[1]);
390 }
391
392 //GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(SUNF, basic,print_func<print_dflt>(&SUNF::print).print_func<FormFormat>(&SUNF::form_print).print_func<FCFormat>(&SUNF::fc_print))
393 GiNaC::registered_class_info & SUNF::get_class_info_static() { return reg_info; }
395 SUNF * SUNF::duplicate() const { SUNF * bp = new SUNF(*this); bp->setflag(GiNaC::status_flags::dynallocated); return bp; }
396 void SUNF::accept(GiNaC::visitor & v) const { if (visitor *p = dynamic_cast<visitor *>(&v)) p->visit(*this); else inherited::accept(v); }
397 const GiNaC::registered_class_info &SUNF::get_class_info() const { return get_class_info_static(); }
398 GiNaC::registered_class_info &SUNF::get_class_info() { return get_class_info_static(); }
399 const char *SUNF::class_name() const { return get_class_info_static().options.get_name(); }
400 //GINAC_IMPLEMENT_REGISTERED_CLASS END
401
402 SUNF::SUNF(ex i, ex j, ex k) : ijk{i,j,k} { }
403 int SUNF::compare_same_type(const basic &other) const {
404 if(!is_a<SUNF>(other)) throw Error("SUNF::compare_same_type");
405 const SUNF &o = static_cast<const SUNF &>(other);
406 for(int i=0; i<3; i++) {
407 auto c = ijk[i].compare(o.ijk[i]);
408 if(c!=0) return c;
409 }
410 return 0;
411 }
412
413 bool SUNF::is_equal_same_type(const basic & other) const {
414 if(!is_a<SUNF>(other)) throw Error("SUNF::is_equal_same_type");
415 const SUNF &o = static_cast<const SUNF &>(other);
416 for(int i=0; i<3; i++) {
417 if(!ijk[i].is_equal(o.ijk[i])) return false;
418 }
419 return true;
420 }
421
422 ex SUNF::eval() const {
423 if(flags & status_flags::evaluated) return *this;
424 if(ijk[0].is_equal(ijk[1]) || ijk[1].is_equal(ijk[2]) || ijk[0].is_equal(ijk[2])) return 0;
425 bool c01 = ex_less(ijk[0],ijk[1]);
426 bool c12 = ex_less(ijk[1],ijk[2]);
427 if(c01 && c12) return this->hold();
428 bool c02 = ex_less(ijk[0],ijk[2]);
429 if(!c01 && c02) return -SUNF(ijk[1],ijk[0],ijk[2]);
430 else if(c02 && !c12) return -SUNF(ijk[0],ijk[2],ijk[1]);
431 else if(!c02 && c01) return SUNF(ijk[2],ijk[0],ijk[1]);
432 else if(c12 && !c02) return SUNF(ijk[1],ijk[2],ijk[0]);
433 else if(!c12 && !c01) return -SUNF(ijk[2],ijk[1],ijk[0]);
434 else return this->hold();
435 }
436
437 void SUNF::print(const print_dflt &c, unsigned) const {
438 c.s << "f(" << ijk[0] << "," << ijk[1] << "," << ijk[2] << ")";
439 }
440
441 void SUNF::form_print(const FormFormat &c, unsigned) const {
442 c << "f(" << ijk[0] << "," << ijk[1] << "," << ijk[2] << ")";
443 }
444
445 void SUNF::fc_print(const FCFormat &c, unsigned) const {
446 c << "SUNF[" << ijk[0] << "," << ijk[1] << "," << ijk[2] << "]";
447 }
448
449 size_t SUNF::nops() const { return 3; }
450 ex SUNF::op(size_t i) const {
451 return ijk[i];
452 }
453 ex & SUNF::let_op(size_t i) {
454 ensure_if_modifiable();
455 return ijk[i];
456 }
457
458 void SUNF::archive(archive_node & n) const {
459 inherited::archive(n);
460 n.add_ex("i", ijk[0]);
461 n.add_ex("j", ijk[1]);
462 n.add_ex("k", ijk[2]);
463 }
464
465 void SUNF::read_archive(const archive_node& n) {
466 inherited::read_archive(n);
467 ex o;
468 n.find_ex("i", o);
469 ijk[0] = ex_to<Index>(o);
470 n.find_ex("j", o);
471 ijk[1] = ex_to<Index>(o);
472 n.find_ex("k", o);
473 ijk[2] = ex_to<Index>(o);
474 }
475
481 ex SUNF::derivative(const symbol & s) const {
482 return 0;
483 }
484
485 //GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(SUNF4, basic,print_func<print_dflt>(&SUNF4::print).print_func<FormFormat>(&SUNF4::form_print).print_func<FCFormat>(&SUNF4::fc_print))
486 GiNaC::registered_class_info & SUNF4::get_class_info_static() { return reg_info; }
488 SUNF4 * SUNF4::duplicate() const { SUNF4 * bp = new SUNF4(*this); bp->setflag(GiNaC::status_flags::dynallocated); return bp; }
489 void SUNF4::accept(GiNaC::visitor & v) const { if (visitor *p = dynamic_cast<visitor *>(&v)) p->visit(*this); else inherited::accept(v); }
490 const GiNaC::registered_class_info &SUNF4::get_class_info() const { return get_class_info_static(); }
491 GiNaC::registered_class_info &SUNF4::get_class_info() { return get_class_info_static(); }
492 const char *SUNF4::class_name() const { return get_class_info_static().options.get_name(); }
493 //GINAC_IMPLEMENT_REGISTERED_CLASS END
494
495 SUNF4::SUNF4(ex i, ex j, ex k, ex l) : ijkl{i,j,k,l} { }
496 int SUNF4::compare_same_type(const basic &other) const {
497 if(!is_a<SUNF4>(other)) throw Error("SUNF4::compare_same_type");
498 const SUNF4 &o = static_cast<const SUNF4 &>(other);
499 for(int i=0; i<4; i++) {
500 auto c = ijkl[i].compare(o.ijkl[i]);
501 if(c!=0) return c;
502 }
503 return 0;
504 }
505
506 bool SUNF4::is_equal_same_type(const basic & other) const {
507 if(!is_a<SUNF4>(other)) throw Error("SUNF4::is_equal_same_type");
508 const SUNF4 &o = static_cast<const SUNF4 &>(other);
509 for(int i=0; i<4; i++) {
510 if(!ijkl[i].is_equal(o.ijkl[i])) return false;
511 }
512 return true;
513 }
514
519 ex SUNF4::eval() const {
520 if(flags & status_flags::evaluated) return *this;
521 if(ijkl[0].is_equal(ijkl[1]) || ijkl[2].is_equal(ijkl[3])) return 0;
522 bool c01 = ex_less(ijkl[0],ijkl[1]);
523 bool c23 = ex_less(ijkl[2],ijkl[3]);
524 if(c01 && c23) return this->hold(); // 01-23
525 else if(!c01 && c23) return -SUNF4(ijkl[1],ijkl[0],ijkl[2],ijkl[3]);
526 else if(!c01 && !c23) return SUNF4(ijkl[1],ijkl[0],ijkl[3],ijkl[2]);
527 else if(c01 && !c23) return -SUNF4(ijkl[0],ijkl[1],ijkl[3],ijkl[2]);
528 else return this->hold();
529 }
530
536 void SUNF4::print(const print_dflt &c, unsigned o) const {
537 c.s << "f(" << ijkl[0] << "," << ijkl[1] << "," << ijkl[2] << "," << ijkl[3] << ")";
538 }
539
540 void SUNF4::fc_print(const FCFormat &c, unsigned o) const {
541 c << "SUNF[" << ijkl[0] << "," << ijkl[1] << "," << ijkl[2] << "," << ijkl[3] << "]";
542 }
543
549 void SUNF4::form_print(const FormFormat &c, unsigned o) const {
550 c << "f4(" << ijkl[0] << "," << ijkl[1] << "," << ijkl[2] << "," << ijkl[3] << ")";
551 }
552
553 size_t SUNF4::nops() const { return 4; }
554 ex SUNF4::op(size_t i) const {
555 return ijkl[i];
556 }
557 ex & SUNF4::let_op(size_t i) {
558 ensure_if_modifiable();
559 return ijkl[i];
560 }
561
566 void SUNF4::archive(archive_node & n) const {
567 inherited::archive(n);
568 n.add_ex("i", ijkl[0]);
569 n.add_ex("j", ijkl[1]);
570 n.add_ex("k", ijkl[2]);
571 n.add_ex("l", ijkl[3]);
572 }
573
578 void SUNF4::read_archive(const archive_node& n) {
579 inherited::read_archive(n);
580 ex o;
581 n.find_ex("i", o);
582 ijkl[0] = ex_to<Index>(o);
583 n.find_ex("j", o);
584 ijkl[1] = ex_to<Index>(o);
585 n.find_ex("k", o);
586 ijkl[2] = ex_to<Index>(o);
587 n.find_ex("l", o);
588 ijkl[3] = ex_to<Index>(o);
589 }
590
596 ex SUNF4::derivative(const symbol & s) const {
597 return 0;
598 }
599
600 ex ncmul_expand(const ex & expr) {
601 MapFunction inner_expand([](const ex & e, MapFunction & self)->ex{
602 if(is_a<add>(e)) {
603 ex res = 0;
604 for(auto ei : e) res += ncmul_expand(ei);
605 return res;
606 } else if(is_a<mul>(e) || is_a<ncmul>(e)) {
607 lst res = lst{ 1 };
608 for(auto ei : e) {
609 ex rei = ncmul_expand(ei);
610 if(!is_a<add>(rei)) rei = lst{ rei };
611 lst ores = res;
612 res = lst{ };
613 for(auto oi : ores) for(auto ri : rei) res.append(oi * ri);
614 }
615 ex ret = 0;
616 for(auto ri : res) ret += ri;
617 return ret;
618 } else return e.map(self);
619 });
620 return inner_expand(expr);
621 }
622
629 ex GMatContract(const ex & expr_in, bool auto_tr) {
630 if(!expr_in.has(GMat(w1,w2,w3))) return expr_in;
631
632 auto expr = expr_in.subs(pow(GMat(w1,w2,w3),2)==GMat(w1,w2,w3)*GMat(w1,w3,w2));
633 expr = expr.subs(GMat(w1,w2,w2)==TR(w1));
634 auto cv_lst = collect_lst(expr, GMat(w1, w2, w3));
635 expr = 0;
636
637 for(auto cv : cv_lst) {
638 auto e = cv.op(1);
639 if(is_zero(e-1) || e.match(GMat(w1, w2, w3))) {
640 if(e.match(GMat(w1, w2, w2))) expr += cv.op(0) * TR(e.op(0));
641 else expr += cv.op(0) * e;
642 continue;
643 } else if(!is_a<mul>(e)) throw Error("GMatContract: collect error: " + ex2str(e));
644
645 lst mats;
646 for(auto item : e) mats.append(item);
647
648 std::map<ex,int,ex_is_less> to_map, from_map;
649 std::set<int> todo;
650 lst mats_idx;
651
652 start:
653 for(int i=0; i<mats.nops(); i++) {
654 auto item = mats.op(i);
655 if(item.op(0).return_type()==return_types::commutative || item.op(0).is_equal(GAS(1))) {
656 mats_idx.append(lst{item,i});
657 } else {
658 if(!item.match(GMat(w1,w2,w3))) {
659 cout << "item in GMatContract: " << item << endl;
660 throw Error("GMatContract faild!");
661 }
662 if(to_map[item.op(1)]!=0 || from_map[item.op(2)]!=0) {
663 if(!auto_tr) throw Error("GMatContract: index conflict for mats: "+ex2str(mats));
664 lst mats2; // to avoid dead-loop
665 mats2.append(GMatT(item));
666 for(int j=0; j<mats.nops(); j++) if(j!=i) mats2.append(mats.op(j));
667 mats = mats2;
668 to_map.clear();
669 from_map.clear();
670 mats_idx.remove_all();
671 goto start;
672 }
673 to_map[item.op(1)] = i+10; // avoid 0 in map
674 from_map[item.op(2)] = i+10; // avoid 0 in map
675 }
676 todo.insert(i);
677 }
678
679 // update to_map/from_map w.r.t mats_idx
680 bool checked = false;
681 while(true) {
682 lst mats_idx2;
683 bool ok = true; // double check
684 for(int i=0; i<mats_idx.nops(); i++) {
685 auto item = mats_idx.op(i).op(0);
686 int ii = ex_to<numeric>(mats_idx.op(i).op(1)).to_int();
687 if(!checked &&
688 to_map[item.op(1)]==0 && from_map[item.op(2)]==0 &&
689 to_map[item.op(2)]==0 && from_map[item.op(1)]==0) {
690 mats_idx2.append(mats_idx.op(i));
691 continue;
692 }
693 ok = false;
694 checked = false;
695 if(to_map[item.op(1)]==0 && from_map[item.op(2)]==0) {
696 to_map[item.op(1)] = ii+10; // avoid 0 in map
697 from_map[item.op(2)] = ii+10; // avoid 0 in map
698 } else if(to_map[item.op(2)]==0 && from_map[item.op(1)]==0) {
699 to_map[item.op(2)] = ii+10; // avoid 0 in map
700 from_map[item.op(1)] = ii+10; // avoid 0 in map
701 // need to swap the 2nd and 3rd index
702 auto li = get_op(mats, ii, 1);
703 auto ri = get_op(mats, ii, 2);
704 let_op(mats, ii, 1, ri);
705 let_op(mats, ii, 2, li);
706 } else {
707 throw Error("GMatContract: index conflict (2).");
708 }
709 }
710 if(mats_idx2.nops()<1) break;
711 mats_idx = mats_idx2;
712 if(ok) checked=true;
713 }
714
715 ex retMat = 1;
716 while(todo.size()>0) {
717 int c = *(todo.begin());
718 todo.erase(c);
719 ex curMat = mats.op(c).op(0);
720 auto li=mats.op(c).op(1);
721 auto ri=mats.op(c).op(2);
722 while(true) {
723 if(li.is_equal(ri)) {
724 retMat *= TR(curMat);
725 break;
726 }
727 int ti = to_map[ri];
728 int fi = from_map[li];
729 if(ti==0 && fi==0) {
730 retMat *= GMat(curMat, li, ri);
731 break;
732 }
733 if(ti!=0) {
734 auto mat = mats.op(ti-10).op(0);
735 if(curMat.is_equal(GAS(1))) curMat = mat;
736 else if(!mat.is_equal(GAS(1))) curMat = curMat * mat;
737 ri = mats.op(ti-10).op(2);
738 todo.erase(ti-10);
739 continue;
740 }
741 if(fi!=0) {
742 auto mat = mats.op(fi-10).op(0);
743 if(curMat.is_equal(GAS(1))) curMat = mat;
744 else if(!mat.is_equal(GAS(1))) curMat = mat * curMat;
745 li = mats.op(fi-10).op(1);
746 todo.erase(fi-10);
747 continue;
748 }
749 }
750 }
751 expr += cv.op(0) * retMat;
752 }
753
754 return expr;
755 }
756
757 ex Contract(const ex & ei) {
758 static exmap cache;
759 if(GMat_using_cache) {
760 auto found = cache.find(ei);
761 if(found!=cache.end()) return found->second;
762 }
763 lst idx_lst;
764 MapFunction get_idx([&idx_lst](const ex & e, MapFunction & self)->ex{
765 if(!Index::has(e) || !Pair::has(e) || e.match(GMat(w1,w2,w3))) return 1; // skip GMat object
766 else if(is_a<Pair>(e)) {
767 if(is_a<Index>(e.op(0)) || is_a<Index>(e.op(1))) idx_lst.append(e);
768 return 1;
769 } else return e.map(self);
770 });
771 get_idx(ei);
772 idx_lst.sort();
773 idx_lst.unique();
774 if(idx_lst.nops()==0) return ei;
775 auto cvs = collect_lst(ei, idx_lst);
776 ex res = 0;
777 for(auto cv : cvs) {
778 auto c = cv.op(0);
779 auto v = form(cv.op(1)); // contract on itself
780 if(!Index::has(v)) res += c * v;
781 else {
782 if(!is_a<mul>(v)) v = lst{ v };
783 exmap repl;
784 ex r = 1; // uncontracted remained index
785 for(auto vi : v) {
786 if(!is_a<Pair>(vi)) r *= vi; // contract may result in a non-Pair object
787 else if(is_a<Index>(vi.op(1)) && c.has(vi.op(1))) repl[vi.op(1)] = vi.op(0);
788 else if(is_a<Index>(vi.op(0)) && c.has(vi.op(0))) repl[vi.op(0)] = vi.op(1);
789 else r *= vi;
790 }
791 res += r * c.subs(repl);
792 }
793 }
794 res = res.subs(SP_map);
795 if(GMat_using_cache) cache[ei] = res;
796 return res;
797 }
798
799 ex GMatOut(const ex & expr_in) {
800 MapFunction inner_out([&](const ex & e, MapFunction & self)->ex {
801 if(e.match(GMat(w1,w2,w3))) {
802 auto e0 = e.op(0);
803 if(is_a<mul>(e0)) {
804 ex c = 1, v = 1;
805 for(auto item : e0) {
806 if(item.return_type()==return_types::commutative) c *= item;
807 else {
808 if(!v.is_equal(1)) {
809 cout << "c=" << c << ", " << "v=" << v << endl;
810 throw Error("GMatOut: v != 1"); // make sure only one non-commutative object
811 }
812 v = item;
813 }
814 }
815 if(v.is_equal(1)) v = GAS(1);
816 return c * GMatOut(GMat(v, e.op(1), e.op(2)));
817 } else return e;
818 } else return e.map(self);
819 });
820 return inner_out(expr_in);
821 }
822
823 ex GMatExpand(const ex & expr_in) {
824 static exmap cache;
825 ex key = expr_in;
826 if(GMat_using_cache) {
827 auto found = cache.find(key);
828 if(found!=cache.end()) return found->second;
829 }
830 MapFunction inner_expand([&](const ex & e, MapFunction & self)->ex {
831 if(!e.has(GMat(w1,w2,w3))) return e;
832 else if(e.match(GMat(w1,w2,w3))) {
833 auto e0 = e.op(0);
834 if(is_a<add>(e0)) {
835 ex res = 0;
836 for(auto item : e0) res += GMatExpand(GMat(item, e.op(1), e.op(2)));
837 return res;
838 } else if(is_a<mul>(e0)) {
839 ex c = 1, v = 1;
840 for(auto item : e0) {
841 if(item.return_type()==return_types::commutative) c *= item;
842 else {
843 if(!v.is_equal(1)) {
844 cout << "c=" << c << ", " << "v=" << v << endl;
845 throw Error("GMatExpand: v != 1"); // make sure only one non-commutative object
846 }
847 v = item;
848 }
849 }
850 if(v.is_equal(1)) v = GAS(1);
851 return c * GMatExpand(GMat(v, e.op(1), e.op(2)));
852 } else if(is_a<ncmul>(e0)) { // expand ncmul
853 ex res;
854 bool first = true;
855 for(auto item : e0) {
856 if(first) {
857 res = item;
858 first = false;
859 continue;
860 }
861 ex ncL = res; // previous result
862 if(!is_a<add>(ncL)) ncL = lst{ ncL };
863 ex ncR = item;
864 if(!is_a<add>(ncR)) ncR = lst{ ncR };
865 res = 0; // current result
866 for(auto iL : ncL) for(auto iR : ncR) res += iL * iR;
867 }
868 ex rs = res;
869 res = 0;
870 if(!is_a<add>(rs)) rs = lst{ rs };
871 for(auto item : rs) { // pull out commutative coefficient
872 ex c = 1, v = 1;
873 if(is_a<mul>(item)) {
874 if(item.nops()==1) throw Error("GMatExpand: item.nops == 1"); // make sure
875 for(auto it : item) {
876 if(it.return_type()==return_types::commutative) c *= it;
877 else {
878 if(!v.is_equal(1)) throw Error("GMatExpand: v != 1"); // make sure only one non-commutative object
879 v = it;
880 }
881 }
882 } else v = item;
883
884 while(true) { // recursive replace ɣ.P * ɣ.P -> P^2 and ɣ.mu * ɣ.mu -> d @ v
885 bool to_exit = true;
886 if(is_a<ncmul>(v)) {
887 bool first = true;
888 ex last = 1, vv = 1;
889 for(auto vi : v) {
890 if(first) {
891 first = false;
892 last = vi;
893 } else {
894 if(last==vi && is_a<DGamma>(vi)) {
895 first = true;
896 last = 1;
897 if(is_a<Vector>(vi.op(0))) c *= SP(vi.op(0));
898 else if(is_a<Index>(vi.op(0))) c *= d;
899 else if(vi.op(0).is_equal(1) || vi.op(0).is_equal(5)) c *= 1; // GAS(1)*GAS(1) = GAS(5)*GAS(5) = 1
900 else throw Error("GMatExpand: only GAS(i/p/1/5) supported.");
901 to_exit = false; // need to cycle again
902 } else {
903 if(last!=GAS(1) && !last.is_equal(1)) vv = vv * last;
904 last = vi;
905 }
906 }
907 }
908 if(!last.is_equal(1) && last!=GAS(1)) v = vv * last; // check last item
909 else v = vv;
910 }
911 if(to_exit) break;
912 }
913 if(v.is_equal(1)) v = GAS(1); // identity matrix
914 res += c * GMat(v, e.op(1), e.op(2));
915 }
916 return res;
917 } else return e;
918 } else return e.map(self);
919 });
920 ex res = inner_expand(expr_in);
921 if(GMat_using_cache) cache[key] = res;
922 return res;
923 }
924
925 ex GMatShift(const ex & expr, const ex & g, bool to_right) {
926 if(!expr.has(g)) return expr;
927 static exmap cache;
928 ex key = lst{expr, g, to_right ? 1 : 0};
929 if(GMat_using_cache) {
930 auto found = cache.find(key);
931 if(found!=cache.end()) return found->second;
932 }
933 MapFunction inner_shift([g,to_right](const ex & e, MapFunction & self)->ex{
934 if(!e.has(g) || !e.has(GMat(w1,w2,w3))) return e;
935 else if(e.match(GMat(w1,w2,w3))) {
936 ex eg = e.op(0);
937 if(!is_a<ncmul>(eg)) eg = lst{ eg };
938 int gi = -1;
939 if(to_right) {
940 for(int i=0; i<eg.nops()-1; i++) if(eg.op(i)==g) { gi = i; break; }
941 if(gi==-1) return e;
942 } else {
943 for(int i=eg.nops()-1; i>0; i--) if(eg.op(i)==g) { gi = i; break; }
944 if(gi==-1) return e;
945 }
946 int gj = gi + ( to_right ? 1 : -1 );
947 ex rem = 1, rem2 = 1;
948 for(int i=0; i<eg.nops(); i++) {
949 if(i!=gi && i!=gj) {
950 rem *= eg.op(i);
951 rem2 *= eg.op(i);
952 }
953 if(i==gi) {
954 if(to_right) rem2 *= eg.op(gj)*eg.op(gi);
955 else rem2 *= eg.op(gi)*eg.op(gj);
956 }
957 }
958 if(eg.op(gi).is_equal(eg.op(gj))) {
959 ex ip = eg.op(gi).op(0);
960 if(rem.is_equal(1)) rem = GAS(1);
961 ex res = GMat(rem, e.op(1), e.op(2));
962 res = GMatShift(res, g, to_right);
963 ex c;
964 if(is_a<Vector>(ip)) c = SP(ip);
965 else if(is_a<Index>(ip)) c = d;
966 else if(eg.op(gi).is_equal(GAS(5))) c = 1;
967 else throw Error("GMatShift: only GAS(i/p/5) supproted.");
968 return c * res;
969 }
970 if(rem.is_equal(1)) rem = GAS(1);
971 if(rem2.is_equal(1)) rem2 = GAS(1);
972 ex res = 0;
973 if(!eg.op(gi).is_equal(GAS(5)) && !eg.op(gj).is_equal(GAS(5))) {
974 res = 2*SP(eg.op(gi).op(0), eg.op(gj).op(0)) * GMat(rem, e.op(1), e.op(2));
975 }
976 res = res - GMat(rem2, e.op(1), e.op(2));
977 return GMatShift(res, g, to_right);
978 } else return e.map(self);
979 });
980 ex res = GMatExpand(Contract(expr)); // add Contract & GMatExpand here
981 res = collect_ex(res, GMat(w1,w2,w3));
982 res = inner_shift(res);
983 if(GMat_using_cache) cache[key] = res;
984 return res;
985 }
986
987 ex GMatSimplify(const ex & expr) {
988 ex res = GMatContract(expr);
989 res = GMatShift(res);
990 res = Contract(res);
991 return res;
992 }
993
994 namespace {
995 // 1st and last should be equal
996 lst shift_1st_to_right(const ex & e) { // return a list of {coeff, gammas}
997 if(!is_a<ncmul>(e)) throw Error("input is not a ncmul.");
998 if(e.nops()==2) {
999 ex e0 = e.op(0);
1000 if(e.op(0)!=e.op(1)) throw Error("shift_1st_to_right: the 2 items are not equal!");
1001 else if(is_a<Index>(e0.op(0))) return lst{ lst{ d, 1 }};
1002 else if(is_a<Vector>(e0.op(0))) return lst{ lst{ SP(e.op(0).op(0)), 1 }};
1003 else if(e0.is_equal(GAS(5))) return lst{ lst{ 1, 1 }};
1004 else {
1005 cout << endl << e << endl;
1006 throw Error("shift_1st_to_right: only GAS(i/p/5) supproted.");
1007 }
1008 }
1009 ex rem = 1;
1010 int n = e.nops();
1011 for(int i=2; i<n; i++) rem *= e.op(i);
1012 lst res = shift_1st_to_right(e.op(0)*rem);
1013 n = res.nops();
1014 for(int i=0; i<n; i++) {
1015 res[i][0] = -res.op(i).op(0);
1016 res[i][1] = e.op(1) * res.op(i).op(1);
1017 }
1018 if(!e.op(0).is_equal(GAS(5)) && !e.op(1).is_equal(GAS(5))) {
1019 if(!is_a<Index>(e.op(0).op(0)) && !is_a<Vector>(e.op(0).op(0))) {
1020 cout << e << endl;
1021 throw Error("shift_12_right: not a Vector or Index");
1022 }
1023 if(!is_a<Index>(e.op(1).op(0)) && !is_a<Vector>(e.op(1).op(0))) {
1024 cout << e << endl;
1025 throw Error("shift_12_right: not a Vector or Index");
1026 }
1027 res.append(lst{ 2*SP(e.op(0).op(0), e.op(1).op(0)), rem });
1028 }
1029 return res;
1030 }
1031 }
1032 ex GMatShift(const ex & expr) {
1033 static exmap cache;
1034 if(GMat_using_cache) {
1035 auto found = cache.find(expr);
1036 if(found!=cache.end()) return found->second;
1037 }
1038
1039 MapFunction inner_shift([](const ex & e, MapFunction & self)->ex{
1040 if(e.match(GMat(w1,w2,w3))) {
1041 ex eg = e.op(0);
1042 if(!is_a<ncmul>(eg)) eg = lst{ eg };
1043
1044 int gi = -1, gj = -1;
1045 for(int i=0; i<eg.nops(); i++) for(int j=i+1; j<eg.nops(); j++) {
1046 if(eg.op(i).is_equal(eg.op(j))) {
1047 gi = i;
1048 gj = j;
1049 goto done;
1050 }
1051 }
1052 return e;
1053 done: ;
1054
1055 ex exL = 1, exM=1, exR = 1;
1056 for(int i=0; i<eg.nops(); i++) {
1057 if(i<gi) exL *= eg.op(i);
1058 else if(i>gj) exR *= eg.op(i);
1059 else exM *= eg.op(i);
1060 }
1061 lst cvs = shift_1st_to_right(exM);
1062
1063 ex res = 0;
1064 for(auto cv : cvs) {
1065 ex item = exL*cv.op(1)*exR;
1066 if(item.is_equal(1)) item = GAS(1);
1067 res += cv.op(0)*GMatShift(GMat(item, e.op(1), e.op(2)));
1068 }
1069
1070 return res;
1071 } else return e.map(self);
1072 });
1073 ex res = GMatExpand(Contract(expr)); // add Contract & GMatExpand here
1074 res = collect_ex(res, GMat(w1,w2,w3));
1075 res = inner_shift(res);
1076 if(GMat_using_cache) cache[expr] = res;
1077 return res;
1078 }
1079
1080 ex GMatECC(const ex & expr, int sign) {
1081 if(!expr.has(DGamma::C)) return expr;
1082 MapFunction inner_ecc([sign](const ex & e, MapFunction & self)->ex{
1083 if(!e.has(DGamma::C)) return e;
1084 else if(e.match(GMat(w1,w2,w3)) || e.match(TR(w))) {
1085 ex eg = e.op(0), cc = 1;
1086 if(is_a<mul>(e.op(0))) {
1087 eg = 1;
1088 for(auto item : e.op(0)) {
1089 if(item.return_type()==return_types::noncommutative) {
1090 if(eg.is_equal(1)) eg = item;
1091 else throw Error("GMatECC:: 2 more noncommutative objects found.");
1092 } else cc *= item;
1093 }
1094 if(eg.is_equal(1)) throw Error("GMatECC:: eg is 1, NOT expected.");
1095 }
1096 if(!is_a<ncmul>(eg)) eg = lst{ eg }; // only one item
1097 int ci = -1;
1098 for(int i=0; i<eg.nops(); i++) if(eg.op(i)==DGamma::C) { ci = i; break; }
1099 if(ci==-1) return e; // not found C
1100 int cj = -1;
1101 for(int i=ci+1; i<eg.nops(); i++) if(eg.op(i)==DGamma::C) { cj = i; break; }
1102 if(cj==-1) return e; // not found C
1103 int cnt = 0; // remaining C
1104 for(int i=cj+1; i<eg.nops(); i++) if(eg.op(i)==DGamma::C) { cnt++; }
1105 ex res = 1;
1106 for(int i=0; i<ci; i++) res *= eg.op(i);
1107 ex m = 1;
1108 for(int i=ci+1; i<cj; i++) m *= eg.op(i);
1109 if(sign<0) cc *= -1; // C = sign * C^{-1}
1110 if(!m.is_equal(1)) res *= gamma_transpose(charge_conjugate(m));
1111 for(int i=cj+1; i<eg.nops(); i++) res *= eg.op(i);
1112 if(e.nops()==3) res = cc * GMat(res, e.op(1), e.op(2));
1113 else res = cc * TR(res);
1114 return cnt<2 ? res : self(res);
1115 } else return e.map(self);
1116 });
1117 ex res = collect_ex(expr, GMat(w1,w2,w3));
1118 res = inner_ecc(res);
1119 if(res.has(TR(w))) { // replace TR(e^T) = TR(e)
1120 res = MapFunction([](const ex & e, MapFunction & self)->ex{
1121 if(!e.has(TR(w))) return e;
1122 else if(e.match(TR(w))) {
1123 auto gs = DGamma::all(e.op(0));
1124 bool ok = true;
1125 for(auto item : gs) {
1126 auto gi = ex_to<DGamma>(item);
1127 if(!item.is_equal(DGamma::C) && !gi.isTr) {
1128 ok = false;
1129 break;
1130 }
1131 }
1132 if(ok) return TR(gamma_transpose(e.op(0)));
1133 else return e;
1134 } else return e.map(self);
1135 })(res);
1136 }
1137 return res;
1138 }
1139
1140 ex GMatT(const ex & expr) {
1141 MapFunction inner_transpose([](const ex & e, MapFunction & self)->ex{
1142 if(!e.has(GMat(w1,w2,w3))) return e;
1143 else if(e.match(GMat(w1,w2,w3))) {
1144 return GMat(gamma_transpose(e.op(0)), e.op(2), e.op(1));
1145 } else return e.map(self);
1146 });
1147 ex res = collect_ex(expr, GMat(w1,w2,w3));
1148 res = inner_transpose(res);
1149 return res;
1150 }
1151
1152 namespace {
1153 void GMat_fc_print(const ex &arg1, const ex &arg2, const ex &arg3, const print_context &c0) {
1154 auto c = static_cast<const FCFormat &>(c0);
1155 c << "GMat[" << arg1 << "," << arg2 << "," << arg3 << "]";
1156 }
1157 }
1158
1159 REGISTER_FUNCTION(GMat, do_not_evalf_params().print_func<FCFormat>(&GMat_fc_print).conjugate_func(mat_conj).set_return_type(return_types::commutative))
1160
1161 bool IsZero(const ex & e) {
1162 try {
1163 exset vs;
1164 for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i) {
1165 if(is_a<symbol>(*i) || is_a<Pair>(*i)) vs.insert(*i);
1166 }
1167
1168 int n = 13;
1169 for(int i=0; i<5; i++) {
1170 exmap nsubs;
1171 for(auto item : vs) {
1172 nsubs[item] = ex(1)/n_nth_prime(n);
1173 n++;
1174 }
1175 ex ret = e.subs(nsubs);
1176 if(!normal(e).is_zero()) return false;
1177 }
1178
1179 auto ret = exnormal(e);
1180 return ret.is_zero();
1181 } catch(...) { }
1182 return is_zero(exnormal(e));
1183 }
1184
1185 ex ToCF(const ex & e) {
1186 ex res = e;
1187 bool todo = true;
1188 while(todo) {
1189 todo = false;
1190 auto cvs = collect_lst(res, lst{NF,TF});
1191 res = 0;
1192 for(auto cv : cvs) {
1193 int degTF = cv.op(1).degree(TF);
1194 int degNF = cv.op(1).ldegree(NF);
1195 if(degTF>0 && degNF<0) {
1196 todo = true;
1197 int n = degTF;
1198 if(degTF+degNF>0) n = -degNF;
1199 res += cv.op(0) * cv.op(1) * pow(TF/NF,-n) * pow(TF*NF-CF,n);
1200 } else if(degTF>0 && degNF>1) {
1201 todo = true;
1202 int n = degTF;
1203 if(degTF>degNF/2) n = degNF/2;
1204 res += cv.op(0) * cv.op(1) * pow(TF*NF*NF,-n) * pow(CF*NF+TF,n);
1205 } else res += cv.op(0) * cv.op(1);
1206 }
1207 }
1208 return res;
1209 }
1210
1211 ex ca_neg_pow_sub(const ex & expr) {
1212 static MapFunction ca_map([](const ex & e, MapFunction & self)->ex {
1213 if(!e.has(CA)) return e;
1214 else if(e.match(pow(CA,w)) && e.op(1).info(info_flags::negint)) return pow(CA-2*CF,-e.op(1));
1215 else return e.map(self);
1216 });
1217 return ca_map(expr);
1218 }
1219
1220 ex ToCACF(const ex & e) { // from FeynCalc
1221 ex res = e.subs(lst{NA==(NF*NF-1),CF==(NF*NF-1)/(2*NF),TF==ex(1)/2});
1222 res = exfactor(res);
1223 // SUNN -> CA
1224 res = res.subs(NF==CA);
1225 // (-1+CA^2)->(-2 CA CF)
1226 res = res.subs(lst{ w*(1-CA)*(1+CA)==-w*2*CA*CF, w*(-1+CA)*(1+CA)==w*2*CA*CF });
1227 res = res.subs(lst{ w1*pow(1-CA,w2)*pow(1+CA,w2)==w1*pow(-2*CA*CF,w2), w1*pow(-1+CA,w2)*pow(1+CA,w2)==w1*pow(2*CA*CF,w2) });
1228 // (((2 - CA^2) CF )/CA ) ->(CF (CA - 4 CF))
1229 res = res.subs(lst{ w*(2-CA*CA)*CF/CA==w*CF*(CA-4*CF), w*(-2+CA*CA)*CF/CA==w*CF*(-CA+4*CF) });
1230 // (1-CA^2) -> (-2 CA CF)
1231 res = res.subs(lst{ w*(1-CA)*(1+CA)==-w*2*CA*CF, w*(-1+CA)*(1+CA)==w*2*CA*CF });
1232 res = res.subs(lst{ w1*pow(1-CA,w2)*pow(1+CA,w2)==-w1*pow(2*CA*CF,w2), w1*pow(-1+CA,w2)*pow(1+CA,w2)==w1*pow(2*CA*CF,w2) });
1233 // (1/CA)^n -> (CA - 2 CF)^n
1234 //res = res.subs(lst{ w/CA==w*(CA-2*CF) });
1235 res = ca_neg_pow_sub(res);
1236 // ((1 - CA^2)*(CA - 2*CF)) -> (-2*CF)
1237 res = res.subs(lst{ w*(1-CA)*(1+CA)*(CA-2*CF)==-2*w*CF, w*(-1+CA)*(1+CA)*(CA-2*CF)==2*w*CF });
1238 // (CA (CA-2 CF)) -> 1
1239 res = res.subs(lst{ w*CA*(CA-2*CF)==w, w*CA*(-CA+2*CF)==-w });
1240 res = res.subs(lst{ w1*pow(CA,w2)*pow(CA-2*CF,w2)==w1, w1*pow(CA,w2)*pow(-CA+2*CF,w2)==w*pow(-1,w2) });
1241 // (CA^2+c)(CA-2CF) -> CA+c(CA-2CF)
1242 res = res.subs(lst{ w0*(CA*CA+w1)*(CA-2*CF)==w0*(CA+w1*(CA-2*CF)), w0*(CA*CA+w1)*(-CA+2*CF)==-w0*(CA+w1*(CA-2*CF)) });
1243 res = res.subs(lst{ w0*(CA*CA+w1)*pow(CA-2*CF,w2)==w0*(CA+w1*(CA-2*CF))*pow(CA-2*CF,w2-1), w0*(CA*CA+w1)*pow(-CA+2*CF,w2)==-w0*(CA+w1*(CA-2*CF))*pow(-CA+2*CF,w2-1) });
1244 return res;
1245 }
1246
1247 ex HomCACF(const ex & e) {
1248 ex res = e.subs(lst{NA==(NF*NF-1),CA==NF,CF==(NF*NF-1)/(2*NF),TF==ex(1)/2});
1249 res = exfactor(res);
1250 if(!is_a<mul>(res)) res = lst{ res };
1251 ex c=1, v=1;
1252 for(auto item : res) {
1253 if(item.has(NF)) c *= item;
1254 else v *= item;
1255 }
1256 c = collect_ex(c, NF);
1257 int deg = c.degree(NF);
1258 int ldeg = -c.ldegree(NF);
1259 if(ldeg>deg) deg = ldeg;
1260 if(deg>0) {
1261 lst vars;
1262 ex eqn = c;
1263 for(int i=0; i<=deg; i++) {
1264 symbol xi;
1265 eqn -= xi * pow(NF,i) * pow((NF*NF-1)/(2*NF), deg-i);
1266 vars.append(xi);
1267 }
1268 eqn = collect_ex(eqn, NF);
1269 int nH = eqn.degree(NF);
1270 int nL = eqn.ldegree(NF);
1271 lst eqns;
1272 for(int i=nL; i<=nH; i++) {
1273 ex cc = eqn.coeff(NF, i);
1274 if(cc.is_zero()) continue;
1275 eqns.append(cc==0);
1276 }
1277 auto sol = lsolve(eqns, vars);
1278 if(sol.nops()!=deg+1) {
1279 cout << "c=" << c << endl;
1280 cout << "sol=" << sol << endl;
1281 throw Error("HomCACF: no solution found!");
1282 }
1283 c = 0;
1284 for(int i=0; i<=deg; i++) c += vars.op(i).subs(sol) * pow(CA,i) * pow(CF, deg-i);
1285 }
1286 return c * v;
1287 }
1288
1289 ex DoColor(const ex & expr, const ex & pref, int method) {
1290 auto cvs = collect_lst(expr, [](const ex &e)->bool{ return Index::hasc(e); });
1291 ex res = 0;
1292 for(auto const & cv : cvs) {
1293 auto cc = cv.op(0);
1294 auto vv = cv.op(1);
1295 if(method==0) vv = HomCACF(form(vv)/pref)*pref;
1296 else vv = ToCACF(form(vv)/pref)*pref;
1297 res += cc * vv;
1298 }
1299 return res;
1300 }
1301}
1302
int * a
#define IMPLEMENT_HAS(classname)
Definition BASIC.h:24
#define DEFAULT_CTOR(classname)
Definition BASIC.h:21
#define IMPLEMENT_ALL(classname)
Definition BASIC.h:30
HEP header file.
static lst all(const ex &e)
static ex C
Definition HEP.h:487
static bool has(const ex &e)
class used to wrap error message
Definition BASIC.h:242
class for FCFormat Output
Definition HEP.h:71
static void ncmul_print(const ncmul &p, const FCFormat &c, unsigned level=0)
Definition Basic.cpp:163
FCFormat(ostream &os, unsigned opt=0)
Definition Basic.cpp:66
const FCFormat & operator<<(const T &v) const
Definition HEP.h:77
class for FormFormat Output
Definition HEP.h:44
static void power_print(const power &p, const FormFormat &c, unsigned level=0)
Definition Basic.cpp:55
const FormFormat & operator<<(const T &v) const
Definition HEP.h:50
FormFormat(ostream &os, unsigned opt=0)
Definition Basic.cpp:34
class for index object
Definition HEP.h:104
int compare_same_type(const GiNaC::basic &other) const override
Definition Basic.cpp:181
const char * class_name() const override
Definition Basic.cpp:177
static bool hasv(const ex &e)
Definition Basic.cpp:240
void archive(archive_node &n) const override
Definition Basic.cpp:210
bool is_equal_same_type(const basic &other) const override
Definition Basic.cpp:190
static GiNaC::registered_class_info & get_class_info_static()
Definition Basic.cpp:171
Symbol name
Definition HEP.h:132
void print(const print_context &c, unsigned level=0) const
Definition Basic.cpp:198
void read_archive(const archive_node &n) override
Definition Basic.cpp:216
static bool hasc(const ex &e)
Definition Basic.cpp:234
ex derivative(const symbol &s) const override
Definition Basic.cpp:224
const GiNaC::registered_class_info & get_class_info() const override
Definition Basic.cpp:175
Pair operator()(const Index &i)
Definition Basic.cpp:202
void accept(GiNaC::visitor &v) const override
Definition Basic.cpp:174
Index * duplicate() const override
Definition Basic.cpp:173
static bool has(const ex &e)
class to wrap map_function of GiNaC
Definition BASIC.h:632
class for Pair object
Definition HEP.h:322
static bool has(const ex &e)
class for SUNF4 object
Definition HEP.h:278
void archive(archive_node &n) const override
save to archvie
Definition Basic.cpp:566
static GiNaC::registered_class_info & get_class_info_static()
Definition Basic.cpp:486
void form_print(const FormFormat &c, unsigned level=0) const
print the Form Format
Definition Basic.cpp:549
void accept(GiNaC::visitor &v) const override
Definition Basic.cpp:489
const GiNaC::registered_class_info & get_class_info() const override
Definition Basic.cpp:490
size_t nops() const override
Definition Basic.cpp:553
void fc_print(const FCFormat &c, unsigned level=0) const
Definition Basic.cpp:540
void print(const print_dflt &c, unsigned level=0) const
normal priint
Definition Basic.cpp:536
ex op(size_t i) const override
Definition Basic.cpp:554
ex & let_op(size_t i) override
Definition Basic.cpp:557
SUNF4 * duplicate() const override
Definition Basic.cpp:488
ex ijkl[4]
Definition HEP.h:303
ex eval() const override
automatical evaluation of SUNF4
Definition Basic.cpp:519
ex derivative(const symbol &s) const override
set derivative of SUNF4 to 0
Definition Basic.cpp:596
bool is_equal_same_type(const basic &other) const override
Definition Basic.cpp:506
const char * class_name() const override
Definition Basic.cpp:492
void read_archive(const archive_node &n) override
read from archive
Definition Basic.cpp:578
int compare_same_type(const GiNaC::basic &other) const override
Definition Basic.cpp:496
class for SUNF object
Definition HEP.h:233
void archive(archive_node &n) const override
Definition Basic.cpp:458
void print(const print_dflt &c, unsigned level=0) const
Definition Basic.cpp:437
const char * class_name() const override
Definition Basic.cpp:399
int compare_same_type(const GiNaC::basic &other) const override
Definition Basic.cpp:403
ex derivative(const symbol &s) const override
set derivative of SUNF to 0
Definition Basic.cpp:481
bool is_equal_same_type(const basic &other) const override
Definition Basic.cpp:413
void form_print(const FormFormat &c, unsigned level=0) const
Definition Basic.cpp:441
static GiNaC::registered_class_info & get_class_info_static()
Definition Basic.cpp:393
ex & let_op(size_t i) override
Definition Basic.cpp:453
SUNF * duplicate() const override
Definition Basic.cpp:395
ex eval() const override
Definition Basic.cpp:422
const GiNaC::registered_class_info & get_class_info() const override
Definition Basic.cpp:397
void read_archive(const archive_node &n) override
Definition Basic.cpp:465
size_t nops() const override
Definition Basic.cpp:449
void accept(GiNaC::visitor &v) const override
Definition Basic.cpp:396
void fc_print(const FCFormat &c, unsigned level=0) const
Definition Basic.cpp:445
ex op(size_t i) const override
Definition Basic.cpp:450
ex ijk[3]
Definition HEP.h:258
class for SUNT object
Definition HEP.h:189
static GiNaC::registered_class_info & get_class_info_static()
Definition Basic.cpp:305
bool is_equal_same_type(const basic &other) const override
Definition Basic.cpp:325
ex & let_op(size_t i) override
Definition Basic.cpp:357
void archive(archive_node &n) const override
Definition Basic.cpp:362
size_t nops() const override
Definition Basic.cpp:353
void form_print(const FormFormat &c, unsigned level=0) const
Definition Basic.cpp:334
ex op(size_t i) const override
Definition Basic.cpp:354
const GiNaC::registered_class_info & get_class_info() const override
Definition Basic.cpp:309
void fc_print(const FCFormat &c, unsigned level=0) const
Definition Basic.cpp:345
void print(const print_dflt &c, unsigned level=0) const
Definition Basic.cpp:349
int compare_same_type(const GiNaC::basic &other) const override
Definition Basic.cpp:315
void read_archive(const archive_node &n) override
Definition Basic.cpp:369
SUNT * duplicate() const override
Definition Basic.cpp:307
const char * class_name() const override
Definition Basic.cpp:311
void accept(GiNaC::visitor &v) const override
Definition Basic.cpp:308
ex derivative(const symbol &s) const override
Definition Basic.cpp:380
ex aij[3]
Definition HEP.h:214
ex conjugate() const override
Definition Basic.cpp:384
class extended to GiNaC symbol class, represent a positive symbol
Definition BASIC.h:113
class for vector object
Definition HEP.h:149
Vector * duplicate() const override
Definition Basic.cpp:252
bool is_equal_same_type(const basic &other) const override
Definition Basic.cpp:269
Pair operator()(const Vector &p)
Definition Basic.cpp:279
const char * class_name() const override
Definition Basic.cpp:256
ex derivative(const symbol &s) const override
Definition Basic.cpp:300
void print(const print_context &c, unsigned level=0) const
Definition Basic.cpp:275
static GiNaC::registered_class_info & get_class_info_static()
Definition Basic.cpp:250
int compare_same_type(const GiNaC::basic &other) const override
Definition Basic.cpp:260
void read_archive(const archive_node &n) override
Definition Basic.cpp:292
void accept(GiNaC::visitor &v) const override
Definition Basic.cpp:253
void archive(archive_node &n) const override
Definition Basic.cpp:287
const GiNaC::registered_class_info & get_class_info() const override
Definition Basic.cpp:254
Symbol name
Definition HEP.h:176
HepLib namespace.
Definition BASIC.cpp:17
ex exfactor(const ex &expr_in, int opt)
factorize a expression
Definition BASIC.cpp:1854
const Symbol NA
bool ex_less(const ex &a, const ex &b)
Definition Sort.cpp:10
ex DoColor(const ex &expr, const ex &pref, int method)
Definition Basic.cpp:1289
void let_op(ex &ex_in, int index1, int index2, const ex item)
update index1-th.index2-th of expression with item
Definition BASIC.cpp:1560
ex HomCACF(const ex &e)
Definition Basic.cpp:1247
const Symbol gs
ex w0
Definition BASIC.h:499
const Symbol NF
ex GMatShift(const ex &expr, const ex &g, bool to_right)
Definition Basic.cpp:925
exmap SP_map
Definition Init.cpp:184
ex collect_ex(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica
Definition BASIC.cpp:1202
ex GMatSimplify(const ex &expr)
Definition Basic.cpp:987
ex GAS(const ex &expr, unsigned rl)
function similar to GAD/GSD in FeynClac
Definition DGamma.cpp:280
ex GMatExpand(const ex &expr_in)
Definition Basic.cpp:823
ex GMatContract(const ex &expr_in, bool auto_tr)
make contract on matrix, i.e., GMat(a,i1,i2)*GMat(b,i2,i3) -> GMat(a*b,i1,i3)
Definition Basic.cpp:629
ex Contract(const ex &ei)
Definition Basic.cpp:757
ex w
Definition Init.cpp:93
const Symbol vs
const Symbol d
ex get_op(const ex ex_in, int index1, int index2)
return index1-th.index2-th of expression
Definition BASIC.cpp:1606
ex ToCACF(const ex &e)
Definition Basic.cpp:1220
ex GMatECC(const ex &expr, int sign)
Definition Basic.cpp:1080
ex GMatOut(const ex &expr_in)
Definition Basic.cpp:799
const Symbol CA
const Symbol TF
ex exnormal(const ex &expr, int opt)
normalize a expression
Definition BASIC.cpp:1916
ex gamma_transpose(const ex &expr)
make the transpose operaton M --> M^T
Definition Form.cpp:685
const Symbol CF
const Symbol as
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:639
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:1222
ex ToCF(const ex &e)
Definition Basic.cpp:1185
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition BASIC.cpp:715
bool IsZero(const ex &e)
const Symbol nL
const Symbol nH
ex w1
Definition BASIC.h:499
ex w3
Definition BASIC.h:499
ex ncmul_expand(const ex &expr)
Definition Basic.cpp:600
REGISTER_FUNCTION(GMat, do_not_evalf_params().print_func< FCFormat >(&GMat_fc_print).conjugate_func(mat_conj).set_return_type(return_types::commutative)) bool IsZero(const ex &e)
Definition Basic.cpp:1159
ex w2
Definition BASIC.h:499
bool GMat_using_cache
Definition Init.cpp:160
ex GMatT(const ex &expr)
Definition Basic.cpp:1140
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:565
ex SP(const ex &a, bool use_map=false)
Definition Pair.cpp:166
ex ca_neg_pow_sub(const ex &expr)
Definition Basic.cpp:1211