Advanced Computing Platform for Theoretical Physics

Commit 32cf7dfb authored by Feng Feng's avatar Feng Feng
Browse files

need to remap the dummy index in TIR

parent 1de3238f
......@@ -104,88 +104,100 @@ namespace HepLib::FC {
for(int j=0; j<er; j++) bi *= SP(item.op(j), is.op(j));
for(int j=er; j<visn; j=j+2) bi *= SP(is.op(j), is.op(j+1));
bi = bi.symmetrize(is);
bi = bi.subs(SP_map);
if(!is_zero(bi)) bis.append(bi);
bis.append(bi);
}}
ex res = 0;
int n = bis.nops();
if(n>0) {
lst mat;
for(auto bi : bis) {
for(auto bj : bis) mat.append(bi*bj);
lst mat;
for(auto bi : bis) {
for(auto bj : bis) mat.append(bi*bj);
}
for(auto bj : bis) mat.append(eqL*bj);
// we need to remap the dummy index, to avoid SP_map set the index object to 0
if(true) {
lst isu = is;
isu.sort();
isu.unique();
exmap i2u, u2i;
int c=0;
for(auto item : isu) {
auto ii = Index("_TIR"+to_string(c++));
i2u[item] = ii;
u2i[ii] = item;
}
for(auto bj : bis) mat.append(eqL*bj);
mat = ex_to<lst>(form(mat));
mat = ex_to<lst>(subs(mat,i2u));
mat = ex_to<lst>(form(mat).subs(u2i));
}
lst rep_vs;
ex tree = mat;
for(const_preorder_iterator i = tree.preorder_begin(); i != tree.preorder_end(); ++i) {
auto e = (*i);
if(is_a<symbol>(e) || is_a<Pair>(e) || is_a<Eps>(e)) {
rep_vs.append(e);
}
}
rep_vs.sort();
rep_vs.unique();
sort_lst(rep_vs);
exmap v2f;
symtab st;
int fvi = 0;
for(auto vi : rep_vs) {
auto name = "v" + to_string(fvi);
v2f[vi] = Symbol(name);
st[name] = vi;
fvi++;
lst rep_vs;
ex tree = mat;
for(const_preorder_iterator i = tree.preorder_begin(); i != tree.preorder_end(); ++i) {
auto e = (*i);
if(is_a<symbol>(e) || is_a<Pair>(e) || is_a<Eps>(e)) {
rep_vs.append(e);
}
}
rep_vs.sort();
rep_vs.unique();
sort_lst(rep_vs);
exmap v2f;
symtab st;
int fvi = 0;
for(auto vi : rep_vs) {
auto name = "v" + to_string(fvi);
v2f[vi] = Symbol(name);
st[name] = vi;
fvi++;
}
stringstream ss;
for(int i=0; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
Fermat fermat;
fermat.Init();
fermat.Execute(ss.str());
ss.clear();
ss.str("");
ss << "Array m[" << n << "," << n+1 << "];" << endl;
fermat.Execute(ss.str());
ss.clear();
ss.str("");
ss << "[m]:=[(";
for(auto item : mat) {
ss << item.subs(v2f) << ",";
}
ss << ")];" << endl;
ss << "Redrowech([m]);" << endl;
auto tmp = ss.str();
stringstream ss;
for(int i=0; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
Fermat fermat;
fermat.Init();
fermat.Execute(ss.str());
ss.clear();
ss.str("");
ss << "Array m[" << n << "," << n+1 << "];" << endl;
fermat.Execute(ss.str());
ss.clear();
ss.str("");
ss << "[m]:=[(";
for(auto item : mat) {
ss << item.subs(v2f) << ",";
}
ss << ")];" << endl;
ss << "Redrowech([m]);" << endl;
auto tmp = ss.str();
string_replace_all(tmp,",)]",")]");
fermat.Execute(tmp);
ss.clear();
ss.str("");
string_replace_all(tmp,",)]",")]");
fermat.Execute(tmp);
ss.clear();
ss.str("");
ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
ss << "![m" << endl;
auto ostr = fermat.Execute(ss.str());
fermat.Exit();
// make sure last char is 0
if(ostr[ostr.length()-1]!='0') throw Error("TIR: last char is NOT 0.");
ostr = ostr.substr(0, ostr.length()-1);
string_trim(ostr);
ostr.erase(0, ostr.find(":=")+2);
string_replace_all(ostr, "[", "{");
string_replace_all(ostr, "]", "}");
Parser fp(st);
auto mat2 = fp.Read(ostr);
ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
ss << "![m" << endl;
auto ostr = fermat.Execute(ss.str());
fermat.Exit();
// make sure last char is 0
if(ostr[ostr.length()-1]!='0') throw Error("TIR: last char is NOT 0.");
ostr = ostr.substr(0, ostr.length()-1);
string_trim(ostr);
ostr.erase(0, ostr.find(":=")+2);
string_replace_all(ostr, "[", "{");
string_replace_all(ostr, "]", "}");
Parser fp(st);
auto mat2 = fp.Read(ostr);
for(int i=0; i<n; i++) {
if(is_zero(mat2.op(i).op(i)) && !is_zero(mat2.op(i).op(n))) throw Error("Zero Determinant in TIR.");
res += bis.op(i) * mat2.op(i).op(n);
}
ex res = 0;
for(int i=0; i<n; i++) {
if(is_zero(mat2.op(i).op(i)) && !is_zero(mat2.op(i).op(n))) throw Error("Zero Determinant in TIR.");
res += bis.op(i) * mat2.op(i).op(n);
}
res = res.subs(SP_map);
cache_map[map_key] = res;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment