23 if(!is_polynomial(expr,
vec2lst(xs)))
return expr;
25 vector<exvector> degs_vec;
27 if(is_zero(cv.op(0)))
continue;
29 for(
auto xi : xs) degs.push_back(cv.op(1).degree(xi));
30 degs_vec.push_back(degs);
34 auto nn = degs_vec.size();
36 for(
int i=0; i<nn; i++) flag[i] =
true;
39 for(
int ri=0; ri<nn; ri++)
for(
int ci=0; ci<nx; ci++) mat(ri,ci) = degs_vec[ri][ci];
41 for(
auto vi : vvi)
for(
auto ii : vi) flag[ii] =
false;
43 for(
int i=0; i<nn; i++) flag[i] =
false;
46 for(
int i=0; i<nn; i++) {
48 exvector & ti = degs_vec[i];
49 for(
int j=i+1; j<nn; j++) {
51 exvector & tj = degs_vec[j];
53 for(
int c=0; c<nx; c++) {
55 if(comp>0) { comp=0;
break; }
57 }
else if(ti[c]>tj[c]) {
58 if(comp<0) { comp=0;
break; }
72 for(
int i=0; i<nn; i++) {
75 for(
int c=0; c<nx; c++) res *= pow(xs[c],degs_vec[i][c]);
91 for(
auto kv : map_vec[0]) {
93 if(is_zero(xi-x(-1)))
continue;
95 if(quick) nxs[xi] = xi.subs(x(
w)==ex(
w+1)/ex(
w+13));
99 for(
auto xi : xs) chk_total *= 1/(nxs[xi]+1);
102 for(
auto imap : map_vec) {
104 for(
auto kv : imap) {
106 auto nxi = kv.first.subs(x(
w)==(13+
w));
107 if(is_zero(xi-x(-1))) intg *= kv.second;
108 else intg *= pow(kv.second, nxs[kv.first]);
114 if(is_zero(tmp-intg))
break;
117 intg = intg.subs(pow(y(
w1),
w2)==1/(
w2+1)).subs(y(
w)==1/ex(2));
120 int_total = int_total.normal();
121 return normal(chk_total-int_total).is_zero();
133 if(in_xpols.has(y(
w)))
throw Error(
"SecDecBase::x2y: y(w) found @ " +
ex2str(in_xpols));
134 lst xpols_lst = in_xpols;
137 for(
auto item : xpols_lst) {
138 if(item.match(x(
w)))
continue;
139 else if(item.match(pow(x(
w1),
w2)))
continue;
140 else if(item.match(pow(
w1,
w2))) item = item.op(0);
144 if(!is_a<mul>(xpols)) xpols = lst{ xpols };
146 xpols_lst.remove_all();
147 for(
auto item : xpols) {
148 if(item.match(x(
w)))
continue;
149 else if(item.match(pow(x(
w1),
w2)))
continue;
150 else if(item.match(pow(
w1,
w2))) xpols_lst.append(item.op(0));
151 else xpols_lst.append(item);
157 for(
auto item : xpols_lst) xpols2 *= item;
158 if(is_a<lst>(xpols)) xpols = xpols.op(0);
159 if(is_zero(xpols2-xpols)) {
164 throw Error(
"SecDecBase::x2y: unexpected region reached.");
165 return vector<exmap>();
175 for(
auto &vi : vmap) {
176 auto ft = f.subs(vi);
179 int ysn = ys_tmp.size();
180 vector<int> y_free_indexs;
181 for(
int j=0; j<ys_tmp.size()+10; j++) {
182 if(!ft.has(y(j))) y_free_indexs.push_back(j);
184 if(xs_tmp.size()>y_free_indexs.size())
throw Error(
"IsBad: xsize>ysize");
185 for(
int i=0; i<xs_tmp.size(); i++) {
186 vi[xs_tmp[i]] = y(y_free_indexs[i]);
188 if(xs_tmp.size()>0) ft = ft.subs(vi);
191 ft = collect_common_factors(
expand_ex(ft,y(
w)));
192 if(is_exactly_a<mul>(ft)) {
194 for (
auto item : ft) {
195 if( !(item.match(y(
w)) || item.match(pow(y(
w),
w1))) ) {
200 }
else if ( ft.match(y(
w)) || ft.match(pow(y(
w),
w1)) ) {
205 for(
int pi=1; pi<std::pow(2, ys_tmp.size()); pi++) {
208 for(
int i=0; i<ys_tmp.size(); i++) {
209 yrepl.append(ys_tmp[i] == ((ci%2)==1 ? 1 : 0));
212 if(normal(ft.subs(yrepl)).is_zero())
return true;
224 if(po_ex.nops()>2)
throw Error(
"AutoEnd: Deltas found @ " +
ex2str(po_ex));
225 lst
const exlist = ex_to<lst>(po_ex.op(1));
226 if(!(exlist.op(0)-1).is_zero())
throw Error(
"AutoEnd: (!(exlist.op(0)-1).is_zero())");
230 for(
int nn=0; nn<=nx; nn++) {
231 for(
int pi=0; pi<std::pow(2, nx); pi++) {
232 int cpi = pi, cn1 = 0;
233 for(
int i=0; i<nx; i++) {
234 if((cpi % 2) == 1) cn1++;
237 if(cn1 != nn)
continue;
239 lst polists = lst{ po_ex.op(0) };
242 for(
int i=0; i<nx; i++) {
246 for(
auto item : polists) {
248 auto tmp = ex_to<lst>(
subs(
subs(item, lst{xs[i]==x0*xy}), lst{xy==xs[i]}));
249 tmp.let_op(0) = tmp.op(0) * x0;
250 polists2.append(tmp);
252 tmp = ex_to<lst>(
subs(
subs(item, lst{xs[i]==(x0-1)*xy+1}), lst{xy==xs[i]}));
253 tmp.let_op(0) = tmp.op(0) * (1-x0);
254 polists2.append(tmp);
262 for(
auto polist : polists) {
264 for(
int i=0; i<polist.nops(); i++) {
265 auto tmp = polist.op(i);
266 auto ntmp =
NN(exlist.op(i));
267 if(!tmp.subs(lst{x(w)==0, y(w)==0}).normal().is_zero())
continue;
268 if( (!tmp.has(x(
w)) && !tmp.has(y(
w))) || (is_a<numeric>(ntmp) && ntmp>0) )
continue;
272 vector<exmap> vmap =
SecDec->x2y(sdList);
274 for(
int ni=0; ni<polist.nops(); ni++) {
275 auto po = polist.op(ni);
276 auto expo = exlist.op(ni);
277 if(expo.info(info_flags::posint))
continue;
278 if(
IsBad(po, vmap)) {
288 for(
auto item : polists) res.push_back(lst{ex_to<lst>(item), exlist});
289 return exvector(std::move(res));
293 cerr <<
ErrColor <<
"polynormial list: " << po_ex.op(0) <<
RESET << endl;
294 throw Error(
"AutoEnd Failed @ ALL possible bisections!");
303 exvector SecDec::DS(
const ex po_ex) {
306 lst
const polist = ex_to<lst>(po_ex.op(0));
307 lst
const exlist = ex_to<lst>(po_ex.op(1));
309 for(
int i=0; i<polist.nops(); i++) {
310 auto tmp = polist.op(i);
311 auto ntmp = exlist.op(i);
312 if(!tmp.subs(lst{x(w)==0, y(w)==0}).normal().is_zero())
continue;
314 if( (!tmp.has(x(
w)) && !tmp.has(y(
w))) || (is_a<numeric>(ntmp) && ntmp>0) )
continue;
318 vector<exmap> vmap =
SecDec->
x2y(sdList);
320 for(
auto vi : vmap) cout << vi << endl;
321 throw Error(
"VerifySD Failed!");
325 for(
auto vi : vmap) {
326 auto ypolist = polist.subs(vi);
329 int ysn = ys_tmp.size();
330 vector<int> y_free_indexs;
331 for(
int j=0; j<xs_tmp.size()+ys_tmp.size()+10; j++) {
332 if(!ypolist.has(y(j))) y_free_indexs.push_back(j);
334 if(xs_tmp.size()>y_free_indexs.size())
throw Error(
"DS: xsize>ysize");
335 for(
int i=0; i<xs_tmp.size(); i++) {
336 vi[xs_tmp[i]] = y(y_free_indexs[i]);
338 if(xs_tmp.size()>0) ypolist = ypolist.subs(vi);
339 if(ypolist.has(x(
w)))
throw Error(
"DS: x(w) found @ " +
ex2str(ypolist));
342 auto ft = collect_common_factors(
expand_ex(ypolist.op(1),y(
w)));
344 ex ct = 1, fsgin = 1;
347 for (
auto item : ft) {
348 if( !(item.match(y(
w)) || item.match(pow(y(
w),
w1))) ) {
353 }
else if( ft.match(y(
w)) || ft.match(pow(y(
w),
w1)) ) {
365 eps_map[epi.op(0)] = epn;
369 bool need_contour_deformation = ft.has(PL(
w));
370 if(ft.has(y(
w)) && !need_contour_deformation) {
373 need_contour_deformation =
false;
374 auto first = tmp.op(0).subs(y(
w)==1);
375 for(
auto item : tmp) {
376 auto chk = item.subs(y(
w)==1)*first;
377 auto nchk =
NN(chk,30);
378 if(!is_a<numeric>(nchk)) {
379 throw Error(
"DS: Not a number: " +
ex2str(chk));
382 need_contour_deformation =
true;
392 if(!is_a<numeric>(tmp))
throw Error(
"DS: NOT a numeric with " +
ex2str(tmp));
394 ct = exp(-I * Pi * exlist.op(1));
403 if(is_a<add>(det))
throw Error(
"DS: det is add " +
ex2str(det));
406 for(
int i=0; i<ys.size(); i++) {
407 auto ndeg = det.degree(ys[i]);
408 if(ndeg!=det.ldegree(ys[i]))
throw Error(
"DS: det is not a monomial with det = "+
ex2str(det));
409 ymol[ys[i]] = ymol[ys[i]] + ndeg;
410 det1 *= pow(ys[i], ndeg);
412 if(!(is_a<numeric>(det/det1) && ex_to<numeric>(det/det1).is_integer())) {
413 throw Error(
"DS: det=" +
ex2str(det) +
", det1=" +
ex2str(det1));
420 pol_exp_lst.append(lst{
subs(FTX(ft,ftxlst)*CT(ct*det/det1), y(
w)==x(
w)), 1 });
422 for(
int i=0; i<ypolist.nops(); i++) {
423 auto tmp = ypolist.op(i);
424 auto nexp = exlist.op(i).normal();
425 bool nchk = (!is_a<numeric>(nexp) || (is_a<numeric>(nexp) && nexp<0));
426 if(nexp.has(x(
w)) || nexp.has(y(
w)))
throw Error(
"DS: x or y found in exp: "+
ex2str(nexp));
431 if(!tmp.is_polynomial(ys))
throw Error(
"DS: NOT a polynomial with "+
ex2str(tmp));
434 if(i==1 && fsgin<0) tmp = ex(0)-tmp;
438 tmp = tmp.subs(pow(pow(y(
w1),
w2),
w3)==pow(y(
w1),
w2*
w3));
439 tmp = tmp.subs(pow(sqrt(y(
w1)),
w2)==pow(y(
w1),
w2/2));
440 tmp = tmp.subs(sqrt(pow(y(
w1),
w2))==pow(y(
w1),
w2/2));
441 if((tmp1-tmp).is_zero())
break;
450 if(!is_a<mul>(tmps)) tmps = lst{tmps};
451 for (
auto item : tmps) {
452 if(item.match(y(
w))) {
454 ymol[yi] = ymol[yi] + nexp;
455 }
else if(item.match(pow(y(
w1),
w2))) {
456 auto yi = item.op(0);
457 ymol[yi] = ymol[yi] + item.op(1) * nexp;
458 }
else if(!item.has(y(
w)) && !item.has(x(
w))) {
459 if(is_a<numeric>(nexp) && ex_to<numeric>(nexp).is_integer()) {
461 }
else if(!item.has(PL(
w)) && !item.has(WRA(
w))) {
462 auto tr =
NN(item.subs(
nReplacement).subs(lst{ CV(w1,w2)==w2 }).subs(eps_map),30);
463 if(!is_a<numeric>(tr)) {
464 throw Error(
"DS: not numeric - item: " +
ex2str(tr) +
" ; " +
ex2str(item));
466 auto nitem = ex_to<numeric>(tr);
467 if( nitem.is_real() && nitem<0 ) {
477 if(nchk && item.subs(lst{y(w)==0,iEpsilon==0}).normal().is_zero()) {
478 throw Error(
"DS: zero item: " +
ex2str(item) +
" and exlist.op(i) = " +
ex2str(nexp));
484 ex pnp = rem-(i==1 && (ft!=1 || hasWRA) ?
iEpsilon : ex(0));
485 pnp = pnp.subs(y(
w)==x(
w));
487 pol_exp_lst.append(lst{pnp, nexp});
488 pol_exp_lst.let_op(0) = pol_exp_lst.op(0).subs(CT(
w)==CT(
w*pow(ct,nexp))).
subs(CT(0)==0);
492 for(
auto & kv : ymol) {
493 auto k = kv.first.subs(y(
w)==x(
w));
495 if(is_a<numeric>(v)) {
496 auto nv = ex_to<numeric>(v);
498 throw Error(
"DS: " +
ex2str(kv.first) +
"^(" +
ex2str(nv) +
") found, " +
ex2str(vi));
501 x_n_lst.append(lst{k, v});
504 sd.push_back(lst{x_n_lst, pol_exp_lst});
507 return exvector(std::move(sd));
512 lst SecDec::Normalize(
const ex &input) {
515 lst in_plst = ex_to<lst>(input.op(0));
516 lst in_nlst = ex_to<lst>(input.op(1));
519 int pnN = in_plst.nops();
520 for(
int i=0; i<pnN; i++) {
521 auto pi = in_plst.op(i);
522 auto ni = in_nlst.op(i);
524 if(!pi.is_polynomial(vars)) {
525 if(ni.info(info_flags::integer)) {
527 if(nd.op(0).is_polynomial(vars) && nd.op(1).is_polynomial(vars)) {
528 in_plst.append(nd.op(1));
529 in_nlst.append(ex(0)-ni);
530 in_plst.let_op(i) = nd.op(0);
535 if(nd.op(0).is_polynomial(vars) && nd.op(1).is_polynomial(vars)) {
537 in_plst.append(nd.op(1));
538 in_nlst.append(ex(0)-ni);
539 in_plst.let_op(i) = nd.op(0);
542 in_plst.append(nd.op(0));
544 in_plst.let_op(i) = nd.op(1);
545 in_nlst.let_op(i) = ex(0)-ni;
553 if(!is_a<mul>(pi)) pis = lst{ pi };
556 for(
auto item : pis) {
557 if(item.is_polynomial(vars)) rem *= item;
558 else if(item.match(pow(
w0,
w1))) {
561 if(pp.is_polynomial(vars) &&
xPositive(pp)) {
563 in_nlst.append(nn*ni);
566 if(nd.op(0).is_polynomial(vars) && nd.op(1).is_polynomial(vars)) {
568 in_plst.append(nd.op(0));
569 in_plst.append(nd.op(1));
570 in_nlst.append(nn*ni);
571 in_nlst.append(ex(0)-nn*ni);
572 }
else { cout << item << endl;
goto nok; }
573 }
else { cout << item << endl;
goto nok; }
577 in_plst.let_op(i) = rem;
581 nok:
throw Error(
"SecDec::Normalize, NOT polynomial found: "+
ex2str(pi));
586 for(
int i=0; i<in_plst.nops(); i++) {
587 if(i!=1 && (in_nlst.op(i).is_zero() || in_plst.op(i)==ex(1)))
continue;
588 if(i!=1 && !in_plst.op(i).has(x(
w)) && !in_plst.op(i).has(y(
w)) && !in_plst.op(i).has(z(
w))) {
589 if(in_nlst.op(i).info(info_flags::integer) ||
xPositive(in_plst.op(i))) {
590 const_term *= pow(in_plst.op(i), in_nlst.op(i));
591 }
else const_term *= exp(log(in_plst.op(i)) * in_nlst.op(i));
594 auto ptmp = in_plst.op(i);
595 auto ntmp = in_nlst.op(i);
596 if(!is_a<mul>(ptmp)) ptmp = lst{ ptmp };
601 eps_map[epi.op(0)] = epn;
605 for(
int j=0; j<ptmp.nops(); j++) {
606 auto tmp = ptmp.op(j);
607 if(!tmp.has(x(
w)) && !tmp.has(y(
w)) && !tmp.has(z(
w))) {
608 if(ntmp.info(info_flags::integer)) {
609 const_term *= pow(tmp,ntmp);
610 }
else if((tmp-
vs).is_zero() || tmp.match(pow(
vs,
w))) {
611 const_term *= pow(tmp,ntmp);
612 }
else if(!tmp.has(PL(
w)) && !tmp.has(
vs) && !tmp.has(WRA(
w))) {
613 auto tr =
NN(tmp.subs(
nReplacement).subs(lst{ CV(w1,w2)==w2 }).subs(eps_map),30);
614 if(!is_a<numeric>(tr) || !tr.info(info_flags::real)) {
615 cerr <<
"tmp: " << tmp << endl;
616 cerr <<
"tr: " << tr << endl;
618 throw Error(
"Normalize: tr is NOT numeric with nReplacement.");
622 const_term *= pow(tmp,ntmp);
624 const_term *= pow(-tmp,ntmp);
630 }
else if(tmp.match(pow(
w1,
w2)) && ntmp.info(info_flags::integer) && tmp.op(1).info(info_flags::integer)) {
631 plst.append(tmp.op(0));
632 nlst.append(ntmp * tmp.op(1));
634 plst.append(tmp.op(0));
635 nlst.append(ntmp * tmp.op(1));
636 }
else if(
xSign(tmp)!=0 ||
xSign(tmp.subs(lst{y(w)==x(w),z(w)==x(w)}))!=0) {
637 if(tmp.subs(lst{x(w)==1, y(w)==1, z(w)==1})>0) {
641 plst.append(ex(0)-tmp);
653 }
else if(tmul != 1) {
660 plst.prepend(const_term);
663 lst plst_comb, nlst_comb;
665 for(
int i=0; i<nlst.nops(); i++) {
666 if(i!=1) np[plst[i]] += nlst[i];
668 map<ex,int,ex_is_less> inp;
669 for(
int i=0; i<nlst.nops(); i++) {
671 plst_comb.append(plst[1]);
672 nlst_comb.append(nlst[1]);
674 if(inp[plst[i]]==1)
continue;
675 plst_comb.append(plst[i]);
676 nlst_comb.append(np[plst[i]]);
681 if(nlst_comb.op(0)!=1) {
682 if(nlst_comb.op(0).info(info_flags::integer) ||
xPositive(plst_comb.op(0))) {
683 plst_comb.let_op(0) = pow(plst_comb.op(0),nlst_comb.op(0));
684 }
else plst_comb.let_op(0) = exp(log(plst_comb.op(0))*nlst_comb.op(0));
685 nlst_comb.let_op(0) = 1;
688 return (input.nops()>2) ? lst{plst_comb, nlst_comb, input.op(2)} : lst{plst_comb, nlst_comb};
701 for(
int i=0; i<xs.size(); i++) x2y.append(xs[i]==y(i));
702 fe = ex_to<lst>(
subs(fe, x2y).
subs(y(
w)==x(
w)));
708 for(
int i=0; i<xs.size(); i++) x2y.append(xs[i]==y(i));
709 vint = vint.subs(x2y).subs(y(
w)==x(
w));
716 for(
int ri=0; ri<2; ri++) {
720 for(
auto fe :
FunExp) funexp.push_back(Normalize(fe));
722 funexp =
GiNaC_Parallel(
FunExp.size(), [&](
int idx)->ex { return Normalize(FunExp[idx]); },
"Normalize");
728 for(
auto fe : funexp) {
730 if(fe.nops()>2) key = iWF(fe.op(2));
731 if(fe.op(1).op(0)!=1) {
733 throw Error(
"Normalizes: fe.op(0).op(0) is NOT 1.");
735 for(
int i=1; i<fe.op(0).nops(); i++) key *= pow(fe.op(0).op(i), fe.op(1).op(i));
736 fn[key] += fe.op(0).op(0);
739 map<ex,int,ex_is_less> ifn;
740 for(
auto fe : funexp) {
742 if(fe.nops()>2) key = iWF(fe.op(2));
743 for(
int i=1; i<fe.op(0).nops(); i++) key *= pow(fe.op(0).op(i), fe.op(1).op(i));
744 if(ifn[key]==1)
continue;
746 funs.append(fn[key]);
748 for(
int i=1; i<fe.op(0).nops(); i++) {
749 funs.append(fe.op(0).op(i));
750 exps.append(fe.op(1).op(i));
752 if(fe.nops()>2)
FunExp.push_back(lst{funs, exps, fe.op(2)});
753 else FunExp.push_back(lst{funs, exps});
763 funexp.push_back(fe);
768 map<ex, ex, ex_is_less> fe_cc;
769 for(
auto fe : funexp) {
773 fun.append(fe.op(0).op(1));
774 exp.append(fe.op(1).op(1));
776 for(
int i=0; i<fe.op(1).nops(); i++) {
778 auto expi = fe.op(1).op(i);
779 if(expi.info(info_flags::nonnegint)) {
780 rem *= pow(fe.op(0).op(i), expi);
782 fun.append(fe.op(0).op(i));
783 exp.append(fe.op(1).op(i));
792 for(
auto kv : fe_cc) {
793 lst fe = ex_to<lst>(kv.first);
804 funexp.push_back(fe);
809 for(
auto fe : funexp) {
810 lst fun = lst{fe.op(0).op(0), fe.op(0).op(1)};
811 lst exp = lst{fe.op(1).op(0), fe.op(1).op(1)};
813 for(
int i=2; i<fe.op(1).nops(); i++) {
814 auto expi = fe.op(1).op(i);
815 if(expi.info(info_flags::nonnegint)) {
816 rem *= pow(fe.op(0).op(i), expi);
818 fun.append(fe.op(0).op(i));
819 exp.append(fe.op(1).op(i));
827 for(
auto item : rem) rem_lst.append(item);
832 for(
auto item : rem_lst) {
847 if(
Verbose > 2) cout <<
PRE <<
"\\--Scaleless: " <<
FunExp.size() <<
" :> " << flush;
852 auto funexp = FunExp[idx];
853 if(funexp.nops()<3) return funexp;
855 auto fun = funexp.op(0);
856 auto exp = funexp.op(1);
857 auto deltas = funexp.op(2);
859 for(int di=0; di<deltas.nops(); di++) {
860 auto delta = ex_to<lst>(deltas.op(di));
862 if(delta.nops()<2) continue;
866 for(int j=0; j<delta.nops(); j++) {
867 sRepl.append(delta[j]==delta[j]*s);
872 for(int j=0; j<fun.nops(); j++) {
873 auto tmp = expand(fun.op(j).subs(sRepl));
874 if(tmp.degree(s)!=tmp.ldegree(s)) {
878 n_s += tmp.degree(s) * exp.op(j);
880 if(!is_s || !normal(n_s+delta.nops()).is_zero()) continue;
883 for(size_t i=1; i<ex_to<numeric>(GiNaC::pow(2,delta.nops())).to_long()-1; i++) {
887 for(int j=0; (j<delta.nops() && ci>0); j++) {
889 sRepl.append(delta[j]==delta[j]*s);
896 for(int j=0; j<fun.nops(); j++) {
897 if(exp.op(j).info(info_flags::nonnegint)) continue;
898 auto tmp = collect_ex(fun.op(j).subs(sRepl),s);
899 if(tmp.degree(s)!=tmp.ldegree(s)) {
903 n_s += tmp.degree(s) * exp.op(j);
906 if(!normal(n_s).is_zero()) {
913 if(!is0)
return lst{fun, exp, deltas};
919 FunExp.shrink_to_fit();
920 for(
auto item : sl_res) {
921 lst fed = ex_to<lst>(item);
922 if(fed.nops()<1)
continue;
923 FunExp.push_back(fed);
926 if(
Verbose > 2) cout << FunExp.size() << endl;
927 if(FunExp.size()<1)
IsZero =
true;
930 void SecDec::RemoveDeltas() {
935 exvector funexp = FunExp;
937 FunExp.shrink_to_fit();
938 for(
auto fe : funexp) {
940 FunExp.push_back(fe);
943 auto xs = ex_to<lst>(fe.op(2).op(0));
945 for(
int i=1; i<fe.op(2).nops(); i++) {
946 re_deltas.append(fe.op(2).op(i));
948 if(re_deltas.nops()>0) exit =
false;
951 FunExp.push_back(lst{fe.op(0), fe.op(1), re_deltas});
953 for(
int i=0; i<xs.nops(); i++) {
956 for(
auto xi : xs) jInv += xi;
957 jInv = jInv.subs(xj==1);
960 for(
int j=0; j<xs.nops(); j++) {
962 if(xxj != xj) repl[xxj] = xj*xxj;
966 lst exps = ex_to<lst>(fe.op(1));
968 for(
int j=0; j<fe.op(0).nops(); j++) {
969 auto fun = fe.op(0).op(j);
971 if(!fun.is_polynomial(xj)) {
972 cerr <<
"xj: " << xj << endl;
973 cerr <<
"fun: " << fun << endl;
974 cerr << funexp << endl;
975 throw Error(
"RemoveDeltas: fun is NOT polynormial of xj.");
977 auto expn = fun.degree(xj);
978 fun = pow(xj, -expn) * fun;
980 fun = fun.subs(xj==jInv);
982 expns += expn * exps.op(j);
986 exps.append(ex(0)-xs.nops()-expns);
987 if(re_deltas.nops()>0) FunExp.push_back(lst{funs, exps, re_deltas});
988 else FunExp.push_back(lst{funs, exps});
995 if(use_Normalizes) Normalizes();
998 void SecDec::XEnd() {
999 if(
Verbose > 2) cout <<
PRE <<
"\\--BiSection: " << FunExp.size() <<
" :> " << flush;
1004 auto fe = FunExp[idx];
1006 if(xSign(fe.op(0).op(1))==0 && !fe.has(WRA(w))) {
1007 auto fes = AutoEnd(fe);
1008 for(auto fei : fes) para_res_lst.append(fei);
1010 para_res_lst.append(fe);
1012 return para_res_lst;
1017 FunExp.shrink_to_fit();
1018 for(
auto &item : funexps) {
1019 for(
auto &it : ex_to<lst>(item)) FunExp.push_back(ex_to<lst>(it));
1021 if(
Verbose > 2) cout << FunExp.size() << endl;
1024 void SecDec::BiSection(ex xi, ex x0) {
1027 for(
auto fe : fes) {
1029 auto tmp = ex_to<lst>(
subs(
subs(fe, lst{xi==x0*xy}), lst{xy==xi}));
1030 tmp.let_op(0).let_op(0) = tmp.op(0).op(0) * x0;
1031 FunExp.push_back(tmp);
1033 tmp = ex_to<lst>(
subs(
subs(fe, lst{xi==(x0-1)*xy+1}), lst{xy==xi}));
1034 tmp.let_op(0).let_op(0) = tmp.op(0).op(0) * (1-x0);
1035 FunExp.push_back(tmp);
1040 void SecDec::SDPrepares() {
1042 if(FunExp.size()<1) {
1050 if(CheckEnd) XEnd();
1053 Integrands.shrink_to_fit();
1056 FunExp.shrink_to_fit();
1057 for(
auto &fe : fes) {
1059 for(
auto item : fe.op(0)) {
1060 if(item.is_zero()) {
1065 if(to_add) FunExp.push_back(fe);
1067 if(FunExp.size()<1) {
1073 for(
auto epi : eps_lst) eps_map[epi.op(0)] = 0;
1080 auto fe = FunExp[idx];
1081 auto xns_pns = DS(fe);
1083 for(auto const &item : xns_pns) {
1087 auto ct = item.op(1).op(0).op(0);
1088 ct = ct.subs(lst{ CT(w)==w,FTX(w1,w2)==1 }).subs(pow(vs,vz)==1);
1089 int sNN = vsN - vsRank(ct.subs(pow(vs,vz)==1));
1093 for(int vn=0; vn<=sNN; vn++) zpols.append(vn);
1096 for(auto xn : item.op(0)) {
1097 auto xn_op1 = expand(xn.op(1).normal());
1098 ex c1 = xn_op1.coeff(vz);
1099 ex c0 = xn_op1.subs(vz==0);
1101 if(!is_a<numeric>(c1)) {
1102 cerr << ErrColor <<
"SDPrepares: c1 is not a number: " << c1 << RESET << endl;
1109 ex zp = (pxn-c0)/c1;
1110 ex zpn = zp.subs(eps_map).subs(lst{ep==0,epz==0});
1111 if(!is_a<numeric>(zpn)) {
1112 cerr << ErrColor <<
"SDPrepares: zpn is not a number: " << zpn << RESET << endl;
1125 for(auto zp : zpols) {
1126 auto item2 = item.subs(vz==ss+zp).subs(ss==vz);
1127 para_res_lst.append(item2);
1130 para_res_lst.append(item);
1133 return para_res_lst;
1136 ex min_expn = 1, min_expn2 = 10;
1137 exvector ibp_in_vec;
1138 for(
auto &item : sd_res) {
1139 for(
auto &it : ex_to<lst>(item)) {
1141 for(
auto xn : it.op(0)) {
1142 ex nxn = xn.op(1).subs(eps_map).subs(lst{
ep==0,
vz==0});
1143 if(nxn<-1) expn += nxn+1;
1144 if(min_expn2>nxn) min_expn2 = nxn+1;
1146 if(expn < min_expn) min_expn = expn;
1148 lst xns = ex_to<lst>(it.op(0));
1150 for(
int i=0; i<it.op(1).nops(); i++) {
1151 lst pn = ex_to<lst>(it.op(1).op(i));
1156 if(pn.op(0).is_equal(1))
continue;
1157 if(is_zero(pn.op(1))) {
1158 if(is_zero(pn.op(0)))
throw Error(
"SDPrepares: 0^0 found.");
1163 ibp_in_vec.push_back(lst{xns, pns});
1167 if(
Verbose > 1) cout <<
PRE <<
"\\--" <<
Color_HighLight <<
"Maximum x^-n: All(" << ex(0)-min_expn <<
"+1X), Max(" << (ex(0)-min_expn2) <<
"+1)" <<
RESET << endl;
1170 exvector ibp_res_vec;
1171 while(ibp_in_vec.size()>0) {
1174 spn <<
"XIBP-" << (pn-1);
1176 GiNaC_Parallel(ibp_in_vec.size(), [&eps_map,&ibp_in_vec,
this](
int idx)->ex {
1182 auto xns_pns = ibp_in_vec[idx];
1184 auto xns = xns_pns.op(0);
1185 auto pns = xns_pns.op(1);
1188 pns.op(0).find(FTX(w1,w2), fts);
1189 bool noFT = (fts.size()==1) && ( is_zero((*(fts.begin())).op(0)-1) );
1190 if(!noFT || disable_Contour) noFT = true;
1192 ex pole_requested = -1;
1193 if(noFT || PoleRequested > -1) pole_requested = PoleRequested;
1195 for(int n=0; n<xns.nops(); n++) {
1197 auto expn = xn.op(1).subs(eps_map).subs(lst{ep==0,vz==0,epz==0}).normal();
1198 if(!is_a<numeric>(expn)) throw Error(
"SDPrepares: expn NOT numeric: " + ex2str(expn));
1200 if(ex_to<numeric>(expn) < pole_requested) {
1202 pns.let_op(0).let_op(0) = pns.op(0).op(0) / (xn.op(1)+1);
1205 for(int i=0; i<xns.nops(); i++) {
1206 if(i!=n) xns2.append(xns.op(i));
1208 lst pns2 = ex_to<lst>(pns);
1209 for(int i=0; i<pns.nops(); i++) {
1210 pns2.let_op(i).let_op(0) = pns2.op(i).op(0).subs(xx==1);
1214 xns_pns_lst.append(lst{xns2, pns2});
1216 for(int i=0; i<pns.nops(); i++) {
1217 lst xns3 = ex_to<lst>(xns);
1218 xns3.let_op(n).let_op(1) = xn.op(1)+1;
1220 ex tmp = ex(0)-pns.op(i).op(1)*diff_ex(pns.op(i).op(0),xx);
1221 if(tmp.is_zero()) continue;
1223 auto xs = get_x_from(tmp);
1226 if(tmp.subs(xi==0).is_zero()) {
1232 if(tz) tmp = collect_common_factors(expand_ex(tmp,x(w)));
1233 if(tmp.is_zero()) continue;
1235 if(tz && is_a<mul>(tmp)) {
1237 for(auto ii : tmp) {
1238 if(ii.match(pow(x(w), w2))) {
1240 for(int ij=0; ij<xns3.nops(); ij++) {
1241 if(xns3.op(ij).op(0)==ii.op(0)) {
1242 xns3.let_op(ij).let_op(1) += ii.op(1);
1247 if(t) xns3.append(lst{ii.op(0), ii.op(1)});
1248 } else if(ii.match(x(w))) {
1250 for(int ij=0; ij<xns3.nops(); ij++) {
1251 if(xns3.op(ij).op(0)==ii) {
1252 xns3.let_op(ij).let_op(1) += 1;
1257 if(t) xns3.append(lst{ii, 1});
1265 lst pns3 = ex_to<lst>(pns);
1266 if(is_zero(pns.op(i).op(1)-1)) {
1267 pns3.let_op(i).let_op(0) = tmp;
1269 pns3.let_op(i).let_op(1) = pns.op(i).op(1)-1;
1270 int nn = pns.nops();
1271 if(!(pns.op(nn-1).op(1)-1).is_zero()) pns3.append(lst{ tmp, 1 });
1272 else pns3.let_op(nn-1).let_op(0) = pns.op(nn-1).op(0) * tmp;
1274 xns_pns_lst.append(lst{xns3, pns3});
1276 return lst{1, xns_pns_lst};
1279 return lst{0, xns_pns };
1281 }, spn.str().c_str());
1286 auto ii = ibp_res[idx];
1287 auto check = ii.op(0);
1289 auto items = ii.op(1);
1290 for(auto &it : ex_to<lst>(items)) ret1.append(it);
1293 sub_exp[pow(exp(w1),w2)]=exp(w1*w2);
1294 sub_exp[sqrt(exp(w1))]=exp(w1/2);
1295 sub_exp[exp(w1)*exp(w2)*w0]=exp(w1+w2)*w0;
1296 sub_exp[exp(w1)*exp(w2)]=exp(w1+w2);
1298 auto item = ii.op(1);
1299 ex expr = 1, expr_exp = 1;
1300 for(auto pn : item.op(1)) {
1301 if(pn.op(0).has(CT(w)) || pn.op(0).has(FTX(w1,w2))) {
1302 if(pn.op(1)!=1) throw Error(
"SDPrepares: exponent of CT is NOT 1.");
1304 } else if(xPositive(pn.op(0)) || pn.op(1).info(info_flags::integer)) {
1305 if(pn.op(0).has(exp(w))) expr_exp *= pow(pn.op(0), pn.op(1));
1306 else expr *= pow(pn.op(0), pn.op(1));
1308 expr_exp *= exp(log(pn.op(0)) * pn.op(1));
1312 auto expo = expr_exp.subs(sub_exp);
1313 if(is_zero(expo-expr_exp)) break;
1317 ret2.append(lst{ item.op(0), expr });
1319 return lst{ret1, ret2};
1322 ibp_in_vec.shrink_to_fit();
1323 for(
auto rets : rets_vec) {
1324 for(
auto item : rets.op(0)) ibp_in_vec.push_back(item);
1325 for(
auto item : rets.op(1)) ibp_res_vec.push_back(item);
1330 GiNaC_Parallel(ibp_res_vec.size(), [&ibp_res_vec,&eps_map](
int idx)->ex {
1333 auto xns_expr = ibp_res_vec[idx];
1335 auto xns = xns_expr.op(0);
1336 auto expr = xns_expr.op(1);
1337 lst exprs = { expr };
1339 for(
auto xn : xns) {
1340 auto expn = xn.op(1).subs(eps_map).subs(lst{ep==0,vz==0,epz==0}).normal();
1341 if(!is_a<numeric>(expn))
throw Error(
"SDPrepares: Not a number with expn = " +
ex2str(expn));
1344 for(
auto it : exprs) {
1345 ex rem = pow(xn.op(0), xn.op(1)) * it;
1346 if(ex_to<numeric>(expn)<=-1) {
1348 ex dit0 = dit.subs(xn.op(0)==0);
1350 if(!is_zero(dit0)) {
1351 rem -= pow(xn.op(0), xn.op(1)) * dit0 / ifact;
1352 exprs2.append(dit0/(xn.op(1)+1)/ifact);
1354 for(
int i=1; i+expn<0; i++) {
1356 if(is_zero(dit))
break;
1357 dit0 = dit.subs(xn.op(0)==0);
1359 if(!is_zero(dit0)) {
1360 rem -= pow(xn.op(0), xn.op(1)+i) * dit0 / ifact;
1361 exprs2.append(dit0/(xn.op(1)+i+1)/ifact);
1370 for(
auto const &it : exprs) {
1371 if(it.is_zero())
continue;
1374 for(
int i=0; i<xs.size(); i++) x2y.append(xs[i]==y(i));
1375 para_res_lst.append(it.subs(x2y).subs(y(
w)==x(
w)));
1378 return para_res_lst;
1382 bool zResides =
false;
1384 for(
auto &item : res) {
1385 for(
auto it : ex_to<lst>(item)) {
1386 if(!it.is_zero()) ints.push_back(it);
1387 if(!zResides && it.has(
vz)) zResides =
true;
1394 auto item = ints[idx];
1400 for(auto ct : cts) {
1404 if(is_a<mul>(ct.op(0))) {
1405 for(auto ii : ct.op(0)) cls.append(ii);
1407 cls.append(ct.op(0));
1409 for(auto cl : cls) {
1410 if(cl.has(vz)) cc *= cl;
1413 if(cc!=1) repl.append(ct==cc*CT(ll));
1416 it = series_ex(it,vz,-1);
1417 it = ex(0)-it.coeff(vz, -1);
1427 void SecDec::EpsExpands() {
1429 if(Integrands.size()<1) {
1436 if(
Verbose > 1) cout <<
PRE <<
"\\--Collecting: " << Integrands.size() <<
" :> " << flush;
1438 for(
auto &item : Integrands) {
1439 if(item.is_zero())
continue;
1441 item.find(CT(
w), cts);
1442 if(cts.size() != 1) {
1443 cerr <<
"cts: " << cts << endl;
1444 cerr <<
"item: " << item << endl;
1445 throw Error(
"EpsExpands: CT size is NOT 1: ");
1447 ex ct = (*(cts.begin())).subs(CT(
w)==
w);
1448 auto it = item.subs(CT(
w)==1);
1449 int_map[it] = int_map[it]+ct;
1453 Integrands.shrink_to_fit();
1454 for(
auto kv : int_map) {
1455 if(kv.second.is_zero())
continue;
1456 Integrands.push_back(CT(kv.second) * kv.first);
1458 if(
Verbose > 1) cout << Integrands.size() << endl;
1466 auto item = Integrands[idx];
1467 if(item.is_zero()) return lst{ lst{0, 0} };
1469 item.find(CT(
w), cts);
1470 if(cts.size() != 1) {
1471 cerr <<
"cts: " << cts << endl;
1472 cerr <<
"item: " << item << endl;
1473 throw Error(
"EpsExpands: CT size is NOT 1: ");
1475 ex ct = (*(cts.begin())).subs(CT(
w)==
w);
1476 auto it = item.subs(CT(
w)==1);
1480 ex ct0 = 1, ct1 = 1;
1482 if(ci.has(epz)) ct1 *= ci;
1496 for(
auto cv : cvs) {
1497 auto cc = cv.op(1) * ct;
1499 for(
auto epi : eps_lst) {
1500 auto epis = epi.op(0);
1501 if(!cc.has(epis) && !vv.has(epis))
continue;
1502 int iN =
ex2int(epi.op(1));
1505 int vN = vv.ldegree(epis);
1510 for(
auto item : eps_lst) pat_lst.append(item.op(0));
1512 for(
auto ic : ics) {
1514 ex pref = ic.op(1)*cc;
1515 ex n1 = intg.subs(x(
w)==ex(1)/7);
1517 n1 = intg.subs(x(
w)==ex(1)/13);
1518 if(is_zero(n1)) intg =
exnormal(intg);
1520 for(
auto epi : eps_lst) pref =
series_ex(pref, epi.op(0),
ex2int(epi.op(1)));
1521 para_res_lst.append(lst{pref, intg});
1525 return para_res_lst;
1529 if(
Verbose > 1) cout <<
PRE <<
"\\--Collecting: ";
1530 map<ex, ex, ex_is_less> int_pref;
1531 size_t ncollect = 0;
1533 for(
auto &item : res) {
1534 ncollect += item.nops();
1535 for(
auto &kv : ex_to<lst>(item)) {
1536 if(!kv.op(1).has(x(
w))) expr_nox += kv.op(0) * kv.op(1);
1537 else int_pref[kv.op(1)] += kv.op(0);
1541 if(
Verbose > 1) cout << ncollect <<
" :> " << flush;
1543 expResult.shrink_to_fit();
1544 for(
auto kv : int_pref) {
1545 if(kv.second.is_zero())
continue;
1546 expResult.push_back(lst{kv.second, kv.first});
1548 if(!is_zero(expr_nox)) {
1549 expr_nox = expr_nox.subs(FTX(
w1,
w2)==1);
1550 expResult.push_back(lst{expr_nox, 1});
1552 if(
Verbose > 1) cout << expResult.size() << endl;
class used to wrap error message
virtual vector< exmap > x2y(const ex &xpol)=0
static bool VerifySD(vector< exmap > vmap, bool quick=true)
Verify the Sector Decompostion is valid.
static ex XMonomials(const ex &expr)
XMonomials.
SecDec by geometric method.
static vector< vector< int > > RunQHull(const matrix &pts)
SecDec the main class to use Sector Decompostion method.
bool IsBad(ex f, vector< exmap > vmap)
Check the x-END.
static bool VerifySD(vector< exmap > vmap, bool quick=true)
exvector AutoEnd(ex po_ex)
Auto BiSection.
namespace for Numerical integration with Sector Decomposition method
exvector get_xy_from(ex pol)
ex xyz_pow_simplify(const ex expr_in)
exvector get_y_from(ex pol)
exvector get_x_from(ex pol)
int epsRank(ex expr_in, ex epi)
ex expand_ex(ex const &expr_in, std::function< bool(const ex &)> has_func)
the expand like Mathematica
const char * Color_HighLight
bool xPositive(ex const expr)
check the expr is xPositive, i.e., each x-monomial item is postive
ex NN(ex expr, int digits)
the nuerical evaluation
map< string, bool > GiNaC_Parallel_RM
ex collect_ex(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica
ex series_ex(ex const &expr_in, ex const &s0, int sn0)
the series like Mathematica, include s^n
string now(bool use_date)
date/time string
ex exnormal(const ex &expr, int opt)
normalize a expression
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}, ... }
ex diff_ex(ex const expr, ex const xp, unsigned nth, bool expand)
the differential like Mathematica
void let_op_append(ex &ex_in, const ex item)
append item into expression
lst vec2lst(const exvector &ev)
convert exvector to lst
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
ex exfactor(const ex &expr, int opt)
factorize a expression
int xSign(ex const expr)
the always sign for expr
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
int ex2int(ex num)
ex to integer
ex subs(const ex &e, init_list sl)