10 return n * UINT64_C(0x4f1bbcdd);
33 const symbol & get_symbol(
const string & s) {
34 static map<string, symbol> dict;
36 if (dict.find(key) == dict.end()) dict[key] = symbol(s,0);
40 const Symbol & get_Symbol(
const string & s) {
41 static map<string, Symbol> dict;
43 if (dict.find(key) == dict.end()) dict[key] = Symbol(s);
47 const iSymbol & get_iSymbol(
const string & s) {
48 static map<string, iSymbol> dict;
50 if (dict.find(key) == dict.end()) dict[key] = iSymbol(s);
59 GiNaC::registered_class_info &
Symbol::get_class_info_static() {
return reg_info; }
76 if(!is_a<Symbol>(other))
throw Error(
"Symbol::compare_same_type");;
78 int ret = get_name().compare(o.get_name());
80 else if(ret<0)
return -1;
89 inherited::archive(n);
97 inherited::read_archive(n);
106 static std::hash<std::string> hs;
107 unsigned seed = hs(get_name()+
typeid(*this).name());
109 setflag(status_flags::hash_calculated);
114 if(!is_a<Symbol>(other))
throw Error(
"Symbol::is_equal_same_type");
116 return get_name()==o->get_name();
120 if(!is_exactly_a<Symbol>(s))
return 0;
127 const ex point = r.rhs();
129 if (is_exactly_a<Symbol>(r.lhs()) && this->is_equal_same_type(ex_to<Symbol>(r.lhs()))) {
130 if (order>0 && !point.is_zero()) seq.emplace_back(expair(point, 0));
131 if (order>1) seq.emplace_back(expair(1, 1));
132 else seq.emplace_back(expair(Order(1), numeric(order)));
133 }
else seq.emplace_back(expair(*
this, 0));
134 return pseries(r, std::move(seq));
150 ex v2 = v1.subs(
vmap);
151 if(v2.is_equal(v1))
break;
164 GiNaC::registered_class_info &
iSymbol::get_class_info_static() {
return reg_info; }
180 if(!is_a<iSymbol>(other))
throw Error(
"iSymbol::compare_same_type");
182 int ret = get_name().compare(o.get_name());
184 else if(ret<0)
return -1;
189 if(!is_a<iSymbol>(other))
throw Error(
"iSymbol::is_equal_same_type");
191 return serial==o->serial;
195 if(!is_exactly_a<iSymbol>(s))
return 0;
202 const ex point = r.rhs();
204 if (is_exactly_a<iSymbol>(r.lhs()) && this->is_equal_same_type(ex_to<iSymbol>(r.lhs()))) {
205 if (order>0 && !point.is_zero()) seq.emplace_back(expair(point, 0));
206 if (order>1) seq.emplace_back(expair(1, 1));
207 else seq.emplace_back(expair(Order(1), numeric(order)));
208 }
else seq.emplace_back(expair(*
this, 0));
209 return pseries(r, std::move(seq));
215 hashvalue = get_symbol(get_name()).gethash();
216 setflag(status_flags::hash_calculated);
225 inherited::archive(n);
233 inherited::read_archive(n);
250 static int GiNaC_Parallel_Level = 0;
260 exvector
GiNaC_Parallel(
int ntotal, std::function<ex(
int)> f,
const string & key) {
261 if(ntotal<1)
return exvector();
265 int fork_retried = 0;
280 for(
int i=0; i<ntotal; i++) ovec.push_back(f(i));
281 return exvector(std::move(ovec));
284 int para_max_run = nproc<0 ? omp_get_num_procs() : nproc;
285 if(para_max_run<0) para_max_run = 0;
286 if(para_max_run>omp_get_num_procs()) para_max_run = omp_get_num_procs();
289 auto ppid = getpid();
291 cmd <<
"mkdir -p " << ppid;
292 if(!
dir_exists(to_string(ppid))) rc = system(cmd.str().c_str());
296 if(nbatch<=0) nbatch = ntotal/para_max_run/5;
297 else if(nbatch > ntotal/para_max_run) nbatch = ntotal/para_max_run;
298 if(nbatch<1) nbatch = 1;
299 int btotal = ntotal/nbatch + ((ntotal%nbatch)==0 ? 0 : 1);
301 bool nst = (GiNaC_Parallel_Level>0);
303 if(getpgid(0)!=ppid) {
306 if(npid < 0)
throw Error(
"GiNaC_Parallel: Error (1) @ fork()");
307 if(npid!=0)
goto wait_label;
310 if(setpgid(0,0))
throw Error(
"GiNaC_Parallel: setpgid(0,0) Failed.");
315 while(fork_retried<5) {
319 for(
int bi=0; bi<btotal; bi++) {
321 for(
int i=0; i<btotal; i++)
if(waitpid(-pgid,NULL,WNOHANG)>0) nactive--;
322 if(verb > 1 && !nst) {
323 cout <<
"\r \r" << pre;
324 cout <<
"\\--Evaluating ";
326 cout <<
Color_HighLight << nbatch <<
"x" <<
RESET <<
"[" << (bi+1) <<
"/" << btotal <<
"] @ " <<
now(
false) << exstr << flush;
331 if(key ==
"") garfn << ppid <<
"/" << bi <<
".gar";
332 else garfn << ppid <<
"/" << bi <<
"." << key <<
".gar";
338 if(getpid()!=pgid) exit(0);
340 if(ec>3*btotal)
throw Error(
"GiNaC_Parallel: Error (2) @ fork()");
341 exstr =
" [fork(" + to_string(ec) +
")]";
342 if(waitpid(-pgid,NULL,0)>0) nactive--;
345 if(pid==0 && setpgid(0, pgid)!=0) {
346 if(setpgid(0, pgid))
throw Error(
"GiNaC_Parallel: setpgid(0, pgid) Failed.");
350 if(nactive >= para_max_run)
if(waitpid(-pgid,NULL,0)>0) nactive--;
356 GiNaC_Parallel_Level++;
359 for(
int ri=0; ri<nbatch; ri++) {
360 int i = bi*nbatch + ri;
361 if(i<ntotal) res_lst.append(f(i));
365 if(key ==
"") garfn << ppid <<
"/" << bi <<
".gar";
366 else garfn << ppid <<
"/" << bi <<
"." << key <<
".gar";
368 }
catch(exception &p) {
369 cout <<
ErrColor <<
"Failed in GiNaC_Parallel!" <<
RESET << endl;
375 if(getpid()!=pgid) exit(0);
376 while (waitpid(-pgid,NULL,0) != -1) { }
377 if(!nst && verb > 1) cout << endl;
378 if(getpid()!=ppid) exit(0);
382 bool all_gar_exists =
true;
383 for(
int bi=0; bi<btotal; bi++) {
385 if(key ==
"") garfn << ppid <<
"/" << bi <<
".gar";
386 else garfn << ppid <<
"/" << bi <<
"." << key <<
".gar";
388 all_gar_exists =
false;
392 if(all_gar_exists)
break;
393 para_max_run = para_max_run/2;
394 if(para_max_run<2)
break;
399 if(npid!=0) waitpid(npid,NULL,0);
408 for(
int bi=0; bi<btotal; bi++) {
409 if(verb > 1 && !nst) {
411 cout <<
"\r \r" << pre;
412 cout <<
"\\--Reading *.gar [" << (bi+1) <<
"/" << btotal <<
"] @ " <<
now(
false) << flush;
414 cout <<
"\r \r" << pre;
415 cout <<
"\\--Reading *." <<
Color_HighLight << key <<
RESET <<
".gar [" << (bi+1) <<
"/" << btotal <<
"] @ " <<
now(
false) << flush;
420 if(key ==
"") garfn << ppid <<
"/" << bi <<
".gar";
421 else garfn << ppid <<
"/" << bi <<
"." << key <<
".gar";
425 res_lst = ex_to<lst>(
garRead(garfn.str()));
426 remove(garfn.str().c_str());
429 }
catch(exception &p) { }
430 if(verb > 1 && !nst) cout <<
" - ReTry" << endl;
432 res_lst.remove_all();
433 for(
int ri=0; ri<nbatch; ri++) {
434 int i = bi*nbatch + ri;
435 if(i<ntotal) res_lst.append(f(i));
438 }
catch(exception &p) {
439 cout <<
ErrColor <<
"Failed in GiNaC_Parallel!" <<
RESET << endl;
441 throw Error(
"GiNaC_Parallel_ReTRY: "+
string(p.what()));
446 for(
auto res : res_lst) ovec_tmp.push_back(res);
448 for(
auto res : res_lst) ovec.push_back(res);
451 for(
auto res : res_lst) ovec.push_back(res);
459 cout <<
"\r \r" << pre;
460 cout <<
"\\--ReWR" << flush;
462 cout <<
"\r \r" << pre;
469 if(key ==
"") garfn << ppid <<
"/ReWR.gar";
470 else garfn << ppid <<
"/ReWR." << key <<
".gar";
477 if(!nst && verb>1) cout << endl;
483 cmd <<
"rm -fr " << ppid;
484 rc = system(cmd.str().c_str());
486 if(ovec.size() != ntotal) {
487 throw Error(
"GiNaC_Parallel: The output size is wrong!");
490 return exvector(std::move(ovec));
499 std::vector<std::string>
split(
const std::string& s,
char delimiter) {
500 std::vector<std::string> tokens;
502 std::istringstream tokenStream(s);
503 while (std::getline(tokenStream, token, delimiter)) {
504 tokens.push_back(token);
516 for(
int i=0; i<m.nops(); i++) ret += m.op(i);
525 string now(
bool use_date) {
529 if(use_date) strftime(tmp,
sizeof(tmp),
"%Y-%m-%d %H:%M:%S",localtime(&timep) );
530 else strftime(tmp,
sizeof(tmp),
"%H:%M:%S",localtime(&timep) );
541 for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i) {
542 if(is_a<symbol>(*i)) ss.insert(*i);
545 for(
auto item : ss) sym_lst.append(item);
557 for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i) {
558 if(is_a<symbol>(*i)) ss.insert(*i);
562 for(
auto item : ss) sym_lst.append(item);
575 if( (fp = popen(cmd.c_str(),
"r")) != NULL) {
576 while(fgets(buf, 128, fp) != NULL) {
591 void garRead(
const string &garfn, map<string, ex> &resMap) {
596 for(
int i=0; i<ar.num_expressions(); i++) {
598 ex res = ar.unarchive_ex(name, i);
609 ex
garRead(
const string &garfn,
const char* key) {
614 auto res = ar.unarchive_ex(key);
628 auto c = ar.unarchive_ex(
"c");
629 auto res = ar.unarchive_ex(
"res");
630 if(c!=19790923)
throw Error(
"garRead: check faild for file: " + garfn);
639 void garWrite(
const string &garfn,
const map<string, ex> &resMap) {
641 for(
const auto & item : resMap) {
642 ar.archive_ex(item.second, item.first.c_str());
654 void garWrite(
const string &garfn,
const ex & res) {
656 ar.archive_ex(res,
"res");
657 ar.archive_ex(19790923,
"c");
670 ex
str2ex(
const string &expr, symtab stab) {
672 ex ret = par.
Read(expr);
683 ex ret = par.
Read(expr);
693 lst
str2lst(
const string &expr, symtab stab) {
695 ex ret = par.
Read(expr);
696 return ex_to<lst>(ret);
706 ex ret = par.
Read(expr);
707 return ex_to<lst>(ret);
739 int nc = ls.op(0).nops();
741 for(
int r=0; r<nr; r++) {
743 for(
int c=0; c<nc; c++) mat(r,c) = row.op(c);
776 ifstream ifs(filename);
777 string ostr((istreambuf_iterator<char>(ifs)), (istreambuf_iterator<char>()));
787 void str2file(
const string & ostr,
string filename) {
789 ofs.open(filename, ios::out);
797 int n = nstr.length();
799 strcpy(nbuff, nstr.c_str());
800 FILE * f = fmemopen(nbuff, n,
"r");
810 ifstream fs(filename);
813 while (std::getline(fs, line)) {
814 if(!skip_empty || line.size()>0) ovec.push_back(line);
845 void ex2file(
const ex & expr,
string filename) {
847 ofs.open(filename, ios::out);
858 void ex2file(
string filename,
const ex & expr) {
869 quadmath_snprintf(buffer,
sizeof buffer,
"%.36QG", num);
882 nss << num.evalf() << endl;
883 __float128 ret = strtoflt128(nss.str().c_str(), NULL);
894 return ex_to<numeric>(num).to_int();
903 return GiNaC::lst(ev.begin(), ev.end());
912 exvector ret(alst.begin(), alst.end());
913 return exvector(std::move(ret));
922 if(!is_a<add>(expr))
return lst{expr};
924 for(
auto item : expr) ret.append(item);
934 if(!is_a<mul>(expr))
return lst{expr};
936 for(
auto item : expr) ret.append(item);
948 for(
int i=bi; i<=ei; i++) ret.append(x(i));
968 ex
series_ex(ex
const & expr_in, ex
const &s0,
int sn0) {
970 return series_ex(expr_in.subs(s0==ss),ss,sn0).subs(ss==s0);
980 ex
series_ex(ex
const & expr_in, symbol
const &s0,
int sn0) {
982 if(!expr.has(s0))
return (sn0>=0 ? expr : 0);
985 expr.find(pow(s0,
w), sset);
987 exmap repl1st, repl2nd;
988 for(
auto pi : sset) {
989 auto sn = expand(pi.op(1));
990 if(!is_a<numeric>(sn)) {
993 if(!is_a<numeric>(item)) ex1 += item;
998 repl1st[pi] = sss * pow(s0, sn);
999 repl2nd[sss] = pow(s0, ex1);
1001 if(!(is_a<numeric>(sn) && ex_to<numeric>(sn).is_rational())) {
1002 cerr <<
"s = " << s0 << endl;
1003 cerr <<
"expr_in = " << expr_in << endl;
1004 throw Error(
"series_ex: Not rational sn = " +
ex2str(sn));
1006 sn_lcm = lcm(sn_lcm, ex_to<numeric>(sn).denom());
1008 if(expr.has(sqrt(s0))) sn_lcm = lcm(sn_lcm, numeric(2));
1009 expr = expr.subs(repl1st);
1012 if(!sn_lcm.is_integer())
throw Error(
"series_ex: Not integer with " +
ex2str(sn_lcm));
1013 if(sn_lcm<0) sn_lcm = numeric(0)-sn_lcm;
1014 int sn = sn0 * sn_lcm.to_int();
1015 expr = expr.subs(pow(s0,
w)==pow(s,
w*sn_lcm)).subs(sqrt(s0)==pow(s,sn_lcm/2)).subs(s0==pow(s,sn_lcm));
1019 for(
auto cv : cvs) {
1023 expr = cv.op(1) + pow(s,sn+exN+2)+pow(s,sn+exN+3);
1024 expr = expr.series(s,sn+exN);
1026 for(
int i=0; i<expr.nops(); i++) {
1027 if(is_order_function(expr.op(i))) {
1032 if(!is_order_function(ot)) {
1033 cerr <<
"expr = " << expr << endl;
1034 throw Error(
"series_ex: Not an Order term with " +
ex2str(ot));
1036 if(ot.op(0).degree(s)>sn) {
1037 expr = series_to_poly(expr);
1039 for(
int i=expr.ldegree(s); (i<=expr.degree(s) && i<=sn); i++) {
1040 ret += cv.op(0) * expr.coeff(s,i) * pow(s0,ex(i)/sn_lcm);
1047 if(!ok)
throw Error(
"series_ex seems not working!");
1049 ret = ret.subs(s==pow(s0,ex(1)/sn_lcm));
1050 ret = ret.subs(repl2nd);
1063 ex
diff_ex(ex
const expr, ex
const xp,
unsigned nth,
bool expand) {
1065 ex res = expr.subs(xp==s);
1066 if(!expand)
return res.diff(s,nth).subs(s==xp);
1069 for(
auto cv : cvs) {
1070 res += cv.op(0) * cv.op(1).diff(s,nth).subs(s==xp);
1078 if(n==1)
return co_epv;
1080 ex co = pow(co_epv.first,2);
1083 auto epv = co_epv.second;
1084 if(epv.size()>10000) cout <<
"Warning: power_expand_2: epv.size > 10000" << endl;
1085 for(
int i=0; i<epv.size(); i++) {
1086 pcmap[epv[i].first] += 2*co_epv.first*epv[i].second;
1087 pcmap[pow(epv[i].first,2)] += pow(epv[i].second,2);
1088 for(
int j=i+1; j<epv.size(); j++) pcmap[epv[i].first*epv[j].first] += 2*epv[i].second*epv[j].second;
1092 epv.reserve(pcmap.size());
1093 for(
auto kv : pcmap) {
1094 if(is_zero(kv.first) || is_zero(kv.second))
continue;
1095 epv.push_back(make_pair(kv.first, kv.second));
1097 return make_pair(co,
epvec_t(std::move(epv)));
1102 if((n % 2)==0)
return co_epvec_t(std::move(co_epv2));
1104 ex co = co_epv.first * co_epv2.first;
1106 for(
auto ep1 : co_epv.second) {
1107 pcmap[ep1.first] += ep1.second * co_epv2.first;
1108 for(
auto ep2 : co_epv2.second) pcmap[ep2.first * ep1.first] += ep1.second * ep2.second;
1110 for(
auto ep2 : co_epv2.second) pcmap[ep2.first] += ep2.second * co_epv.first;
1113 epv.reserve(pcmap.size());
1114 for(
auto kv : pcmap) {
1115 if(is_zero(kv.first) || is_zero(kv.second))
continue;
1116 epv.push_back(make_pair(kv.first, kv.second));
1118 return make_pair(co,
epvec_t(std::move(epv)));
1122 if(!has_func(expr_in)) {
1126 return make_pair(co,epv);
1127 }
if(is_a<add>(expr_in)) {
1130 for(
auto item : expr_in) {
1131 if(has_func(item)) {
1134 for(
auto ep : co_epv.second) pcmap[
ep.first] +=
ep.second;
1138 epv.reserve(pcmap.size());
1139 for(
auto kv : pcmap) {
1140 if(is_zero(kv.first) || is_zero(kv.second))
continue;
1141 epv.push_back(make_pair(kv.first, kv.second));
1143 return make_pair(co,
epvec_t(std::move(epv)));
1144 }
else if(is_a<mul>(expr_in)) {
1147 for(
auto item : expr_in) {
1148 if(!has_func(item)) {
1149 for(
auto &
ep : epv)
ep.second *= item;
1154 for(
auto ep2 : co_epv.second) {
1155 pcmap[ep2.first] += ep2.second * co;
1156 for(
auto ep1 : epv) pcmap[ep1.first * ep2.first] += ep1.second * ep2.second;
1158 for(
auto ep1 : epv) pcmap[ep1.first] += ep1.second * co_epv.first;
1161 epv.reserve(pcmap.size());
1162 for(
auto kv : pcmap) {
1163 if(is_zero(kv.first) || is_zero(kv.second))
continue;
1164 epv.push_back(make_pair(kv.first, kv.second));
1168 return make_pair(co,
epvec_t(std::move(epv)));
1169 }
else if(is_a<power>(expr_in) && expr_in.op(1).info(info_flags::nonnegint)) {
1170 int n = ex_to<numeric>(expr_in.op(1)).to_int();
1176 epv.push_back(make_pair(expr_in, 1));
1177 return make_pair(co,
epvec_t(std::move(epv)));
1179 throw Error(
"inner_expand_collect unexpected region reached.");
1188 ex
expand_ex(ex
const &expr_in, std::function<
bool(
const ex &)> has_func) {
1190 ex ret = co_epv.first;
1191 for(
auto ep : co_epv.second) ret +=
ep.second *
ep.first;
1202 ex
collect_ex(ex
const &expr_in, std::function<
bool(
const ex &)> has_func,
int opt) {
1205 res_vec.reserve(cvs.nops());
1206 for(
auto cv : cvs) {
1209 if(cc.is_zero() || vv.is_zero())
continue;
1210 res_vec.push_back(cc * vv);
1212 return add(res_vec);
1222 lst
collect_lst(ex
const &expr_in, std::function<
bool(
const ex &)> has_func,
int opt) {
1224 ex cf = co_epv.first;
1226 if(!is_zero(cf)) co_epv.second.push_back(make_pair(1,cf));
1227 for(
auto ep : co_epv.second) {
1235 if(!is_zero(cc)) res_lst.append(lst{cc, vv});
1248 expr.find(zeta(
w), zs);
1249 expr.find(zeta(
w,
w), zs);
1252 for(
auto zi : zs) repl.append(zi==
NN(zi));
1253 return expr.subs(repl);
1257 static symbol sPi(
"dPi"), sEuler(
"dEuler"), siEpsilon(
"diEpsilon");
1258 lst repl = lst{Pi==sPi, Euler==sEuler,
iEpsilon==siEpsilon};
1259 return EvalF(expr.subs(repl));
1262 static symbol sPi(
"qPi"), sEuler(
"qEuler"), siEpsilon(
"qiEpsilon");
1263 lst repl = lst{Pi==sPi, Euler==sEuler,
iEpsilon==siEpsilon};
1264 return EvalF(expr.subs(repl));
1267 static symbol sPi(
"mpPi"), sEuler(
"mpEuler"), siEpsilon(
"mpiEpsilon");
1268 lst repl = lst{Pi==sPi, Euler==sEuler,
iEpsilon==siEpsilon};
1269 return EvalF(expr.subs(repl));
1278 ex
NN(ex expr,
int digits) {
1280 auto nexpr = evalf(expr);
1291 auto tmp = expand(expr);
1292 if(tmp.is_zero())
return false;
1294 if(is_a<add>(tmp)) {
1295 for(
auto item : tmp) {
1296 auto nit =
NN(item.subs(x(
w)==1).normal());
1297 if(!is_a<numeric>(nit) || nit<0)
return false;
1301 auto ntmp =
NN(tmp.subs(x(
w)==1).normal());
1302 ret = (is_a<numeric>(ntmp) && ntmp>0);
1313 if(is_zero(expr))
return 0;
1315 else if(
xPositive(ex(0)-expr))
return -1;
1325 auto tmp = ex_to<lst>(ex_in);
1336 auto tmp = ex_to<lst>(ex_in);
1346 auto tmp = ex_to<lst>(ex_in);
1356 auto tmp = ex_to<lst>(ex_in);
1368 auto tmp = ex_to<lst>(ex_in.op(index));
1370 ex_in.let_op(index) = tmp;
1380 auto tmp = ex_to<lst>(ex_in.op(index));
1382 ex_in.let_op(index) = tmp;
1393 auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1395 ex_in.let_op(index1).let_op(index2) = tmp;
1406 auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1408 ex_in.let_op(index1).let_op(index2) = tmp;
1418 auto tmp = ex_to<lst>(ex_in.op(index));
1420 ex_in.let_op(index) = tmp;
1430 auto tmp = ex_to<lst>(ex_in.op(index));
1432 ex_in.let_op(index) = tmp;
1443 auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1445 ex_in.let_op(index1).let_op(index2) = tmp;
1456 auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1458 ex_in.let_op(index1).let_op(index2) = tmp;
1467 auto tmp = ex_to<lst>(ex_in.op(index));
1469 ex_in.let_op(index) = tmp;
1478 auto tmp = ex_to<lst>(ex_in.op(index));
1480 ex_in.let_op(index) = tmp;
1490 auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1492 ex_in.let_op(index1).let_op(index2) = tmp;
1502 auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1504 ex_in.let_op(index1).let_op(index2) = tmp;
1513 auto tmp = ex_to<lst>(ex_in.op(index));
1515 ex_in.let_op(index) = tmp;
1524 auto tmp = ex_to<lst>(ex_in.op(index));
1526 ex_in.let_op(index) = tmp;
1536 auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1538 ex_in.let_op(index1).let_op(index2) = tmp;
1548 auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1550 ex_in.let_op(index1).let_op(index2) = tmp;
1560 void let_op(ex &ex_in,
int index1,
int index2,
const ex item) {
1561 ex_in.let_op(index1).let_op(index2) = item;
1571 void let_op(lst &ex_in,
int index1,
int index2,
const ex item) {
1572 ex_in.let_op(index1).let_op(index2) = item;
1583 void let_op(ex &ex_in,
int index1,
int index2,
int index3,
const ex item) {
1584 ex_in.let_op(index1).let_op(index2).let_op(index3) = item;
1595 void let_op(lst &ex_in,
int index1,
int index2,
int index3,
const ex item) {
1596 ex_in.let_op(index1).let_op(index2).let_op(index3) = item;
1606 ex
get_op(
const ex ex_in,
int index1,
int index2) {
1607 return ex_in.op(index1).op(index2);
1617 ex
get_op(
const lst ex_in,
int index1,
int index2) {
1618 return ex_in.op(index1).op(index2);
1629 ex
get_op(
const ex ex_in,
int index1,
int index2,
int index3) {
1630 return ex_in.op(index1).op(index2).op(index3);
1641 ex
get_op(
const lst ex_in,
int index1,
int index2,
int index3) {
1642 return ex_in.op(index1).op(index2).op(index3);
1663 if(!is_a<XIntegral>(other))
throw Error(
"XIntegral::compare_same_type");
1684 else if(i==2)
return Deltas;
1685 else throw Error(
"XIntegral::op, It is required that i<3.");
1688 ensure_if_modifiable();
1691 else if(i==2)
return Deltas;
1692 else throw Error(
"XIntegral::let_op, It is required that i<3.");
1696 inherited::archive(n);
1699 n.add_ex(
"Deltas",
Deltas);
1703 inherited::read_archive(n);
1706 n.find_ex(
"Deltas",
Deltas);
1712 if(fed.nops()>2)
Deltas = fed.op(2);
1722 for(
int i=0; i<ps.nops(); i++) {
1723 if(is_zero(ns.op(i)))
continue;
1726 for(
auto li : loops) {
1739 }
else if(is_a<numeric>(ns.op(i)) && (ns.op(i)<0)) {
1740 throw Error(
"XIntegral, negative powers not supported yet.");
1744 for(
auto li : loops) {
1745 if(!is_a<numeric>(cpi.coeff(li,2)))
continue;
1746 numeric nm = ex_to<numeric>(cpi.coeff(li,2));
1747 if(nm.is_zero())
continue;
1748 sgn = nm>0 ? -1 : 1;
1753 if(!is_a<numeric>(cpi.coeff(
iEpsilon)))
continue;
1754 numeric nm = ex_to<numeric>(cpi.coeff(
iEpsilon));
1755 if(!nm.is_zero()) sgn = nm>0 ? -1 : 1;
1760 if(is_a<numeric>(cpi) && ex_to<numeric>(cpi)>0) sgn = -1;
1765 if(sgn==-1) pre *= exp(I * Pi * ns.op(i));
1770 if(!is_zero(ns.op(i)-1)) {
1772 exps.append(ns.op(i)-1);
1774 pre /= tgamma(ns.op(i));
1778 ex dl2 = loops.nops()*(4-2*
ep)/2;
1779 pre *= tgamma(
a-dl2) * pow(I,loops.nops()) * pow(Pi,dl2) * pow(2*Pi,loops.nops()*(2*
ep-4));
1783 for(
int i=0; i<loops.nops(); i++) {
1784 ex t2 = rem.coeff(loops.op(i),2);
1785 ex t1 = rem.coeff(loops.op(i),1);
1786 ex t0 = rem.coeff(loops.op(i),0);
1793 rem = (t0 - pow(t1,2)/(4*t2)).expand();
1797 rem = (rem * u).normal();
1800 exps.prepend(dl2-
a);
1802 exps.prepend(
a-dl2-(4-2*
ep)/2);
1803 auto num_den = numer_denom(pre);
1804 funs.append(num_den.op(0));
1806 if(!is_zero(num_den.op(1))) {
1807 funs.append(num_den.op(1));
1821 return omp_get_num_procs();
1833 return nd.op(0)/nd.op(1);
1843 return collect_common_factors(expr);
1855 if(opt==
o_none)
return expr;
1868 if(opt==1)
return expand(expr);
1869 else return expand(expr);
1873 exvector ev_fermat_pern(
const exvector & ev,
int np) {
1874 if(ev.size()<np)
return ev;
1880 for(
int i=0; i<np; i++) {
1888 return exvector(std::move(evs));
1891 ex normal_fermat_pern(
const ex & e,
int np) {
1893 exvector evs(e.begin(), e.end());
1894 while(evs.size()>np) evs = ev_fermat_pern(evs,np);
1907 if(opt<0)
return normal_fermat_pern(expr,-opt);
1908 else if(opt==
o_normal)
return normal(expr);
1922 ex
exnd(
const ex & expr,
int opt) {
1925 else return numer_denom(expr);
1929 static map<pid_t, Form> form_map;
1930 auto pid = getpid();
1931 if((form_map.find(pid)==form_map.end())) {
1933 ss <<
"AutoDeclare Symbols fv;" << endl;
1934 form_map[pid].Init(
"form");
1935 form_map[pid].Execute(ss.str());
1937 Form &fprc = form_map[pid];
1939 auto expr_in = expr;
1941 expr_in = expr_in.to_polynomial(map_rat);
1944 for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
1946 if(is_a<symbol>(e))
vs.insert(e);
1954 string cvv =
"fv"+to_string(cid);
1956 st[cvv] = ss.subs(map_rat);
1959 expr_in = expr_in.subs(s2v);
1960 oss <<
"Local ff = " << expr_in <<
";" << endl;
1961 oss <<
".sort" << endl;
1963 auto ostr = fprc.
Execute(oss.str(),
"ff");
1965 ex ret = fp.
Read(ostr);
1967 }
catch(
Error& err) {
1968 cout << oss.str() << endl;
1969 form_map.erase(pid);
1975 if(is_a<mul>(expr)) {
1980 static map<pid_t, Form> form_map;
1981 auto pid = getpid();
1982 if((form_map.find(pid)==form_map.end())) {
1984 ss <<
"AutoDeclare Symbols fv;" << endl;
1985 form_map[pid].Init(
"form");
1986 form_map[pid].Execute(ss.str());
1988 Form &fprc = form_map[pid];
1990 auto expr_in = expr;
1992 expr_in = expr_in.to_polynomial(map_rat);
1995 for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
1997 if(is_a<symbol>(e))
vs.insert(e);
2005 string cvv =
"fv"+to_string(cid);
2007 st[cvv] = ss.subs(map_rat);
2010 expr_in = expr_in.subs(s2v);
2011 oss <<
"Local ff = " << expr_in <<
";" << endl;
2012 oss <<
"Factorize ff;" << endl;
2013 oss <<
".sort" << endl;
2015 auto ostr = fprc.
Execute(oss.str(),
"ff");
2021 ex ret = fp.
Read(ostr);
2024 for(
auto item :
add2lst(ret)) res *= item.subs(sfac==1);
2026 }
catch(
Error& err) {
2027 cout << oss.str() << endl;
2028 form_map.erase(pid);
2050 GINAC_IMPLEMENT_PRINT_CONTEXT(MMAFormat, print_dflt)
2071 int nr = mat.rows();
2072 int nc = mat.cols();
2073 for(
int r=0; r<nr; r++) {
2075 for(
int c=0; c<nc; c++) {
2076 mat(r,c).print(*
this);
2077 if(c+1!=nc) s <<
",";
2080 if(r+1!=nr) s <<
",";
2088 auto vend = e.end();
2089 if (i==vend) { s <<
"{}";
return *
this; }
2103 auto send = e.end();
2104 if (i==send) { s <<
"{}";
return *
this; }
2118 auto mend = e.end();
2119 if (i==mend) { s <<
"{}";
return *
this; }
2122 i->first.print(*
this);
2124 i->second.print(*
this);
2135 ar.archive_ex(exv.size(),
"size");
2136 for(
int i=0; i<exv.size(); i++) {
2137 ar.archive_ex(exv[i], to_string(i).c_str());
2139 ofstream out(garfn);
2150 map<string, ex> dict;
2151 for(
int i=0; i<ar.num_expressions(); i++) {
2153 ex res = ar.unarchive_ex(name, i);
2157 auto size = ex_to<numeric>(dict[
"size"]).to_int();
2158 if(exv.size()>0 && size != exv.size())
throw Error(
"garRead: exvector size>0 & not match!");
2159 if(exv.size()<1) exv.resize(size);
2160 for(
int i=0; i<size; i++) exv[i] = dict[to_string(i)];
2164 auto cvs_vec =
GiNaC_Parallel(exv.size(), [&exv,pats](
int idx)->ex {
2165 return collect_lst(exv[idx], pats);
2168 for(
auto cvs : cvs_vec)
for(
auto cv : cvs) res_map[cv.op(1)] += cv.op(0);
2170 for(
auto kv : res_map) res_vec.push_back(lst{kv.second, kv.first});
2171 res_vec =
GiNaC_Parallel(res_vec.size(), [&res_vec,opt](
int idx)->ex {
2172 return exnormal(res_vec[idx].op(0), opt) * res_vec[idx].op(1);
2174 return add(res_vec);
2178 auto cvs_vec =
GiNaC_Parallel(exv.size(), [&exv,pats](
int idx)->ex {
2179 return collect_lst(exv[idx], pats);
2182 for(
auto cvs : cvs_vec)
for(
auto cv : cvs) res_map[cv.op(1)] += cv.op(0);
2184 for(
auto kv : res_map) res_vec.push_back(lst{kv.second, kv.first});
2185 res_vec =
GiNaC_Parallel(res_vec.size(), [&res_vec,opt](
int idx)->ex {
2186 return exnormal(res_vec[idx].op(0),opt) * res_vec[idx].op(1);
2188 return add(res_vec);
2192 auto cvs_vec =
GiNaC_Parallel(exv.size(), [&exv,pats](
int idx)->ex {
2193 return collect_lst(exv[idx], pats);
2196 for(
auto cvs : cvs_vec)
for(
auto cv : cvs) res_map[cv.op(1)] += cv.op(0);
2198 for(
auto kv : res_map) res_vec.push_back(lst{kv.second, kv.first});
2199 res_vec =
GiNaC_Parallel(res_vec.size(), [&res_vec,opt](
int idx)->ex {
2200 return exnormal(res_vec[idx].op(0),opt) * res_vec[idx].op(1);
2202 return add(res_vec);
2206 if(!is_a<add>(e))
throw Error(
"add_collect_normal: input is NOT a add class.");
2207 exvector exv(e.begin(), e.end());
2212 if(!is_a<add>(e))
throw Error(
"add_collect_normal: input is NOT a add class.");
2213 exvector exv(e.begin(), e.end());
2218 if(!is_a<add>(e))
throw Error(
"add_collect_normal: input is NOT a add class.");
2219 exvector exv(e.begin(), e.end());
2224 for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i)
if(is_a<wildcard>(*i))
return true;
2230 for(
auto kv : r)
if(is_a<mul>(kv.first) && !
has_w(kv.first)) repl[
w*kv.first]=
w*kv.second;
2235 for(
auto item : r)
if(is_a<mul>(item) && !
has_w(item)) repl.append(
w*item.op(0)==
w*item.op(1));
2242 ar.archive_ex(e,
"e");
2247 for(
auto const & e : es) ar.archive_ex(e,
"e");
2252 ar.archive_ex(e1,
"e");
2253 ar.archive_ex(e2,
"e");
2256 void ReShare(
const ex & e1,
const ex & e2,
const ex & e3) {
2258 ar.archive_ex(e1,
"e");
2259 ar.archive_ex(e2,
"e");
2260 ar.archive_ex(e3,
"e");
2265 for(
auto & e : ev) ar.archive_ex(e,
"e");
2268 void ReShare(
const exvector & ev1,
const exvector & ev2) {
2270 for(
auto & e : ev1) ar.archive_ex(e,
"e");
2271 for(
auto & e : ev2) ar.archive_ex(e,
"e");
2275 auto v = cln::the<cln::cl_I>(ex_to<numeric>(n).to_cl_N());
2276 return numeric(cln::nextprobprime(v));
2279 auto v = cln::the<cln::cl_I>(numeric(n).to_cl_N());
2280 return numeric(cln::nextprobprime(v));
2286 if(is_a<numeric>(e)) {
2287 auto ne = ex_to<numeric>(e);
2288 if(ne.is_crational())
return e;
2289 auto zz = ne.to_cl_N();
2290 auto re = cln::rationalize(cln::realpart(zz));
2291 auto im = cln::rationalize(cln::imagpart(zz));
2292 return numeric(cln::complex(re,im));
2293 }
else return e.map(
self);
2310 cln::default_float_format = cln::float_format(prec);
2318 cln::default_float_format = prec;
2324 return cln::default_float_format;
2327 void arg2map(
int argc,
char** argv,
const char *optstring, std::map<char,std::string> & kv) {
2328 auto o_opterr = opterr;
2331 int opt = getopt(argc, argv, optstring);
2334 if(optarg) kv[opt] = optarg;
2344 for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i)
if(is_a<symbol>(*i))
return true;
unsigned golden_ratio_hash(uintptr_t n)
#define IMPLEMENT_HAS(classname)
#define DEFAULT_CTOR(classname)
#define IMPLEMENT_ALL(classname)
class used to wrap error message
const char * what() const
Error(string _msg)
Error constructor,.
class to wrap map_function of GiNaC
class to parse for string or file, helpful with Symbol class
ex Read(const string &instr, bool s2S=true)
class extended to GiNaC symbol class, represent a positive symbol
ex series(const relational &s, int order, unsigned options=0) const override
void archive(archive_node &n) const override
Symbol archive.
ex imag_part() const override
Symbol * duplicate() const override
bool is_equal_same_type(const basic &other) const override
ex conjugate() const override
static ex set_all(const ex &expr)
const char * class_name() const override
void read_archive(const archive_node &n) override
Symbol read_archive.
static GiNaC::registered_class_info & get_class_info_static()
ex derivative(const symbol &s) const override
unsigned calchash() const override
ex evalf() const override
void set(const ex &v) const
static std::map< std::string, ex > Table
void accept(GiNaC::visitor &v) const override
ex real_part() const override
int compare_same_type(const GiNaC::basic &other) const override
const GiNaC::registered_class_info & get_class_info() const override
XIntegral Class, preface to SecDec.
ex op(size_t i) const override
int compare_same_type(const GiNaC::basic &other) const override
void accept(GiNaC::visitor &v) const override
XIntegral * duplicate() const override
void archive(archive_node &n) const override
ex & let_op(size_t i) override
void print(const print_dflt &c, unsigned level=0) const
void read_archive(const archive_node &n) override
const char * class_name() const override
static GiNaC::registered_class_info & get_class_info_static()
const GiNaC::registered_class_info & get_class_info() const override
size_t nops() const override
class extended to GiNaC symbol class, pure imaginary symbol
bool is_equal_same_type(const basic &other) const override
ex evalf() const override
static std::map< std::string, ex > Table
void archive(archive_node &n) const override
Symbol archive.
iSymbol * duplicate() const override
ex real_part() const override
ex derivative(const symbol &s) const override
static GiNaC::registered_class_info & get_class_info_static()
ex conjugate() const override
const char * class_name() const override
ex series(const relational &s, int order, unsigned options=0) const override
void accept(GiNaC::visitor &v) const override
ex imag_part() const override
const GiNaC::registered_class_info & get_class_info() const override
int compare_same_type(const GiNaC::basic &other) const override
void read_archive(const archive_node &n) override
Symbol read_archive.
unsigned calchash() const override
static ex sum(lst m)
sum all elements in the list
ex expand_ex(ex const &expr_in, std::function< bool(const ex &)> has_func)
the expand like Mathematica
const char * Color_HighLight
ex collect_factors(const ex &expr)
a wrapper for collect_common_factors, catch errors
void arg2map(int argc, char **argv, const char *optstring, std::map< char, std::string > &kv)
__float128 ex2q(ex num)
ex of numeric to __float128
ex EvalF(ex expr)
the nuerical evaluation, Digits=100 will be used
lst str2lst(const string &expr, symtab stab)
convert string to the lst, using Parser internally
ex exnd(const ex &expr, int opt)
num_den a expression
ex form_factor(const ex &expr, bool nd=true)
lst gather_symbols(const ex &e)
get all symbols from input expression
void let_op_prepend(ex &ex_in, const ex item)
preppend item into expression
lst xlst(int bi, int ei)
return a lst: x(bi), x(bi+1), ..., x(ei)
int GiNaC_Parallel_Process
void let_op(ex &ex_in, int index1, int index2, const ex item)
update index1-th.index2-th of expression with item
std::initializer_list< ex > init_list
ex q2ex(__float128 num)
__float128 to ex
bool xPositive(ex const expr)
check the expr is xPositive, i.e., each x-monomial item is postive
void ex2file(const ex &expr, string filename)
export expression file
ex normal_fermat(const ex &expr, bool dfactor)
return the normalizatied expression, using fermat_numer_denom
ex NN(ex expr, int digits)
the nuerical evaluation
map< string, bool > GiNaC_Parallel_RM
ex factor_flint(const ex &e, bool nd=true)
ex collect_ex(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica
void garRead(const string &garfn, map< string, ex > &resMap)
garRead from file, and output in a map
map< string, int > GiNaC_Parallel_NB
void ReShare(const ex &e)
bool file_exists(string fn)
ex factor_form(const ex &expr, bool nd)
factorize a expression using FORM
ex series_ex(ex const &expr_in, ex const &s0, int sn0)
the series like Mathematica, include s^n
ex str2ex(const string &expr, symtab stab)
convert string to ex expression, using Parser internally
ex form_eval(const ex &expr)
string now(bool use_date)
date/time string
vector< pair< ex, ex > > epvec_t
ex numer_denom_fermat(const ex &expr, bool dfactor=false)
return the numerator and denominator after normalization
string RunOS(const string &cmd)
run symtem command and return the output as string
lst add2lst(const ex &expr)
convert add to lst
ex fermat_normal(const ex &expr, bool dfactor=false)
int CpuCores()
return the cpu cores using OpenMP
map< string, int > GiNaC_Parallel_NP
map< string, int > GiNaC_Parallel_Verb
map< string, bool > GiNaC_Parallel_ReWR
exvector lst2vec(const lst &alst)
convert lst to exvector
ex get_op(const ex ex_in, int index1, int index2)
return index1-th.index2-th of expression
vector< string > file2strvec(string filename, bool skip_empty)
read file content to string vector
ex file2ex(string filename)
read file content to ex
ex exnormal(const ex &expr, int opt)
normalize a expression
const int o_normal_fermat
ex nextprime(const ex &n)
ex normal_flint(const ex &expr, int opt=o_flint)
ex inner_factor_form(const ex &expr)
lst mul2lst(const ex &expr)
convert mul to lst
string file2str(string filename)
read file content to string
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
ex Rationalize(const ex &e, int dn)
void string_replace_all(string &str, const string &from, const string &to)
ex exexpand(const ex &expr, int opt)
factorize a expression
void garWrite(const string &garfn, const map< string, ex > &resMap)
garWrite to write the string-key map to the archive
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
void let_op_remove_first(ex &ex_in)
remove first from expression
ex exfactor(const ex &expr, int opt)
factorize a expression
matrix lst2mat(const lst &ls)
convert 2Dim lst to matrix
bool dir_exists(string dir)
pair< ex, epvec_t > inner_expand_collect(ex const &expr_in, std::function< bool(const ex &)> has_func, int depth=0)
void str2file(const string &ostr, string filename)
export string to a file
co_epvec_t power_expand_2(const co_epvec_t &co_epv, int n)
void set_precision(long prec, bool push)
ex add_collect_normal(const exvector &exv, ex const &pats, int opt)
map< string, string > GiNaC_Parallel_PRE
int xSign(ex const expr)
the always sign for expr
pair< ex, epvec_t > co_epvec_t
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
void let_op_remove_last(ex &ex_in)
remove last from expression
std::vector< std::string > split(const std::string &s, char delimiter)
split the string into serveral part, separated by the delimiter
bool has_symbol(const ex &e)
std::stack< cln::float_format_t > cln_prec_stack
std::stack< long > digits_stack
ex numer_fermat(const ex &expr)
ex fermat_numer_denom(const ex &expr, bool dfactor=false)
int ex2int(ex num)
ex to integer
ex subs(const ex &e, init_list sl)