7#include "flint/ulong_extras.h"
36 GINAC_IMPLEMENT_PRINT_CONTEXT(FormFormat, print_dflt)
57 c <<
"((" << p.op(0) <<
")*(" << p.op(0) <<
"))";
59 c <<
"((" << p.op(0) <<
")^(" << p.op(1) <<
"))";
68 GINAC_IMPLEMENT_PRINT_CONTEXT(FCFormat, print_dflt)
91 for(
int r=0; r<nr; r++) {
93 for(
int c=0; c<nc; c++) {
94 mat(r,c).print(*
this);
107 if (i==vend) { s <<
"{}";
return *
this; }
122 if (i==send) { s <<
"{}";
return *
this; }
137 if (i==mend) { s <<
"{}";
return *
this; }
140 i->first.print(*
this);
142 i->second.print(*
this);
152 class ncmul_hack :
public ncmul {
154 ncmul_hack(ncmul _nm) : ncmul(_nm){ }
155 void print(
const FCFormat & c,
unsigned level) {
156 printseq(c,
'(',
'.',
')', precedence(), level);
159 ex mat_conj(
const ex & e1,
const ex & e2,
const ex & e3) {
160 return GMat(e1.conjugate(), e3, e2);
164 ncmul_hack(nm).print(c, level);
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());
186 else if(ret>0)
return 1;
187 else return dim.compare(o.
dim);
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);
203 return Pair(*
this, i);
207 return Pair(p, *
this);
211 inherited::archive(n);
212 n.add_string(
"name",
name.get_name());
213 n.add_ex(
"dim",
dim);
217 inherited::read_archive(n);
219 n.find_string(
"name", nstr);
221 n.find_ex(
"dim",
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;
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;
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;
261 if(!is_a<Vector>(other))
throw Error(
"Vector::compare_same_type");
263 auto ret =
name.get_name().compare(o.
name.get_name());
265 else if(ret<0)
return -1;
270 if(!is_a<Vector>(other))
throw Error(
"Vector::is_equal_same_type");
272 return (
name.get_name() == o.
name.get_name());
280 return Pair(*
this, p);
284 return Pair(*
this, i);
288 inherited::archive(n);
289 n.add_string(
"name",
name.get_name());
293 inherited::read_archive(n);
296 n.find_string(
"name", nstr);
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]);
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;
335 if(is_a<lst>(
aij[0])) {
337 for(
auto item :
aij[0]) {
338 if(first) { first=
false; c <<
"T(" << item; }
339 else c <<
"," << item;
341 }
else c <<
"T(" <<
aij[0];
342 c <<
"," <<
aij[1] <<
"," <<
aij[2] <<
")";
346 c <<
"SUNTF[" <<
aij[0] <<
"," <<
aij[1] <<
"," <<
aij[2] <<
"]";
350 c.s <<
"T(" <<
aij[0] <<
"," <<
aij[1] <<
"," <<
aij[2] <<
")";
358 ensure_if_modifiable();
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]);
370 inherited::read_archive(n);
373 aij[0] = ex_to<Index>(o);
375 aij[1] = ex_to<Index>(o);
377 aij[2] = ex_to<Index>(o);
386 lst argv = ex_to<lst>(
aij[0]);
388 for(
auto it=argv.rbegin(); it!=argv.rend(); ++it)
as.append(*it);
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]);
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;
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;
427 if(c01 && c12)
return this->hold();
434 else return this->hold();
438 c.s <<
"f(" <<
ijk[0] <<
"," <<
ijk[1] <<
"," <<
ijk[2] <<
")";
442 c <<
"f(" <<
ijk[0] <<
"," <<
ijk[1] <<
"," <<
ijk[2] <<
")";
446 c <<
"SUNF[" <<
ijk[0] <<
"," <<
ijk[1] <<
"," <<
ijk[2] <<
"]";
454 ensure_if_modifiable();
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]);
466 inherited::read_archive(n);
469 ijk[0] = ex_to<Index>(o);
471 ijk[1] = ex_to<Index>(o);
473 ijk[2] = ex_to<Index>(o);
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]);
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;
520 if(flags & status_flags::evaluated)
return *
this;
524 if(c01 && c23)
return this->hold();
528 else return this->hold();
537 c.s <<
"f(" <<
ijkl[0] <<
"," <<
ijkl[1] <<
"," <<
ijkl[2] <<
"," <<
ijkl[3] <<
")";
541 c <<
"SUNF[" <<
ijkl[0] <<
"," <<
ijkl[1] <<
"," <<
ijkl[2] <<
"," <<
ijkl[3] <<
"]";
550 c <<
"f4(" <<
ijkl[0] <<
"," <<
ijkl[1] <<
"," <<
ijkl[2] <<
"," <<
ijkl[3] <<
")";
558 ensure_if_modifiable();
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]);
579 inherited::read_archive(n);
582 ijkl[0] = ex_to<Index>(o);
584 ijkl[1] = ex_to<Index>(o);
586 ijkl[2] = ex_to<Index>(o);
588 ijkl[3] = ex_to<Index>(o);
606 }
else if(is_a<mul>(e) || is_a<ncmul>(e)) {
610 if(!is_a<add>(rei)) rei = lst{ rei };
613 for(
auto oi : ores)
for(
auto ri : rei) res.append(oi * ri);
616 for(
auto ri : res) ret += ri;
618 }
else return e.map(self);
620 return inner_expand(expr);
630 if(!expr_in.has(GMat(
w1,
w2,
w3)))
return expr_in;
633 expr = expr.subs(GMat(
w1,
w2,
w2)==TR(
w1));
637 for(
auto cv : cv_lst) {
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;
643 }
else if(!is_a<mul>(e))
throw Error(
"GMatContract: collect error: " +
ex2str(e));
646 for(
auto item : e) mats.append(item);
648 std::map<ex,int,ex_is_less> to_map, from_map;
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});
658 if(!item.match(GMat(
w1,
w2,
w3))) {
659 cout <<
"item in GMatContract: " << item << endl;
660 throw Error(
"GMatContract faild!");
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));
665 mats2.append(
GMatT(item));
666 for(
int j=0; j<mats.nops(); j++)
if(j!=i) mats2.append(mats.op(j));
670 mats_idx.remove_all();
673 to_map[item.op(1)] = i+10;
674 from_map[item.op(2)] = i+10;
680 bool checked =
false;
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();
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));
695 if(to_map[item.op(1)]==0 && from_map[item.op(2)]==0) {
696 to_map[item.op(1)] = ii+10;
697 from_map[item.op(2)] = ii+10;
698 }
else if(to_map[item.op(2)]==0 && from_map[item.op(1)]==0) {
699 to_map[item.op(2)] = ii+10;
700 from_map[item.op(1)] = ii+10;
702 auto li =
get_op(mats, ii, 1);
703 auto ri =
get_op(mats, ii, 2);
707 throw Error(
"GMatContract: index conflict (2).");
710 if(mats_idx2.nops()<1)
break;
711 mats_idx = mats_idx2;
716 while(todo.size()>0) {
717 int c = *(todo.begin());
719 ex curMat = mats.op(c).op(0);
720 auto li=mats.op(c).op(1);
721 auto ri=mats.op(c).op(2);
723 if(li.is_equal(ri)) {
724 retMat *= TR(curMat);
728 int fi = from_map[li];
730 retMat *= GMat(curMat, li, ri);
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);
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);
751 expr += cv.op(0) * retMat;
760 auto found = cache.find(ei);
761 if(found!=cache.end())
return found->second;
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);
769 }
else return e.map(self);
774 if(idx_lst.nops()==0)
return ei;
779 auto v =
form(cv.op(1));
782 if(!is_a<mul>(v)) v = lst{ v };
786 if(!is_a<Pair>(vi)) r *= vi;
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);
791 res += r * c.subs(repl);
801 if(e.match(GMat(
w1,
w2,
w3))) {
805 for(
auto item : e0) {
806 if(item.return_type()==return_types::commutative) c *= item;
809 cout <<
"c=" << c <<
", " <<
"v=" << v << endl;
810 throw Error(
"GMatOut: v != 1");
815 if(v.is_equal(1)) v =
GAS(1);
816 return c *
GMatOut(GMat(v, e.op(1), e.op(2)));
818 }
else return e.map(self);
820 return inner_out(expr_in);
827 auto found = cache.find(key);
828 if(found!=cache.end())
return found->second;
831 if(!e.has(GMat(
w1,
w2,
w3)))
return e;
832 else if(e.match(GMat(
w1,
w2,
w3))) {
836 for(
auto item : e0) res +=
GMatExpand(GMat(item, e.op(1), e.op(2)));
838 }
else if(is_a<mul>(e0)) {
840 for(
auto item : e0) {
841 if(item.return_type()==return_types::commutative) c *= item;
844 cout <<
"c=" << c <<
", " <<
"v=" << v << endl;
845 throw Error(
"GMatExpand: v != 1");
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)) {
855 for(
auto item : e0) {
862 if(!is_a<add>(ncL)) ncL = lst{ ncL };
864 if(!is_a<add>(ncR)) ncR = lst{ ncR };
866 for(
auto iL : ncL)
for(
auto iR : ncR) res += iL * iR;
870 if(!is_a<add>(rs)) rs = lst{ rs };
871 for(
auto item : rs) {
873 if(is_a<mul>(item)) {
874 if(item.nops()==1)
throw Error(
"GMatExpand: item.nops == 1");
875 for(
auto it : item) {
876 if(it.return_type()==return_types::commutative) c *= it;
878 if(!v.is_equal(1))
throw Error(
"GMatExpand: v != 1");
894 if(last==vi && is_a<DGamma>(vi)) {
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;
900 else throw Error(
"GMatExpand: only GAS(i/p/1/5) supported.");
903 if(last!=
GAS(1) && !last.is_equal(1)) vv = vv * last;
908 if(!last.is_equal(1) && last!=
GAS(1)) v = vv * last;
913 if(v.is_equal(1)) v =
GAS(1);
914 res += c * GMat(v, e.op(1), e.op(2));
918 }
else return e.map(self);
920 ex res = inner_expand(expr_in);
925 ex
GMatShift(
const ex & expr,
const ex & g,
bool to_right) {
926 if(!expr.has(g))
return expr;
928 ex key = lst{expr, g, to_right ? 1 : 0};
930 auto found = cache.find(key);
931 if(found!=cache.end())
return found->second;
934 if(!e.has(g) || !e.has(GMat(
w1,
w2,
w3)))
return e;
935 else if(e.match(GMat(
w1,
w2,
w3))) {
937 if(!is_a<ncmul>(eg)) eg = lst{ eg };
940 for(
int i=0; i<eg.nops()-1; i++)
if(eg.op(i)==g) { gi = i; break; }
943 for(
int i=eg.nops()-1; i>0; i--)
if(eg.op(i)==g) { gi = i;
break; }
946 int gj = gi + ( to_right ? 1 : -1 );
947 ex rem = 1, rem2 = 1;
948 for(
int i=0; i<eg.nops(); i++) {
954 if(to_right) rem2 *= eg.op(gj)*eg.op(gi);
955 else rem2 *= eg.op(gi)*eg.op(gj);
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));
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.");
970 if(rem.is_equal(1)) rem =
GAS(1);
971 if(rem2.is_equal(1)) rem2 =
GAS(1);
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));
976 res = res - GMat(rem2, e.op(1), e.op(2));
978 }
else return e.map(self);
982 res = inner_shift(res);
996 lst shift_1st_to_right(
const ex & e) {
997 if(!is_a<ncmul>(e))
throw Error(
"input is not a ncmul.");
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 }};
1005 cout << endl << e << endl;
1006 throw Error(
"shift_1st_to_right: only GAS(i/p/5) supproted.");
1011 for(
int i=2; i<n; i++) rem *= e.op(i);
1012 lst res = shift_1st_to_right(e.op(0)*rem);
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);
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))) {
1021 throw Error(
"shift_12_right: not a Vector or Index");
1023 if(!is_a<Index>(e.op(1).op(0)) && !is_a<Vector>(e.op(1).op(0))) {
1025 throw Error(
"shift_12_right: not a Vector or Index");
1027 res.append(lst{ 2*
SP(e.op(0).op(0), e.op(1).op(0)), rem });
1035 auto found = cache.find(expr);
1036 if(found!=cache.end())
return found->second;
1040 if(e.match(GMat(
w1,
w2,
w3))) {
1042 if(!is_a<ncmul>(eg)) eg = lst{ eg };
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))) {
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);
1061 lst cvs = shift_1st_to_right(exM);
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)));
1071 }
else return e.map(self);
1075 res = inner_shift(res);
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))) {
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.");
1094 if(eg.is_equal(1))
throw Error(
"GMatECC:: eg is 1, NOT expected.");
1096 if(!is_a<ncmul>(eg)) eg = lst{ eg };
1098 for(
int i=0; i<eg.nops(); i++)
if(eg.op(i)==
DGamma::C) { ci = i;
break; }
1099 if(ci==-1)
return e;
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;
1104 for(
int i=cj+1; i<eg.nops(); i++)
if(eg.op(i)==
DGamma::C) { cnt++; }
1106 for(
int i=0; i<ci; i++) res *= eg.op(i);
1108 for(
int i=ci+1; i<cj; i++) m *= eg.op(i);
1109 if(sign<0) cc *= -1;
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);
1118 res = inner_ecc(res);
1119 if(res.has(TR(
w))) {
1121 if(!e.has(TR(
w)))
return e;
1122 else if(e.match(TR(
w))) {
1125 for(
auto item :
gs) {
1126 auto gi = ex_to<DGamma>(item);
1127 if(!item.is_equal(
DGamma::C) && !gi.isTr) {
1134 }
else return e.map(self);
1142 if(!e.has(GMat(
w1,
w2,
w3)))
return e;
1143 else if(e.match(GMat(
w1,
w2,
w3))) {
1145 }
else return e.map(self);
1148 res = inner_transpose(res);
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 <<
"]";
1159 REGISTER_FUNCTION(GMat, do_not_evalf_params().print_func<FCFormat>(&GMat_fc_print).conjugate_func(mat_conj).set_return_type(return_types::commutative))
1161 bool IsZero(
const ex & e) {
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);
1169 for(
int i=0; i<5; i++) {
1171 for(
auto item :
vs) {
1172 nsubs[item] = ex(1)/n_nth_prime(n);
1175 ex ret = e.subs(nsubs);
1176 if(!normal(e).is_zero())
return false;
1180 return ret.is_zero();
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) {
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) {
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);
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);
1217 return ca_map(expr);
1224 res = res.subs(
NF==
CA);
1250 if(!is_a<mul>(res)) res = lst{ res };
1252 for(
auto item : res) {
1253 if(item.has(
NF)) c *= item;
1257 int deg = c.degree(
NF);
1258 int ldeg = -c.ldegree(
NF);
1259 if(ldeg>deg) deg = ldeg;
1263 for(
int i=0; i<=deg; i++) {
1265 eqn -= xi * pow(
NF,i) * pow((
NF*
NF-1)/(2*
NF), deg-i);
1269 int nH = eqn.degree(
NF);
1270 int nL = eqn.ldegree(
NF);
1272 for(
int i=
nL; i<=
nH; i++) {
1273 ex cc = eqn.coeff(
NF, i);
1274 if(cc.is_zero())
continue;
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!");
1284 for(
int i=0; i<=deg; i++) c += vars.op(i).subs(sol) * pow(
CA,i) * pow(
CF, deg-i);
1289 ex
DoColor(
const ex & expr,
const ex & pref,
int method) {
1292 for(
auto const & cv : cvs) {
#define IMPLEMENT_HAS(classname)
#define DEFAULT_CTOR(classname)
#define IMPLEMENT_ALL(classname)
static lst all(const ex &e)
static bool has(const ex &e)
class used to wrap error message
int compare_same_type(const GiNaC::basic &other) const override
const char * class_name() const override
static bool hasv(const ex &e)
void archive(archive_node &n) const override
bool is_equal_same_type(const basic &other) const override
static GiNaC::registered_class_info & get_class_info_static()
void print(const print_context &c, unsigned level=0) const
void read_archive(const archive_node &n) override
static bool hasc(const ex &e)
ex derivative(const symbol &s) const override
const GiNaC::registered_class_info & get_class_info() const override
Pair operator()(const Index &i)
void accept(GiNaC::visitor &v) const override
Index * duplicate() const override
static bool has(const ex &e)
class to wrap map_function of GiNaC
static bool has(const ex &e)
void archive(archive_node &n) const override
save to archvie
static GiNaC::registered_class_info & get_class_info_static()
void form_print(const FormFormat &c, unsigned level=0) const
print the Form Format
void accept(GiNaC::visitor &v) const override
const GiNaC::registered_class_info & get_class_info() const override
size_t nops() const override
void fc_print(const FCFormat &c, unsigned level=0) const
void print(const print_dflt &c, unsigned level=0) const
normal priint
ex op(size_t i) const override
ex & let_op(size_t i) override
SUNF4 * duplicate() const override
ex eval() const override
automatical evaluation of SUNF4
ex derivative(const symbol &s) const override
set derivative of SUNF4 to 0
bool is_equal_same_type(const basic &other) const override
const char * class_name() const override
void read_archive(const archive_node &n) override
read from archive
int compare_same_type(const GiNaC::basic &other) const override
void archive(archive_node &n) const override
void print(const print_dflt &c, unsigned level=0) const
const char * class_name() const override
int compare_same_type(const GiNaC::basic &other) const override
ex derivative(const symbol &s) const override
set derivative of SUNF to 0
bool is_equal_same_type(const basic &other) const override
void form_print(const FormFormat &c, unsigned level=0) const
static GiNaC::registered_class_info & get_class_info_static()
ex & let_op(size_t i) override
SUNF * duplicate() const override
const GiNaC::registered_class_info & get_class_info() const override
void read_archive(const archive_node &n) override
size_t nops() const override
void accept(GiNaC::visitor &v) const override
void fc_print(const FCFormat &c, unsigned level=0) const
ex op(size_t i) const override
static GiNaC::registered_class_info & get_class_info_static()
bool is_equal_same_type(const basic &other) const override
ex & let_op(size_t i) override
void archive(archive_node &n) const override
size_t nops() const override
void form_print(const FormFormat &c, unsigned level=0) const
ex op(size_t i) const override
const GiNaC::registered_class_info & get_class_info() const override
void fc_print(const FCFormat &c, unsigned level=0) const
void print(const print_dflt &c, unsigned level=0) const
int compare_same_type(const GiNaC::basic &other) const override
void read_archive(const archive_node &n) override
SUNT * duplicate() const override
const char * class_name() const override
void accept(GiNaC::visitor &v) const override
ex derivative(const symbol &s) const override
ex conjugate() const override
class extended to GiNaC symbol class, represent a positive symbol
Vector * duplicate() const override
bool is_equal_same_type(const basic &other) const override
Pair operator()(const Vector &p)
const char * class_name() const override
ex derivative(const symbol &s) const override
void print(const print_context &c, unsigned level=0) const
static GiNaC::registered_class_info & get_class_info_static()
int compare_same_type(const GiNaC::basic &other) const override
void read_archive(const archive_node &n) override
void accept(GiNaC::visitor &v) const override
void archive(archive_node &n) const override
const GiNaC::registered_class_info & get_class_info() const override
ex exfactor(const ex &expr_in, int opt)
factorize a expression
bool ex_less(const ex &a, const ex &b)
ex DoColor(const ex &expr, const ex &pref, int method)
void let_op(ex &ex_in, int index1, int index2, const ex item)
update index1-th.index2-th of expression with item
ex GMatShift(const ex &expr, const ex &g, bool to_right)
ex collect_ex(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica
ex GMatSimplify(const ex &expr)
ex GAS(const ex &expr, unsigned rl)
function similar to GAD/GSD in FeynClac
ex GMatExpand(const ex &expr_in)
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)
ex Contract(const ex &ei)
ex get_op(const ex ex_in, int index1, int index2)
return index1-th.index2-th of expression
ex GMatECC(const ex &expr, int sign)
ex GMatOut(const ex &expr_in)
ex exnormal(const ex &expr, int opt)
normalize a expression
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}, ... }
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
ex ncmul_expand(const ex &expr)
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)
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)
ex ca_neg_pow_sub(const ex &expr)