Advanced Computing Platform for Theoretical Physics

Commit dda4a772 authored by Feng Feng's avatar Feng Feng
Browse files

remove some normal operations in Apart

parent c9556078
......@@ -725,7 +725,7 @@ namespace HepLib {
* @param exvec input exvector
* @return lst
*/
lst exvec2lst(const exvector & exvec) {
lst vec2lst(const exvector & exvec) {
lst ret;
for(auto item : exvec) ret.append(item);
return ret;
......@@ -1931,7 +1931,7 @@ namespace HepLib {
* @param factor true for factorize on the denominator
* @return a list of { numer, denom }
*/
ex fermat_numer_denom(const ex & expr, bool factor) {
ex numer_denom_fermat(const ex & expr, bool factor) {
static map<pid_t, Fermat> fermat_map;
static int v_max = 0;
......@@ -1988,7 +1988,7 @@ namespace HepLib {
if(!is_a<mul>(expr_in)) expr_in = lst{expr_in};
for(auto item : expr_in) {
if(!is_a<add>(item)) item = lst{item};
if(fermat_use_array) ss << "Array m[" << item.nops() << "];" << endl;
if(fermat_using_array) ss << "Array m[" << item.nops() << "];" << endl;
else ss << "res:=0;" << endl;
fermat.Execute(ss.str());
ss.clear();
......@@ -1998,15 +1998,15 @@ namespace HepLib {
for(int i=0; i<item.nops(); i++) {
ex tt = item.op(i).subs(v2f);
nn_chk2 += tt.subs(nn_map);
if(fermat_use_array) ss << "m[" << (i+1) << "]:=";
if(fermat_using_array) ss << "m[" << (i+1) << "]:=";
else ss << "item:=";
ss << tt << ";" << endl;
if(!fermat_use_array) ss << "res:=res+item;" << endl;
if(!fermat_using_array) ss << "res:=res+item;" << endl;
fermat.Execute(ss.str());
ss.clear();
ss.str("");
}
if(fermat_use_array) {
if(fermat_using_array) {
ss << "res:=Sumup([m]);" << endl;
fermat.Execute(ss.str());
ss.clear();
......@@ -2021,9 +2021,9 @@ namespace HepLib {
ss.clear();
ss.str("");
// note the order, before exfactor (fermat_normal will be called again here)
// note the order,(normal_fermat will be called again in factor_form)
ss << "&(U=0);" << endl; // disable ugly printing
if(fermat_use_array) ss << "@(res,[m]);" << endl;
if(fermat_using_array) ss << "@(res,[m]);" << endl;
else ss << "@(res,item);" << endl;
ss << "&_G;" << endl;
fermat.Execute(ss.str());
......@@ -2045,7 +2045,7 @@ namespace HepLib {
Parser fp(st);
auto ret = fp.Read(ostr);
num *= ret.op(0);
if(factor) den *= exfactor(ret.op(1));
if(factor) den *= factor_form(ret.op(1));
else den *= ret.op(1);
}
//fermat.Exit();
......@@ -2067,7 +2067,7 @@ namespace HepLib {
* @param factor true for factorize on the denominator
* @return the normalized expression: numer/denom
*/
ex fermat_normal(const ex & expr, bool factor) {
ex normal_fermat(const ex & expr, bool factor) {
auto nd = fermat_numer_denom(expr, factor);
return nd.op(0)/nd.op(1);
}
......@@ -2087,15 +2087,37 @@ namespace HepLib {
/**
* @brief factorize a expression
* @param expr the input expression
* @param fm FactorMethod::FORM to use FORM, otherwise using GiNaC for factorization
* @param opt 1 to use FORM, otherwise using GiNaC for factorization
* @return factorized result
*/
ex exfactor(const ex & expr, FactorMethod fm) {
if(fm == FactorMethod::FORM) return form_factor(expr);
ex exfactor(const ex & expr, int opt) {
if(opt==1) return factor_form(expr);
else return ginac_factor(expr);
}
ex inner_form_factor(const ex & expr) {
/**
* @brief normalize a expression
* @param expr the input expression
* @param opt 1 to use Fermat, otherwise using GiNaC for normalization
* @return factorized result
*/
ex exnormal(const ex & expr, int opt) {
if(opt==1) return normal_fermat(expr);
else return normal(expr);
}
/**
* @brief num_den a expression
* @param expr the input expression
* @param opt 1 to use Fermat, otherwise using GiNaC for numer_denom
* @return lst of { num, den }
*/
ex exnd(const ex & expr, int opt) {
if(opt==1) return numer_denom_fermat(expr);
else return numer_denom(expr);
}
ex inner_factor_form(const ex & expr) {
static map<pid_t, Form> form_map;
auto pid = getpid();
if((form_map.find(pid)==form_map.end())) { // init section
......@@ -2155,10 +2177,10 @@ namespace HepLib {
* @param expr the input expression
* @return factorized result
*/
ex form_factor(const ex & expr) {
ex factor_form(const ex & expr) {
auto num_den = fermat_numer_denom(expr);
if(is_zero(num_den.op(1)-1)) return inner_form_factor(num_den.op(0));
return inner_form_factor(num_den.op(0))/inner_form_factor(num_den.op(1));
if(is_zero(num_den.op(1)-1)) return inner_factor_form(num_den.op(0));
return inner_factor_form(num_den.op(0))/inner_factor_form(num_den.op(1));
}
//-----------------------------------------------------------
......
......@@ -54,23 +54,23 @@ namespace HepLib {
/*-----------------------------------------------------*/
// Terminal Color
/*-----------------------------------------------------*/
#define RESET "\033[0m"
#define BLACK "\033[30m" /* Black */
#define RED "\033[31m" /* Red */
#define GREEN "\033[32m" /* Green */
#define YELLOW "\033[33m" /* Yellow */
#define BLUE "\033[34m" /* Blue */
#define MAGENTA "\033[35m" /* Magenta */
#define CYAN "\033[36m" /* Cyan */
#define WHITE "\033[37m" /* White */
#define BOLDBLACK "\033[1m\033[30m" /* Bold Black */
#define BOLDRED "\033[1m\033[31m" /* Bold Red */
#define BOLDGREEN "\033[1m\033[32m" /* Bold Green */
#define BOLDYELLOW "\033[1m\033[33m" /* Bold Yellow */
#define BOLDBLUE "\033[1m\033[34m" /* Bold Blue */
#define BOLDMAGENTA "\033[1m\033[35m" /* Bold Magenta */
#define BOLDCYAN "\033[1m\033[36m" /* Bold Cyan */
#define BOLDWHITE "\033[1m\033[37m" /* Bold White */
#define RESET "\033[0m"
#define BLACK "\033[30m"
#define RED "\033[31m"
#define GREEN "\033[32m"
#define YELLOW "\033[33m"
#define BLUE "\033[34m"
#define MAGENTA "\033[35m"
#define CYAN "\033[36m"
#define WHITE "\033[37m"
#define BOLDBLACK "\033[1m\033[30m"
#define BOLDRED "\033[1m\033[31m"
#define BOLDGREEN "\033[1m\033[32m"
#define BOLDYELLOW "\033[1m\033[33m"
#define BOLDBLUE "\033[1m\033[34m"
#define BOLDMAGENTA "\033[1m\033[35m"
#define BOLDCYAN "\033[1m\033[36m"
#define BOLDWHITE "\033[1m\033[37m"
/**
* @brief class extended to GiNaC symbol class, represent a positive symbol
......@@ -268,7 +268,9 @@ namespace HepLib {
ex garRead(const string &garfn, const char* key);
ex garRead(const string &garfn);
void garWrite(const string &garfn, const map<string, ex> &resMap);
inline void garWrite(const map<string, ex> &resMap, const string &garfn) { garWrite(garfn,resMap); }
void garWrite(const string &garfn, const ex & res);
inline void garWrite(const ex & res, const string &garfn) { garWrite(garfn,res); }
ex str2ex(const string &expr, symtab stab);
ex str2ex(const string &expr);
lst str2lst(const string &expr, symtab stab);
......@@ -283,7 +285,7 @@ namespace HepLib {
string ex2str(const ex &expr);
ex q2ex(__float128);
__float128 ex2q(ex);
lst exvec2lst(const exvector & exvec);
lst vec2lst(const exvector & exvec);
exvector lst2vec(const lst & alst);
lst add2lst(const ex & expr);
lst mul2lst(const ex & expr);
......@@ -344,14 +346,19 @@ namespace HepLib {
ex mma_diff(ex const expr, ex const xp, unsigned nth=1, bool expand=false);
extern bool fermat_use_array;
ex fermat_numer_denom(const ex & expr, bool factor=false);
ex fermat_normal(const ex & expr, bool factor=false);
extern bool fermat_using_array;
ex numer_denom_fermat(const ex & expr, bool factor=false);
inline ex fermat_numer_denom(const ex & expr, bool factor=false) { return numer_denom_fermat(expr,factor); }
ex form_factor(const ex & expr);
ex normal_fermat(const ex & expr, bool factor=false);
inline ex fermat_normal(const ex & expr, bool factor=false) { return normal_fermat(expr,factor); }
enum FactorMethod { GiNaC, FORM };
ex exfactor(const ex & expr, FactorMethod fm = FORM);
ex factor_form(const ex & expr);
inline ex form_factor(const ex & expr) { return factor_form(expr); }
ex exfactor(const ex & expr, int opt = 1);
ex exnormal(const ex & expr, int opt = 1);
ex exnd(const ex & expr, int opt = 1);
ex collect_factors(const ex & expr);
......
......@@ -97,11 +97,11 @@ expr expr::expand() {
}
expr expr::normal() {
return expr(GiNaC::ex(HepLib::fermat_normal(_expr)));
return expr(GiNaC::ex(HepLib::normal_fermat(_expr)));
}
expr expr::factor() {
return expr(GiNaC::ex(HepLib::form_factor(_expr)));
return expr(GiNaC::ex(HepLib::factor_form(_expr)));
}
expr expr::series(const expr &s, int o) {
......@@ -114,11 +114,11 @@ expr expand(const expr &e) {
}
expr normal(const expr &e) {
return expr(HepLib::fermat_normal(e._expr));
return expr(HepLib::normal_fermat(e._expr));
}
expr factor(const expr &e) {
return expr(HepLib::form_factor(e._expr));
return expr(HepLib::factor_form(e._expr));
}
expr subs(const expr &e, const std::vector<expr> &ev) {
......@@ -256,8 +256,8 @@ void Integral::Replacements(const std::vector<expr> &ev, int lt) {
void Integral::Evaluate() {
if(isX) {
HepLib::SD::XIntegrand xint;
xint.Functions = HepLib::exvec2lst(_Propagators_Functions);
xint.Exponents = HepLib::exvec2lst(_Exponents);
xint.Functions = HepLib::vec2lst(_Propagators_Functions);
xint.Exponents = HepLib::vec2lst(_Exponents);
HepLib::SD::SecDec secdec;
HepLib::Verbose = verb;
......@@ -267,10 +267,10 @@ void Integral::Evaluate() {
_Result = secdec.ResultError;
} else {
HepLib::SD::FeynmanParameter fp;
fp.Propagators = HepLib::exvec2lst(_Propagators_Functions);
fp.Exponents = HepLib::exvec2lst(_Exponents);
fp.LoopMomenta = HepLib::exvec2lst(_Loops);
fp.tLoopMomenta = HepLib::exvec2lst(_tLoops);
fp.Propagators = HepLib::vec2lst(_Propagators_Functions);
fp.Exponents = HepLib::vec2lst(_Exponents);
fp.LoopMomenta = HepLib::vec2lst(_Loops);
fp.tLoopMomenta = HepLib::vec2lst(_tLoops);
for(auto item : _lReplacements) fp.lReplacements[item.op(0)] = item.op(1);
for(auto item : _lReplacements) fp.tReplacements[item.op(0)] = item.op(1);
......
......@@ -122,7 +122,7 @@ namespace HepLib {
// null vector
static exmap null_cache;
if(null_cache.find(sub_matrix(mat,0,nrow,0,ncol))==null_cache.end()){
if(null_cache.find(sub_matrix(mat,0,nrow,0,ncol))==null_cache.end()) {
if(Apart_using_fermat) {
static map<pid_t, Fermat> fermat_map;
static int v_max = 0;
......@@ -177,7 +177,7 @@ namespace HepLib {
ss << "[m]:=[(";
for(int c=0; c<ncol; c++) {
for(int r=0; r<nrow; r++) {
ss << mat(r,c).subs(v2f) << ",";
ss << mat(r,c).subs(iEpsilon==0).subs(v2f) << ",";
}
}
for(int r=0; r<nrow; r++) ss << "0,";
......@@ -196,7 +196,7 @@ namespace HepLib {
ss.str("");
//fermat.Exit();
// note the order, before exfactor (fermat_normal will be called again here)
// note the order, before exfactor (normal_fermat will be called again here)
ss << "&(U=0);" << endl; // disable ugly printing
ss << "@([m]);" << endl;
ss << "&_G;" << endl;
......@@ -237,7 +237,7 @@ namespace HepLib {
}
// Solve M*V = 0
matrix zero(nrow, 1);
matrix s = ex_to<matrix>(sub_matrix(mat,0,nrow,0,ncol)).solve(v,zero);
matrix s = ex_to<matrix>(sub_matrix(mat,0,nrow,0,ncol).subs(iEpsilon==0)).solve(v,zero);
for(int r=0; r<ncol; r++) null_vec.append(s(r,0).subs(sRepl));
}
null_cache[sub_matrix(mat,0,nrow,0,ncol)] = null_vec;
......@@ -421,7 +421,7 @@ namespace HepLib {
}
nrepl[iEpsilon]=0;
ex chk = ApartIR2ex(subs(res,nrepl))-subs(expr_ino,nrepl);
chk = normal(chk);
chk = normal_fermat(chk);
if(!is_zero(chk)) throw Error("Apart@1 random check Failed.");
return res;
}
......@@ -443,7 +443,7 @@ namespace HepLib {
ex expr = expr_in.subs(map1);
if(!is_a<mul>(expr)) expr = lst{expr};
// check only, try fermat_normal if faild
// check only, try normal_fermat if faild
bool ok = true;
for(auto item : expr) {
bool has_var=false;
......@@ -462,12 +462,13 @@ namespace HepLib {
}
if(!ok) {
expr = expr_in.subs(map1);
expr = fermat_normal(expr,true); // need option factor=true
expr = normal_fermat(expr,true); // need option factor=true, factor denominator
if(!is_a<mul>(expr)) expr = lst{expr};
}
lst pnlst;
ex pref = 1;
map<ex,int,ex_is_less> count_ip;
for(auto item : expr) {
bool has_var=false;
for(auto v : vars) {
......@@ -519,9 +520,39 @@ namespace HepLib {
break;
}
}
ex key = pc.subs(iEpsilon==0).expand();
count_ip[key] = count_ip[key]+1;
pnlst.append(lst{ pc, nc });
} else pref *= item;
}
bool needs_again = false;
for(auto kv : count_ip) {
if(kv.second>1) { needs_again = true; break; }
}
if(needs_again) {
exmap imap;
for(auto pn : pnlst) {
auto key = pn.op(0);
if(key.has(iEpsilon)) imap[key.subs(iEpsilon==0)] = key;
}
exmap p2n;
for(auto pn : pnlst) {
auto key = pn.op(0);
if(!key.has(iEpsilon)) key = key.subs(imap);
p2n[key] = p2n[key] + pn.op(1);
}
pnlst.remove_all();
for(auto kv : p2n) {
if(is_zero(kv.second)) continue;
auto k = kv.first;
auto v = kv.second;
if(is_a<numeric>(v) && ex_to<numeric>(v)>0) {
k = k.subs(iEpsilon==0);
}
pnlst.append(lst{k, v});
}
}
sort_lst(pnlst);
if(pnlst.nops()==0) return pref * ApartIR(1,vars_in);
......@@ -532,12 +563,12 @@ namespace HepLib {
matrix mat(nrow+2, ncol);
for(int c=0; c<ncol; c++) {
ex pn = pnlst.op(c);
mat(nrow+1,c) = pn.op(1);
mat(nrow+1,c) = normal(pn.op(1));
auto tmp = pn.op(0);
for(int r=0; r<nrow; r++) {
mat(r,c) = tmp.coeff(vars.op(r));
mat(r,c) = normal_fermat(tmp.coeff(vars.op(r)));
}
mat(nrow,c) = tmp.subs(vars0);
mat(nrow,c) = normal_fermat(tmp.subs(vars0));
}
ex ret = Apart(mat);
......@@ -664,6 +695,8 @@ namespace HepLib {
std::function<lst(const Base &, const ex &)> uf) {
if(loops_exts.nops()<2) throw Error("loops_exts size() < 2;");
lst lmom = ex_to<lst>(loops_exts.op(0));
lst emom = ex_to<lst>(loops_exts.op(1));
string wdir = to_string(getpid());
if(IBPmethod==1) wdir = wdir + "_FIRE";
......@@ -678,13 +711,9 @@ namespace HepLib {
}
auto air_intg =
GiNaC_Parallel(air_vec.size(), [aparted,air_vec,loops_exts] (int idx) {
GiNaC_Parallel(air_vec.size(), [aparted,air_vec,lmom,emom] (int idx) {
auto air = air_vec[idx];
if(!aparted) {
lst lmom = ex_to<lst>(loops_exts.op(0));
lst emom = ex_to<lst>(loops_exts.op(1));
air = Apart(air,lmom,emom);
}
if(!aparted) air = Apart(air,lmom,emom);
air = air.subs(SP_map);
air = ApartIRC(air);
exset intg;
......@@ -724,14 +753,8 @@ namespace HepLib {
for(auto kv : sps) repls.append(kv.first == kv.second);
lst loops, exts;
for(auto li : loops_exts.op(0)) {
if(is_a<Vector>(li)) loops.append(ex_to<Vector>(li).name);
else loops.append(li);
}
for(auto li : loops_exts.op(1)) {
if(is_a<Vector>(li)) exts.append(ex_to<Vector>(li).name);
else exts.append(li);
}
for(auto li : lmom) loops.append(ex_to<Vector>(li).name);
for(auto li : emom) exts.append(ex_to<Vector>(li).name);
exmap IR2F;
std::map<ex, IBP::Base*, ex_is_less> p2IBP;
......@@ -753,7 +776,7 @@ namespace HepLib {
if(cut_props.nops()>0) {
ex cuts = cut_props;
cuts = cuts.subs(SP_map);
for(auto cut : cut_props) pns.prepend(lst{ SP2sp(cut), 1 });
for(auto cut : cuts) pns.prepend(lst{ SP2sp(cut), 1 });
vars.append(cuts);
}
......@@ -778,7 +801,7 @@ namespace HepLib {
ibp->Internal = loops;
ibp->External = exts;
ibp->Replacements = repls;
ibp->Pairs = ex_to<lst>(SP2sp(Pair::all(vars)));
//ibp->Pairs = ex_to<lst>(SP2sp(Pair::all(vars)));
ibp->WorkingDir = wdir;
ibp->ProblemNumber = pn++;
if(cut_props.nops()>0) {
......@@ -839,7 +862,7 @@ namespace HepLib {
air = _F2ex(air);
auto cv_lst = mma_collect_lst(air, F(w1,w2));
air = 0;
for(auto cv : cv_lst) air += cv.op(1) * fermat_normal(cv.op(0));
for(auto cv : cv_lst) air += cv.op(1) * (cv.op(0));
return air;
}, "A2F");
......@@ -883,7 +906,7 @@ namespace HepLib {
rr = subs_naive(rr,mi_rules.first);
auto cv_lst = mma_collect_lst(rr, F(w1,w2));
rr = 0;
for(auto cv : cv_lst) rr += cv.op(1) * fermat_normal(cv.op(0));
for(auto cv : cv_lst) rr += cv.op(1) * (cv.op(0));
rules.append(item.op(0)==rr);
}
return rules;
......@@ -904,7 +927,7 @@ namespace HepLib {
air = _F2ex(air);
auto cv_lst = mma_collect_lst(air, F(w1,w2));
air = 0;
for(auto cv : cv_lst) air += cv.op(1) * fermat_normal(cv.op(0));
for(auto cv : cv_lst) air += cv.op(1) * (cv.op(0));
return air;
}, "A2F");
......@@ -912,5 +935,261 @@ namespace HepLib {
for(int i=0; i<air_vec.size(); i++) air_vec[i] = air_res[i];
}
/**
* @brief perform IBP reduction on the Aparted input
* @param IBPmethod ibp method used, 0-No IBP, 1-FIRE, 2-KIRA
* @param air_vec vector contains aparted input, ApartIRC will be call internally
* @param aip AIOption for ApartIBP input
* @return nothing returned, the input air_vec will be updated
*/
void ApartIBP(int IBPmethod, exvector &air_vec, AIOption aip) {
string wdir = to_string(getpid());
if(IBPmethod==1) wdir = wdir + "_FIRE";
else if(IBPmethod==2) wdir = wdir + "_KIRA";
bool aparted = false;
for(auto air : air_vec) {
if(air.has(ApartIR(w1,w2))) {
aparted = true;
break;
}
}
if(!aparted) {
exset vset;
lst lmom = ex_to<lst>(aip.AInternal);
lst emom = ex_to<lst>(aip.AExternal);
for(auto &air : air_vec) {
air = mma_collect(air,lmom,false,true);
find(air,coVF(w),vset);
}
exvector vvec;
for(auto item : vset) vvec.push_back(item);
auto ret = GiNaC_Parallel(vvec.size(), [vvec,lmom,emom] (int idx) {
auto air = vvec[idx].op(0);
air = Apart(air,lmom,emom);
air = air.subs(SP_map);
air = ApartIRC(air);
return air;
}, "Apart");
exmap v2v;
for(int i=0; i<vvec.size(); i++) v2v[vvec[i]] = ret[i];
for(auto &air : air_vec) air = air.subs(v2v);
}
exset intg;
for(auto &air : air_vec) {
find(air, ApartIR(w1, w2), intg);
}
exvector intg_sort;
for(auto item : intg) intg_sort.push_back(item);
sort_vec(intg_sort);
for(auto sp : aip.CPairs) SP_map.erase(sp);
// from here, Vector will be replaced by its name Symbol
lst repls;
auto sps = sp_map();
for(auto kv : sps) repls.append(kv.first == kv.second);
lst loops, exts;
for(auto li : aip.IInternal) {
if(is_a<Vector>(li)) loops.append(ex_to<Vector>(li).name);
else loops.append(li);
}
for(auto li : aip.IExternal) {
if(is_a<Vector>(li)) exts.append(ex_to<Vector>(li).name);
else exts.append(li);
}
exmap IR2F;
std::map<ex, IBP::Base*, ex_is_less> p2IBP;
vector<IBP::Base*> ibp_vec;
int pn=1;
for(auto ir : intg_sort) {
auto mat = ex_to<matrix>(ir.op(0));
auto vars = ex_to<lst>(ir.op(1));
lst pns;
int nrow = mat.rows();
for(int c=0; c<mat.cols(); c++) {
ex pc = 0;
for(int r=0; r<nrow-2; r++) pc += mat(r,c) * vars.op(r);
pc += mat(nrow-2,c);
pc = SP2sp(pc);
pns.append(lst{ pc, ex(0)-mat(nrow-1,c) }); // note the convension
}
sort_lst(pns); // sort before cuts
if(aip.Cuts.nops()>0) {
ex cuts = aip.Cuts;
cuts = cuts.subs(SP_map);
for(auto cut : cuts) pns.prepend(lst{ SP2sp(cut), 1 });
}
lst props, ns;
for(auto item : pns) {
props.append(item.op(0));
ns.append(item.op(1));
}
if(p2IBP[props]==NULL) {
IBP::Base* ibp;
if(IBPmethod==0) ibp = new Base();
else if(IBPmethod==1) ibp = new FIRE();
else if(IBPmethod==2) ibp = new KIRA();
else {
ibp = new Base();
IBPmethod = 0;
}
p2IBP[props] = ibp;
ibp->Propagators = props;
ibp->Internal = loops;
ibp->External = exts;
ibp->Replacements = repls;
if(aip.IPairs.nops()>0) {
for(auto item : aip.IPairs) ibp->Pairs.append(SP2sp(item));
}
if(aip.I_Internal.nops()>0) {
for(auto item : aip.I_Internal) {
if(is_a<Vector>(item)) ibp->_Internal.append(ex_to<Vector>(item).name);
else ibp->_Internal.append(item);
}
}
ibp->WorkingDir = wdir;
ibp->ProblemNumber = pn++;
if(aip.Cuts.nops()>0) {
for(int i=0; i<aip.Cuts.nops(); i++) ibp->Cuts.append(i+1);
}
ibp_vec.push_back(ibp);
}
IBP::Base* ibp = p2IBP[props];
ibp->Integrals.append(ns);
IR2F[ir] = F(ibp->ProblemNumber, ns);
}
if(Verbose>0) cout << " \\--Total Ints/Pros: " << intg_sort.size() << "/" << ibp_vec.size() << " @ " << now(false) << endl;
vector<Base*> base_vec;
for(auto ibp : ibp_vec) base_vec.push_back(ibp);
auto rules_ints = FindRules(base_vec, false, aip.UF);
map<int,lst> pn_ints_map;
for(auto item : rules_ints.second) {
int pn = ex_to<numeric>(item.op(0)).to_int();
pn_ints_map[pn].append(item.op(1));
}
vector<IBP::Base*> ibp_vec_re;
int nints = 0;
for(auto pi : pn_ints_map) {
auto ibp = ibp_vec[pi.first-1];
ibp->Integrals = pi.second;
nints += ibp->Integrals.nops();
ibp_vec_re.push_back(ibp);
}
if(Verbose>0) cout << " \\--Refined Ints/Pros: " << nints << "/" << ibp_vec_re.size() << " @ " << now(false) << endl;
MapFunction _F2ex([ibp_vec](const ex &e, MapFunction &self)->ex {
if(!e.has(F(w1,w2))) return e;
else if(e.match(F(w1,w2))) {
int idx = ex_to<numeric>(e.op(0)).to_int();
auto pso = ex_to<lst>(ibp_vec[idx-1]->Propagators);
auto nso = ex_to<lst>(e.op(1));
lst ps, ns;
for(int i=0; i<pso.nops(); i++) {
if(nso.op(i).is_zero()) continue;
ps.append(pso.op(i));
ns.append(nso.op(i));
}
return F(ps,ns);
} else return e.map(self);
});
if(IBPmethod==0) {
auto air_res =
GiNaC_Parallel(air_vec.size(), 1, [&](int idx)->ex {
auto air = air_vec[idx];
air = subs_naive(air,IR2F);
air = subs_naive(air,rules_ints.first);
air = _F2ex(air);
auto cv_lst = mma_collect_lst(air, F(w1,w2));
air = 0;
for(auto cv : cv_lst) air += cv.op(1) * (cv.op(0));
return air;
}, "A2F");
for(auto fp : ibp_vec) delete fp;
system(("rm -rf "+wdir).c_str());
for(int i=0; i<air_vec.size(); i++) air_vec[i] = air_res[i];
return;
}
if(IBPmethod==1) {
for(auto ibp : ibp_vec_re) ibp->Export();
auto nproc = 2*CpuCores()/FIRE::Threads;
int cproc = 0;
if(nproc<2) nproc = 2;
#pragma omp parallel for num_threads(nproc) schedule(dynamic, 1)
for(int pi=0; pi<ibp_vec_re.size(); pi++) {
if(Verbose>1) {
#pragma omp critical
cout << "\r \r" << " \\--FIRE Reduction [" << (++cproc) << "/" << ibp_vec_re.size() << "] " << flush;
}
ibp_vec_re[pi]->Run();
}
if(Verbose>1) cout << "@" << now(false) << endl;
for(auto item : ibp_vec_re) item->Import();
system(("rm -rf "+wdir).c_str());
} else if(IBPmethod==2) {
for(auto ibp : ibp_vec_re) ibp->Reduce();
system(("rm -rf "+wdir).c_str());
}
vector<Base*> base_re;
for(auto f : ibp_vec_re) base_re.push_back(f);
auto mi_rules = FindRules(base_re, true, aip.UF);
auto rules_vec =
GiNaC_Parallel(ibp_vec_re.size(), 10, [&](int idx)->ex {
lst rules;
for(auto item : ibp_vec_re[idx]->Rules) {
auto rr = item.op(1);
rr = subs_naive(rr,mi_rules.first);
auto cv_lst = mma_collect_lst(rr, F(w1,w2));
rr = 0;
for(auto cv : cv_lst) rr += cv.op(1) * (cv.op(0));
rules.append(item.op(0)==rr);
}
return rules;
}, "F2F");
exmap F2F;
for(auto rule : rules_vec) {
for(auto lr : rule) F2F[lr.op(0)] = lr.op(1);
}
auto air_res =
GiNaC_Parallel(air_vec.size(), 1, [&](int idx)->ex {
auto air = air_vec[idx];
air = subs_naive(air,IR2F);
air = subs_naive(air,rules_ints.first);
air = subs_naive(air,mi_rules.first);
air = subs_naive(air,F2F);
air = _F2ex(air);
auto cv_lst = mma_collect_lst(air, F(w1,w2));
air = 0;
for(auto cv : cv_lst) air += cv.op(1) * (cv.op(0));
return air;
}, "A2F");
for(auto fp : ibp_vec) delete fp;
for(int i=0; i<air_vec.size(); i++) air_vec[i] = air_res[i];
}
}
......@@ -4,6 +4,7 @@
*/
#include "HEP.h"
#include "cln/cln.h"
namespace HepLib {
......@@ -694,5 +695,25 @@ namespace HepLib {
REGISTER_FUNCTION(Matrix, do_not_evalf_params().conjugate_func(mat_conj).set_return_type(return_types::commutative))
bool isZero(const ex & e) {
exset vs;
for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i) {
if(is_a<symbol>(*i) || is_a<Pair>(*i)) vs.insert(*i);
}
int n = 3;
for(int i=0; i<5; i++) {
exmap nsubs;
for(auto item : vs) {
nsubs[item] = ex(1)/numeric(cln::nextprobprime(n));
n++;
}
ex ret = e.subs(nsubs);
if(!normal(e).is_zero()) return false;
}
auto ret = normal_fermat(e);
return ret.is_zero();
}
}
......@@ -341,6 +341,22 @@ namespace HepLib {
ex ApartIRC(const ex & expr_in);
void ApartIBP(int IBPmethod, exvector &io_vec, const lst & loops_exts=lst{}, const lst & cut_props=lst{}, std::function<lst(const Base &, const ex &)> uf=IBP::LoopUF);
struct AIOption {
lst AInternal; // Internal for Apart
lst AExternal; // External for Apart
lst IInternal; // Internal for IBP
lst I_Internal; // _Internal for IBP
lst IExternal; // External for IBP
lst Cuts; // Cut Propagators. optional
lst CPairs; // Pairs in Cut, to be cleared. optional
lst IPairs; // Pair for IBP. optional
std::function<lst(const Base &, const ex &)> UF = IBP::LoopUF;
};
void ApartIBP(int IBPmethod, exvector &io_vec, AIOption aip);
bool isZero(const ex & e);
#ifndef DOXYGEN_SKIP
class ApartIR1_SERIAL { public: static unsigned serial; };
......
......@@ -74,7 +74,9 @@ namespace HepLib::IBP {
exvector IBPvec;
lst ns0;
for(int i=0; i<pdim; i++) ns0.append(0);
for(auto loop : Internal) {
auto Dloop = Internal;
if(_Internal.nops()>0) Dloop = _Internal;
for(auto loop : Dloop) {
lst dp_lst;
for(int i=0; i<pdim; i++) {
auto s = ex_to<Symbol>(loop);
......
......@@ -25,6 +25,7 @@ namespace HepLib::IBP {
class Base {
public:
lst Internal;
lst _Internal;
lst External;
lst Replacements;
lst Propagators;
......
......@@ -41,7 +41,7 @@ namespace HepLib {
int Verbose = 0;
int GiNaC_Parallel_Process = -1;
const Symbol D("D");
bool fermat_use_array = true;
bool fermat_using_array = true;
int NNDigits = 100;
HepFormat::_init::_init() {
......
......@@ -15,7 +15,7 @@ namespace HepLib::SD {
auto xs = get_x_from(expr);
auto nx = xs.size();
if(nx<1) return expr;
if(!is_polynomial(expr, exvec2lst(xs))) return expr;
if(!is_polynomial(expr, vec2lst(xs))) return expr;
auto cvs = mma_collect_lst(expr,x(w));
vector<exvector> degs_vec;
for(auto cv : cvs) {
......@@ -1375,7 +1375,7 @@ namespace HepLib::SD {
auto pref = mma_series(ct2, ep, epN-di);
if(pref.has(vs)) pref = mma_series(pref, vs, sN);
if(is_zero(intg.subs(x(w)==ex(1)/7)) && is_zero(intg.subs(x(w)==ex(1)/13))) {
intg = fermat_normal(intg);
intg = normal_fermat(intg);
}
para_res_lst.append(lst{pref * pow(ep, di), intg});
}
......@@ -1404,7 +1404,7 @@ namespace HepLib::SD {
auto pref = mma_series(ct2, ep, epN-di);
if(pref.has(vs)) pref = mma_series(pref, vs, sN);
if(is_zero(intg.subs(x(w)==ex(1)/7)) || is_zero(intg.subs(x(w)==ex(1)/13))) {
intg = fermat_normal(intg);
intg = normal_fermat(intg);
}
para_res_lst.append(lst{eps_ci * pref * pow(eps, sdi) * pow(ep, di), intg});
}
......
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