HepLib
BASIC.cpp
Go to the documentation of this file.
1 
6 #include "BASIC.h"
7 #include "cln/cln.h"
8 
9 inline unsigned golden_ratio_hash(uintptr_t n) {
10  return n * UINT64_C(0x4f1bbcdd);
11 }
12 
13 extern "C" {
14  #include <quadmath.h>
15 }
16 
17 namespace HepLib {
18 
23  Error::Error(string _msg) : msg(_msg) { }
24  const char * Error::what() const throw () {
25  return msg.c_str();
26  }
27 
28  /*-----------------------------------------------------*/
29  // Symbol
30  /*-----------------------------------------------------*/
31 
32  namespace {
33  const symbol & get_symbol(const string & s) {
34  static map<string, symbol> dict;
35  string key = s;
36  if (dict.find(key) == dict.end()) dict[key] = symbol(s,0);
37  return dict[key];
38  }
39 
40  const Symbol & get_Symbol(const string & s) {
41  static map<string, Symbol> dict;
42  string key = s;
43  if (dict.find(key) == dict.end()) dict[key] = Symbol(s);
44  return dict[key];
45  }
46 
47  const iSymbol & get_iSymbol(const string & s) {
48  static map<string, iSymbol> dict;
49  string key = s;
50  if (dict.find(key) == dict.end()) dict[key] = iSymbol(s);
51  return dict[key];
52  }
53  }
54 
55  //DEFAULT_CTOR(Symbol)
56  IMPLEMENT_HAS(Symbol)
57  IMPLEMENT_ALL(Symbol)
58  //GINAC_IMPLEMENT_REGISTERED_CLASS(Symbol, symbol)
59  GiNaC::registered_class_info & Symbol::get_class_info_static() { return reg_info; }
61  Symbol * Symbol::duplicate() const { Symbol * bp = new Symbol(*this); bp->setflag(GiNaC::status_flags::dynallocated); return bp; }
62  void Symbol::accept(GiNaC::visitor & v) const { if (visitor *p = dynamic_cast<visitor *>(&v)) p->visit(*this); else inherited::accept(v); }
63  const GiNaC::registered_class_info &Symbol::get_class_info() const { return get_class_info_static(); }
64  GiNaC::registered_class_info &Symbol::get_class_info() { return get_class_info_static(); }
65  const char *Symbol::class_name() const { return get_class_info_static().options.get_name(); }
66  //GINAC_IMPLEMENT_REGISTERED_CLASS END
67 
72  Symbol::Symbol(const string &s) : symbol(get_symbol(s)) { Table[s]=*this; }
73  Symbol::Symbol() : symbol(get_symbol("")) { Table[""]=*this; }
74 
75  int Symbol::compare_same_type(const basic &other) const {
76  if(!is_a<Symbol>(other)) throw Error("Symbol::compare_same_type");;
77  const Symbol &o = static_cast<const Symbol &>(other);
78  int ret = get_name().compare(o.get_name());
79  if(ret==0) return 0;
80  else if(ret<0) return -1;
81  else return 1;
82  }
83 
88  void Symbol::archive(archive_node & n) const {
89  inherited::archive(n);
90  }
91 
96  void Symbol::read_archive(const archive_node& n) {
97  inherited::read_archive(n);
98  }
99 
100  ex Symbol::eval() const { return *this; }
101  ex Symbol::evalf() const { return *this; }
102  ex Symbol::conjugate() const { return *this; }
103  ex Symbol::real_part() const { return *this; }
104  ex Symbol::imag_part() const { return 0; }
105  unsigned Symbol::calchash() const {
106  static std::hash<std::string> hs;
107  unsigned seed = hs(get_name()+typeid(*this).name());
108  hashvalue = golden_ratio_hash(seed);
109  setflag(status_flags::hash_calculated);
110  return hashvalue;
111  }
112 
113  bool Symbol::is_equal_same_type(const basic & other) const {
114  if(!is_a<Symbol>(other)) throw Error("Symbol::is_equal_same_type");
115  const Symbol *o = static_cast<const Symbol *>(&other);
116  return get_name()==o->get_name(); // should be the same as symbol class
117  }
118 
119  ex Symbol::derivative(const symbol & s) const {
120  if(!is_exactly_a<Symbol>(s)) return 0;
121  if(is_equal_same_type(s)) return 1;
122  else return 0;
123  }
124 
125  ex Symbol::series(const relational & r, int order, unsigned options) const {
126  epvector seq;
127  const ex point = r.rhs();
128 
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));
135  }
136 
137  void Symbol::set_name(string n) { throw Error("Symbol can not reset the name!"); }
138 
139  void Symbol::set(const ex & v) const { vmap[*this] = v; }
140  void Symbol::unset() const { vmap.erase(*this); }
141 
142  void Symbol::set(const Symbol & s, const ex & v) { vmap[s] = v; }
143  void Symbol::set(const string & str, const ex & v) { vmap[Symbol(str)] = v; }
144  void Symbol::unset(const Symbol &s) { vmap.erase(s); }
145  void Symbol::unset(const string &str) { vmap.erase(Symbol(str)); }
146  void Symbol::unset_all() { vmap.clear(); }
147  ex Symbol::set_all(const ex & expr) {
148  ex v1 = expr;
149  while(true) {
150  ex v2 = v1.subs(vmap);
151  if(v2.is_equal(v1)) break;
152  v1 = v2;
153  }
154  return v1;
155  }
156 
157  /*-----------------------------------------------------*/
158  // iSymbol
159  /*-----------------------------------------------------*/
160 
161  //DEFAULT_CTOR(iSymbol)
163  //GINAC_IMPLEMENT_REGISTERED_CLASS(iSymbol, symbol)
164  GiNaC::registered_class_info & iSymbol::get_class_info_static() { return reg_info; }
166  iSymbol * iSymbol::duplicate() const { iSymbol * bp = new iSymbol(*this); bp->setflag(GiNaC::status_flags::dynallocated); return bp; }
167  void iSymbol::accept(GiNaC::visitor & v) const { if (visitor *p = dynamic_cast<visitor *>(&v)) p->visit(*this); else inherited::accept(v); }
168  const GiNaC::registered_class_info &iSymbol::get_class_info() const { return get_class_info_static(); }
169  GiNaC::registered_class_info &iSymbol::get_class_info() { return get_class_info_static(); }
170  const char *iSymbol::class_name() const { return get_class_info_static().options.get_name(); }
171  //GINAC_IMPLEMENT_REGISTERED_CLASS END
172 
177  iSymbol::iSymbol(const string &s) : symbol(get_symbol(s)) { Table[s]=*this; }
178  iSymbol::iSymbol() : symbol(get_symbol("")) { Table[""]=*this; }
179  int iSymbol::compare_same_type(const basic &other) const {
180  if(!is_a<iSymbol>(other)) throw Error("iSymbol::compare_same_type");
181  const iSymbol &o = static_cast<const iSymbol &>(other);
182  int ret = get_name().compare(o.get_name());
183  if(ret==0) return 0;
184  else if(ret<0) return -1;
185  else return 1;
186  }
187 
188  bool iSymbol::is_equal_same_type(const basic & other) const {
189  if(!is_a<iSymbol>(other)) throw Error("iSymbol::is_equal_same_type");
190  const iSymbol *o = static_cast<const iSymbol *>(&other);
191  return serial==o->serial; // should be the same as symbol class
192  }
193 
194  ex iSymbol::derivative(const symbol & s) const {
195  if(!is_exactly_a<iSymbol>(s)) return 0;
196  if(is_equal_same_type(s)) return 1;
197  else return 0;
198  }
199 
200  ex iSymbol::series(const relational & r, int order, unsigned options) const {
201  epvector seq;
202  const ex point = r.rhs();
203 
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));
210  }
211 
212  void iSymbol::set_name(string n) { throw Error("iSymbol can not reset the name!"); }
213 
214  unsigned iSymbol::calchash() const {
215  hashvalue = get_symbol(get_name()).gethash();
216  setflag(status_flags::hash_calculated);
217  return hashvalue;
218  }
219 
224  void iSymbol::archive(archive_node & n) const {
225  inherited::archive(n);
226  }
227 
232  void iSymbol::read_archive(const archive_node& n) {
233  inherited::read_archive(n);
234  }
235 
236  ex iSymbol::eval() const { return *this; }
237  ex iSymbol::evalf() const { return *this; }
238  ex iSymbol::conjugate() const { return (*this)*(-1); }
239  ex iSymbol::real_part() const { return ex(0); }
240  ex iSymbol::imag_part() const { return (-I)*(*this); }
241 
242  /*-----------------------------------------------------*/
243  // Global varibales
244  /*-----------------------------------------------------*/
245 
246  const char* ErrColor = RED;
247  const char* WarnColor = MAGENTA;
248  const char* Color_HighLight = WHITE;
249 
250  static int GiNaC_Parallel_Level = 0; // for the internal usage only
260  exvector GiNaC_Parallel(int ntotal, std::function<ex(int)> f, const string & key) {
261  if(ntotal<1) return exvector();
262  int ec = 0;
263  int nactive = 0;
264  string exstr = "";
265  int fork_retried = 0;
266  int nproc = GiNaC_Parallel_Process;
267  if(GiNaC_Parallel_NP.find(key)!=GiNaC_Parallel_NP.end()) nproc = GiNaC_Parallel_NP[key];
268 
269  bool rm = true;
270  if(GiNaC_Parallel_RM.find(key)!=GiNaC_Parallel_RM.end()) rm = GiNaC_Parallel_RM[key];
271  string pre = PRE;
272  if(GiNaC_Parallel_PRE.find(key)!=GiNaC_Parallel_PRE.end()) pre = GiNaC_Parallel_PRE[key];
273 
274  int verb = Verbose;
275  if(GiNaC_Parallel_Verb.find(key)!=GiNaC_Parallel_Verb.end()) verb = GiNaC_Parallel_Verb[key];
276 
277  // nproc=0, non-parallel
278  if(nproc==0) {
279  exvector ovec;
280  for(int i=0; i<ntotal; i++) ovec.push_back(f(i));
281  return exvector(std::move(ovec));
282  }
283 
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();
287 
288  int rc;
289  auto ppid = getpid();
290  ostringstream cmd;
291  cmd << "mkdir -p " << ppid;
292  if(!dir_exists(to_string(ppid))) rc = system(cmd.str().c_str());
293 
294  int nbatch = GiNaC_Parallel_Batch;
295  if(GiNaC_Parallel_NB.find(key)!=GiNaC_Parallel_NB.end()) nbatch = GiNaC_Parallel_NB[key];
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);
300 
301  bool nst = (GiNaC_Parallel_Level>0);
302  pid_t npid=0, pgid;
303  if(getpgid(0)!=ppid) { // not group leader
304  if(nst) {
305  npid = fork();
306  if(npid < 0) throw Error("GiNaC_Parallel: Error (1) @ fork()");
307  if(npid!=0) goto wait_label;
308  } // else - for case in a shell script
309  if(setpgid(0,0)) {
310  if(setpgid(0,0)) throw Error("GiNaC_Parallel: setpgid(0,0) Failed.");
311  }
312  }
313 
314  pgid = getpid(); // should be groud leader @ here
315  while(fork_retried<5) {
316  ec = 0;
317  exstr = "";
318  nactive = 0;
319  for(int bi=0; bi<btotal; bi++) {
320  restart: ;
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 ";
325  if(key != "") cout << Color_HighLight << key << RESET << " ";
326  cout << Color_HighLight << nbatch << "x" << RESET << "[" << (bi+1) << "/" << btotal << "] @ " << now(false) << exstr << flush;
327  }
328 
329  if(fork_retried>0) { // skip when bi.*.gar exists
330  ostringstream garfn;
331  if(key == "") garfn << ppid << "/" << bi << ".gar";
332  else garfn << ppid << "/" << bi << "." << key << ".gar";
333  if(file_exists(garfn.str())) continue;
334  }
335 
336  auto pid = fork();
337  if(pid < 0) {
338  if(getpid()!=pgid) exit(0); // make sure exit child process
339  ec++;
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--;
343  goto restart; // parent process goes next cycle
344  }
345  if(pid==0 && setpgid(0, pgid)!=0) {
346  if(setpgid(0, pgid)) throw Error("GiNaC_Parallel: setpgid(0, pgid) Failed.");
347  }
348  if(pid>0) {
349  nactive++;
350  if(nactive >= para_max_run) if(waitpid(-pgid,NULL,0)>0) nactive--;
351  continue; // parent process goes next cycle
352  }
353 
354  // pid = 0, child process
355  In_GiNaC_Parallel = true;
356  GiNaC_Parallel_Level++;
357  try {
358  lst res_lst;
359  for(int ri=0; ri<nbatch; ri++) {
360  int i = bi*nbatch + ri;
361  if(i<ntotal) res_lst.append(f(i));
362  else break;
363  }
364  ostringstream garfn;
365  if(key == "") garfn << ppid << "/" << bi << ".gar";
366  else garfn << ppid << "/" << bi << "." << key << ".gar";
367  garWrite(garfn.str(), res_lst);
368  } catch(exception &p) { // NOT use Error
369  cout << ErrColor << "Failed in GiNaC_Parallel!" << RESET << endl;
370  cout << ErrColor << p.what() << RESET << endl;
371  }
372  exit(0);
373  }
374 
375  if(getpid()!=pgid) exit(0); // make sure child exit
376  while (waitpid(-pgid,NULL,0) != -1) { }
377  if(!nst && verb > 1) cout << endl;
378  if(getpid()!=ppid) exit(0); // make the forked group leader exit
379 
380  if(ec<2) break;
381  // check all *.gar and retry
382  bool all_gar_exists = true;
383  for(int bi=0; bi<btotal; bi++) {
384  ostringstream garfn;
385  if(key == "") garfn << ppid << "/" << bi << ".gar";
386  else garfn << ppid << "/" << bi << "." << key << ".gar";
387  if(!file_exists(garfn.str())) {
388  all_gar_exists = false;
389  break;
390  }
391  }
392  if(all_gar_exists) break;
393  para_max_run = para_max_run/2;
394  if(para_max_run<2) break;
395  fork_retried++;
396  }
397 
398  wait_label:
399  if(npid!=0) waitpid(npid,NULL,0); // wait nest process to exit
400 
401  //#define use_gar_write // use ReWR
402 
403  exvector ovec;
404  if(true) {
405  #ifdef use_gar_write
406  exvector ovec_tmp;
407  #endif
408  for(int bi=0; bi<btotal; bi++) {
409  if(verb > 1 && !nst) {
410  if(key == "") {
411  cout << "\r \r" << pre;
412  cout << "\\--Reading *.gar [" << (bi+1) << "/" << btotal << "] @ " << now(false) << flush;
413  } else {
414  cout << "\r \r" << pre;
415  cout << "\\--Reading *." << Color_HighLight << key << RESET << ".gar [" << (bi+1) << "/" << btotal << "] @ " << now(false) << flush;
416  }
417  }
418 
419  ostringstream garfn;
420  if(key == "") garfn << ppid << "/" << bi << ".gar";
421  else garfn << ppid << "/" << bi << "." << key << ".gar";
422  lst res_lst;
423  try {
424  if(file_exists(garfn.str())) {
425  res_lst = ex_to<lst>(garRead(garfn.str()));
426  remove(garfn.str().c_str());
427  goto done;
428  }
429  } catch(exception &p) { }
430  if(verb > 1 && !nst) cout << " - ReTry" << endl;
431  try {
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));
436  else break;
437  }
438  } catch(exception &p) {
439  cout << ErrColor << "Failed in GiNaC_Parallel!" << RESET << endl;
440  cout << ErrColor << p.what() << RESET << endl;
441  throw Error("GiNaC_Parallel_ReTRY: "+string(p.what()));
442  }
443  done: ;
444  #ifdef use_gar_write
445  if(GiNaC_Parallel_ReWR.find(key)==GiNaC_Parallel_ReWR.end() || GiNaC_Parallel_ReWR[key]) {
446  for(auto res : res_lst) ovec_tmp.push_back(res);
447  } else {
448  for(auto res : res_lst) ovec.push_back(res);
449  }
450  #else
451  for(auto res : res_lst) ovec.push_back(res);
452  #endif
453  }
454 
455  #ifdef use_gar_write
456  if(!nst && verb>1) {
457  cout << endl;
458  if(key == "") {
459  cout << "\r \r" << pre;
460  cout << "\\--ReWR" << flush;
461  } else {
462  cout << "\r \r" << pre;
463  cout << "\\--ReWR " << Color_HighLight << key << RESET << flush;
464  }
465  }
466 
467  if(GiNaC_Parallel_ReWR.find(key)==GiNaC_Parallel_ReWR.end() || GiNaC_Parallel_ReWR[key]) {
468  ostringstream garfn;
469  if(key == "") garfn << ppid << "/ReWR.gar";
470  else garfn << ppid << "/ReWR." << key << ".gar";
471  garWrite(garfn.str(), ovec_tmp);
472  ovec_tmp.clear();
473  garRead(garfn.str(), ovec);
474  }
475  #endif
476 
477  if(!nst && verb>1) cout << endl;
478  }
479 
480  if(rm) {
481  cmd.clear();
482  cmd.str("");
483  cmd << "rm -fr " << ppid;
484  rc = system(cmd.str().c_str());
485  }
486  if(ovec.size() != ntotal) {
487  throw Error("GiNaC_Parallel: The output size is wrong!");
488  }
489 
490  return exvector(std::move(ovec));
491  }
492 
499  std::vector<std::string> split(const std::string& s, char delimiter) {
500  std::vector<std::string> tokens;
501  std::string token;
502  std::istringstream tokenStream(s);
503  while (std::getline(tokenStream, token, delimiter)) {
504  tokens.push_back(token);
505  }
506  return tokens;
507  }
508 
514  ex lstHelper::sum(lst m) {
515  ex ret = 0;
516  for(int i=0; i<m.nops(); i++) ret += m.op(i);
517  return ret;
518  }
519 
525  string now(bool use_date) {
526  time_t timep;
527  time (&timep);
528  char tmp[64];
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) );
531  return tmp;
532  }
533 
539  lst gather_symbols(const ex & e) {
540  exset ss;
541  for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i) {
542  if(is_a<symbol>(*i)) ss.insert(*i);
543  }
544  lst sym_lst;
545  for(auto item : ss) sym_lst.append(item);
546  return sym_lst;
547  }
548 
554  lst gather_symbols(const exvector & ve) {
555  exset ss;
556  for(auto e : ve) {
557  for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i) {
558  if(is_a<symbol>(*i)) ss.insert(*i);
559  }
560  }
561  lst sym_lst;
562  for(auto item : ss) sym_lst.append(item);
563  return sym_lst;
564  }
565 
571  string RunOS(const string & cmd) {
572  char buf[128];
573  ostringstream oss;
574  FILE *fp = NULL;
575  if( (fp = popen(cmd.c_str(), "r")) != NULL) {
576  while(fgets(buf, 128, fp) != NULL) {
577  oss << buf;
578  }
579  pclose(fp);
580  fp = NULL;
581  }
582  return oss.str();
583  }
584 
585 
591  void garRead(const string &garfn, map<string, ex> &resMap) {
592  archive ar;
593  ifstream in(garfn);
594  in >> ar;
595  in.close();
596  for(int i=0; i<ar.num_expressions(); i++) {
597  string name;
598  ex res = ar.unarchive_ex(name, i);
599  resMap[name] = res;
600  }
601  }
602 
609  ex garRead(const string &garfn, const char* key) { // use the const char *, not string
610  archive ar;
611  ifstream in(garfn);
612  in >> ar;
613  in.close();
614  auto res = ar.unarchive_ex(key);
615  return res;
616  }
617 
623  ex garRead(const string &garfn) {
624  archive ar;
625  ifstream in(garfn);
626  in >> ar;
627  in.close();
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);
631  return res;
632  }
633 
639  void garWrite(const string &garfn, const map<string, ex> &resMap) {
640  archive ar;
641  for(const auto & item : resMap) {
642  ar.archive_ex(item.second, item.first.c_str());
643  }
644  ofstream out(garfn);
645  out << ar;
646  out.close();
647  }
648 
654  void garWrite(const string &garfn, const ex & res) {
655  archive ar;
656  ar.archive_ex(res, "res");
657  ar.archive_ex(19790923, "c");
658  ofstream out(garfn);
659  out << ar;
660  out.close();
661  }
662 
663 
670  ex str2ex(const string &expr, symtab stab) {
671  Parser par(stab);
672  ex ret = par.Read(expr);
673  return ret;
674  }
675 
681  ex str2ex(const string &expr) {
682  Parser par;
683  ex ret = par.Read(expr);
684  return ret;
685  }
686 
693  lst str2lst(const string &expr, symtab stab) {
694  Parser par(stab);
695  ex ret = par.Read(expr);
696  return ex_to<lst>(ret);
697  }
698 
704  lst str2lst(const string &expr) {
705  Parser par;
706  ex ret = par.Read(expr);
707  return ex_to<lst>(ret);
708  }
709 
715  string ex2str(const ex &expr) {
716  ostringstream oss;
717  oss << expr;
718  return oss.str();
719  }
720 
726  string ex2str(const exvector &expr) {
727  ostringstream oss;
728  oss << expr;
729  return oss.str();
730  }
731 
737  matrix lst2mat(const lst & ls) {
738  int nr = ls.nops();
739  int nc = ls.op(0).nops();
740  matrix mat(nr,nc);
741  for(int r=0; r<nr; r++) {
742  auto row = ls.op(r);
743  for(int c=0; c<nc; c++) mat(r,c) = row.op(c);
744  }
745  return mat;
746  }
747 
753  string ex2str(const exmap &expr) {
754  ostringstream oss;
755  oss << expr;
756  return oss.str();
757  }
758 
764  string ex2str(const exset &expr) {
765  ostringstream oss;
766  oss << expr;
767  return oss.str();
768  }
769 
775  string file2str(string filename) {
776  ifstream ifs(filename);
777  string ostr((istreambuf_iterator<char>(ifs)), (istreambuf_iterator<char>()));
778  ifs.close();
779  return ostr;
780  }
781 
787  void str2file(const string & ostr, string filename) {
788  std::ofstream ofs;
789  ofs.open(filename, ios::out);
790  ofs << ostr << endl;
791  ofs.close();
792  }
793 
794  // not finished, just for memo
795  void str2file(char * buff, FILE* fh) {
796  string nstr;
797  int n = nstr.length();
798  char nbuff[n+1];
799  strcpy(nbuff, nstr.c_str());
800  FILE * f = fmemopen(nbuff, n, "r");
801  }
802 
809  vector<string> file2strvec(string filename, bool skip_empty) {
810  ifstream fs(filename);
811  vector<string> ovec;
812  std::string line;
813  while (std::getline(fs, line)) {
814  if(!skip_empty || line.size()>0) ovec.push_back(line);
815  }
816  fs.close();
817  return ovec;
818  }
819 
825  ex file2ex(string filename) {
826  return str2ex(file2str(filename));
827  }
828 
835  ex file2ex(string filename, symtab st) {
836  return str2ex(file2str(filename), st);
837  }
838 
845  void ex2file(const ex & expr, string filename) {
846  std::ofstream ofs;
847  ofs.open(filename, ios::out);
848  ofs << expr << endl;
849  ofs.close();
850  }
851 
858  void ex2file(string filename, const ex & expr) {
859  ex2file(expr, filename);
860  }
861 
867  ex q2ex(__float128 num) {
868  char buffer[128];
869  quadmath_snprintf(buffer, sizeof buffer, "%.36QG", num);
870  numeric ret(buffer);
871  return ret;
872  }
873 
879  __float128 ex2q(ex num) {
880  ostringstream nss;
881  set_precision(40);
882  nss << num.evalf() << endl;
883  __float128 ret = strtoflt128(nss.str().c_str(), NULL);
884  reset_precision();
885  return ret;
886  }
887 
893  int ex2int(ex num) {
894  return ex_to<numeric>(num).to_int();
895  }
896 
902  lst vec2lst(const exvector & ev) {
903  return GiNaC::lst(ev.begin(), ev.end());
904  }
905 
911  exvector lst2vec(const lst & alst) {
912  exvector ret(alst.begin(), alst.end());
913  return exvector(std::move(ret));
914  }
915 
921  lst add2lst(const ex & expr) {
922  if(!is_a<add>(expr)) return lst{expr};
923  lst ret;
924  for(auto item : expr) ret.append(item);
925  return ret;
926  }
927 
933  lst mul2lst(const ex & expr) {
934  if(!is_a<mul>(expr)) return lst{expr};
935  lst ret;
936  for(auto item : expr) ret.append(item);
937  return ret;
938  }
939 
946  lst xlst(int bi, int ei) {
947  lst ret;
948  for(int i=bi; i<=ei; i++) ret.append(x(i));
949  return ret;
950  }
951 
957  lst xlst(int ei) {
958  return xlst(0, ei);
959  }
960 
968  ex series_ex(ex const & expr_in, ex const &s0, int sn0) {
969  static symbol ss;
970  return series_ex(expr_in.subs(s0==ss),ss,sn0).subs(ss==s0);
971  }
972 
980  ex series_ex(ex const & expr_in, symbol const &s0, int sn0) {
981  ex expr = expr_in;
982  if(!expr.has(s0)) return (sn0>=0 ? expr : 0);
983 
984  exset sset;
985  expr.find(pow(s0, w), sset);
986  numeric sn_lcm = 1;
987  exmap repl1st, repl2nd;
988  for(auto pi : sset) {
989  auto sn = expand(pi.op(1));
990  if(!is_a<numeric>(sn)) {
991  ex ex1 = 0, ex2 = 0;
992  for(auto item : add2lst(sn)) {
993  if(!is_a<numeric>(item)) ex1 += item;
994  else ex2 += item;
995  }
996  sn = ex2;
997  static symbol sss;
998  repl1st[pi] = sss * pow(s0, sn);
999  repl2nd[sss] = pow(s0, ex1);
1000  }
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));
1005  }
1006  sn_lcm = lcm(sn_lcm, ex_to<numeric>(sn).denom());
1007  }
1008  if(expr.has(sqrt(s0))) sn_lcm = lcm(sn_lcm, numeric(2));
1009  expr = expr.subs(repl1st);
1010 
1011  symbol s("s");
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));
1016 
1017  auto cvs = collect_lst(expr,s);
1018  ex ret = 0;
1019  for(auto cv : cvs) {
1020  bool ok = false;
1021  int exN = 1;
1022  while(exN<10) {
1023  expr = cv.op(1) + pow(s,sn+exN+2)+pow(s,sn+exN+3);
1024  expr = expr.series(s,sn+exN);
1025  ex ot = 0;
1026  for(int i=0; i<expr.nops(); i++) {
1027  if(is_order_function(expr.op(i))) {
1028  ot = expr.op(i);
1029  break;
1030  }
1031  }
1032  if(!is_order_function(ot)) {
1033  cerr << "expr = " << expr << endl;
1034  throw Error("series_ex: Not an Order term with " + ex2str(ot));
1035  }
1036  if(ot.op(0).degree(s)>sn) {
1037  expr = series_to_poly(expr);
1038  expr = collect_ex(expr,s);
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);
1041  }
1042  ok = true;
1043  break;
1044  }
1045  exN++;
1046  }
1047  if(!ok) throw Error("series_ex seems not working!");
1048  }
1049  ret = ret.subs(s==pow(s0,ex(1)/sn_lcm)); // need this for log-terms
1050  ret = ret.subs(repl2nd);
1051  ret = collect_ex(ret,s0);
1052  return ret;
1053  }
1054 
1063  ex diff_ex(ex const expr, ex const xp, unsigned nth, bool expand) {
1064  symbol s;
1065  ex res = expr.subs(xp==s);
1066  if(!expand) return res.diff(s,nth).subs(s==xp);
1067  auto cvs = collect_lst(res, s);
1068  res = 0;
1069  for(auto cv : cvs) {
1070  res += cv.op(0) * cv.op(1).diff(s,nth).subs(s==xp);
1071  }
1072  return res;
1073  }
1074 
1075  typedef vector<pair<ex,ex>> epvec_t; // first-pat, second-coeff
1076  typedef pair<ex,epvec_t> co_epvec_t; // first - overall coeff
1077  co_epvec_t power_expand_2(const co_epvec_t & co_epv, int n) {
1078  if(n==1) return co_epv;
1079  else if(n==2) {
1080  ex co = pow(co_epv.first,2);
1081  exmap pcmap;
1082  if(true) {
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;
1089  }
1090  }
1091  epvec_t epv;
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));
1096  }
1097  return make_pair(co, epvec_t(std::move(epv)));
1098  }
1099  int n2 = n/2;
1100  auto co_epv2 = power_expand_2(co_epv, n2);
1101  co_epv2 = power_expand_2(co_epv2, 2);
1102  if((n % 2)==0) return co_epvec_t(std::move(co_epv2));
1103 
1104  ex co = co_epv.first * co_epv2.first;
1105  exmap pcmap;
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;
1109  }
1110  for(auto ep2 : co_epv2.second) pcmap[ep2.first] += ep2.second * co_epv.first;
1111 
1112  epvec_t epv;
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));
1117  }
1118  return make_pair(co, epvec_t(std::move(epv)));
1119  }
1120 
1121  pair<ex,epvec_t> inner_expand_collect(ex const &expr_in, std::function<bool(const ex &)> has_func, int depth=0) {
1122  if(!has_func(expr_in)) {
1123  ex co;
1124  epvec_t epv;
1125  co = expr_in;
1126  return make_pair(co,epv);
1127  } if(is_a<add>(expr_in)) {
1128  ex co = 0;
1129  exmap pcmap;
1130  for(auto item : expr_in) {
1131  if(has_func(item)) {
1132  auto co_epv = inner_expand_collect(item, has_func, depth+1);
1133  co += co_epv.first;
1134  for(auto ep : co_epv.second) pcmap[ep.first] += ep.second;
1135  } else co += item;
1136  }
1137  epvec_t epv;
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));
1142  }
1143  return make_pair(co, epvec_t(std::move(epv)));
1144  } else if(is_a<mul>(expr_in)) {
1145  ex co = 1;
1146  epvec_t epv;
1147  for(auto item : expr_in) {
1148  if(!has_func(item)) {
1149  for(auto & ep : epv) ep.second *= item;
1150  co *= item;
1151  } else {
1152  exmap pcmap;
1153  auto co_epv = inner_expand_collect(item, has_func, depth+1);
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;
1157  }
1158  for(auto ep1 : epv) pcmap[ep1.first] += ep1.second * co_epv.first;
1159  co *= co_epv.first;
1160  epv.clear();
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));
1165  }
1166  }
1167  }
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();
1171  auto co_epv = inner_expand_collect(expr_in.op(0), has_func, depth+1);
1172  return power_expand_2(co_epv,n);
1173  } else {
1174  ex co = 0;
1175  epvec_t epv;
1176  epv.push_back(make_pair(expr_in, 1));
1177  return make_pair(co,epvec_t(std::move(epv)));
1178  }
1179  throw Error("inner_expand_collect unexpected region reached.");
1180  }
1181 
1188  ex expand_ex(ex const &expr_in, std::function<bool(const ex &)> has_func) {
1189  auto co_epv = inner_expand_collect(expr_in, has_func, 0);
1190  ex ret = co_epv.first;
1191  for(auto ep : co_epv.second) ret += ep.second * ep.first;
1192  return ret;
1193  }
1194 
1202  ex collect_ex(ex const &expr_in, std::function<bool(const ex &)> has_func, int opt) {
1203  auto cvs = collect_lst(expr_in, has_func, opt);
1204  exvector res_vec;
1205  res_vec.reserve(cvs.nops());
1206  for(auto cv : cvs) {
1207  auto cc = cv.op(0);
1208  auto vv = cv.op(1);
1209  if(cc.is_zero() || vv.is_zero()) continue;
1210  res_vec.push_back(cc * vv);
1211  }
1212  return add(res_vec);
1213  }
1214 
1222  lst collect_lst(ex const &expr_in, std::function<bool(const ex &)> has_func, int opt) {
1223  auto co_epv = inner_expand_collect(expr_in, has_func);
1224  ex cf = co_epv.first;
1225  lst res_lst;
1226  if(!is_zero(cf)) co_epv.second.push_back(make_pair(1,cf));
1227  for(auto ep : co_epv.second) {
1228  ex vv = ep.first;
1229  ex cc = ep.second;
1230  if(opt==o_normal || opt==o_fermat || opt==o_fermatfD || opt==o_fermatN || opt==o_flint || opt==o_flintf || opt==o_flintfD) cc = exnormal(cc,opt);
1231  else if(opt==o_form) cc = exfactor(cc,opt);
1232  else if(opt==o_normal_fermat) cc = exnormal(normal(cc),o_fermat);
1233  else if(opt==o_normal_form) cc = form_factor(normal(cc),o_fermat);
1234  else if(opt==o_fermat_form) cc = form_factor(fermat_normal(cc),o_fermat);
1235  if(!is_zero(cc)) res_lst.append(lst{cc, vv});
1236  }
1237  return res_lst;
1238  }
1239 
1245  ex EvalF(ex expr) {
1246  exset zs;
1247  //patterns needing evalf()
1248  expr.find(zeta(w), zs);
1249  expr.find(zeta(w,w), zs);
1250 
1251  lst repl;
1252  for(auto zi : zs) repl.append(zi==NN(zi));
1253  return expr.subs(repl);
1254  }
1255 
1256  ex EvalL(ex expr) {
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));
1260  }
1261  ex EvalQ(ex expr) {
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));
1265  }
1266  ex EvalMP(ex expr) {
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));
1270  }
1271 
1278  ex NN(ex expr, int digits) {
1279  set_precision(digits);
1280  auto nexpr = evalf(expr);
1281  reset_precision();
1282  return nexpr;
1283  }
1284 
1290  bool xPositive(ex const expr) {
1291  auto tmp = expand(expr);
1292  if(tmp.is_zero()) return false;
1293  bool ret = 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;
1298  }
1299  ret = true;
1300  } else {
1301  auto ntmp = NN(tmp.subs(x(w)==1).normal());
1302  ret = (is_a<numeric>(ntmp) && ntmp>0);
1303  }
1304  return ret;
1305  }
1306 
1312  int xSign(ex const expr) {
1313  if(is_zero(expr)) return 0;
1314  if(xPositive(expr)) return 1;
1315  else if(xPositive(ex(0)-expr)) return -1;
1316  else return 0;
1317  }
1318 
1324  void let_op_append(ex & ex_in, const ex item) {
1325  auto tmp = ex_to<lst>(ex_in);
1326  tmp.append(item);
1327  ex_in = tmp;
1328  }
1329 
1335  void let_op_prepend(ex & ex_in, const ex item) {
1336  auto tmp = ex_to<lst>(ex_in);
1337  tmp.prepend(item);
1338  ex_in = tmp;
1339  }
1340 
1345  void let_op_remove_last(ex & ex_in) {
1346  auto tmp = ex_to<lst>(ex_in);
1347  tmp.remove_last();
1348  ex_in = tmp;
1349  }
1350 
1355  void let_op_remove_first(ex & ex_in) {
1356  auto tmp = ex_to<lst>(ex_in);
1357  tmp.remove_first();
1358  ex_in = tmp;
1359  }
1360 
1367  void let_op_append(ex & ex_in, int index, ex const item) {
1368  auto tmp = ex_to<lst>(ex_in.op(index));
1369  tmp.append(item);
1370  ex_in.let_op(index) = tmp;
1371  }
1372 
1379  void let_op_append(lst & ex_in, int index, ex const item) {
1380  auto tmp = ex_to<lst>(ex_in.op(index));
1381  tmp.append(item);
1382  ex_in.let_op(index) = tmp;
1383  }
1384 
1392  void let_op_append(ex & ex_in, int index1, int index2, ex const item) {
1393  auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1394  tmp.append(item);
1395  ex_in.let_op(index1).let_op(index2) = tmp;
1396  }
1397 
1405  void let_op_append(lst & ex_in, int index1, int index2, ex const item) {
1406  auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1407  tmp.append(item);
1408  ex_in.let_op(index1).let_op(index2) = tmp;
1409  }
1410 
1417  void let_op_prepend(ex & ex_in, int index, ex const item) {
1418  auto tmp = ex_to<lst>(ex_in.op(index));
1419  tmp.prepend(item);
1420  ex_in.let_op(index) = tmp;
1421  }
1422 
1429  void let_op_prepend(lst & ex_in, int index, ex const item) {
1430  auto tmp = ex_to<lst>(ex_in.op(index));
1431  tmp.prepend(item);
1432  ex_in.let_op(index) = tmp;
1433  }
1434 
1442  void let_op_prepend(ex & ex_in, int index1, int index2, ex const item) {
1443  auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1444  tmp.prepend(item);
1445  ex_in.let_op(index1).let_op(index2) = tmp;
1446  }
1447 
1455  void let_op_prepend(lst & ex_in, int index1, int index2, ex const item) {
1456  auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1457  tmp.prepend(item);
1458  ex_in.let_op(index1).let_op(index2) = tmp;
1459  }
1460 
1466  void let_op_remove_last(ex & ex_in, int index) {
1467  auto tmp = ex_to<lst>(ex_in.op(index));
1468  tmp.remove_last();
1469  ex_in.let_op(index) = tmp;
1470  }
1471 
1477  void let_op_remove_last(lst & ex_in, int index) {
1478  auto tmp = ex_to<lst>(ex_in.op(index));
1479  tmp.remove_last();
1480  ex_in.let_op(index) = tmp;
1481  }
1482 
1489  void let_op_remove_last(ex & ex_in, int index1, int index2) {
1490  auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1491  tmp.remove_last();
1492  ex_in.let_op(index1).let_op(index2) = tmp;
1493  }
1494 
1501  void let_op_remove_last(lst & ex_in, int index1, int index2) {
1502  auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1503  tmp.remove_last();
1504  ex_in.let_op(index1).let_op(index2) = tmp;
1505  }
1506 
1512  void let_op_remove_first(ex & ex_in, int index) {
1513  auto tmp = ex_to<lst>(ex_in.op(index));
1514  tmp.remove_first();
1515  ex_in.let_op(index) = tmp;
1516  }
1517 
1523  void let_op_remove_first(lst & ex_in, int index) {
1524  auto tmp = ex_to<lst>(ex_in.op(index));
1525  tmp.remove_first();
1526  ex_in.let_op(index) = tmp;
1527  }
1528 
1535  void let_op_remove_first(ex & ex_in, int index1, int index2) {
1536  auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1537  tmp.remove_first();
1538  ex_in.let_op(index1).let_op(index2) = tmp;
1539  }
1540 
1547  void let_op_remove_first(lst & ex_in, int index1, int index2) {
1548  auto tmp = ex_to<lst>(ex_in.op(index1).op(index2));
1549  tmp.remove_first();
1550  ex_in.let_op(index1).let_op(index2) = tmp;
1551  }
1552 
1560  void let_op(ex &ex_in, int index1, int index2, const ex item) {
1561  ex_in.let_op(index1).let_op(index2) = item;
1562  }
1563 
1571  void let_op(lst &ex_in, int index1, int index2, const ex item) {
1572  ex_in.let_op(index1).let_op(index2) = item;
1573  }
1574 
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;
1585  }
1586 
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;
1597  }
1598 
1606  ex get_op(const ex ex_in, int index1, int index2) {
1607  return ex_in.op(index1).op(index2);
1608  }
1609 
1617  ex get_op(const lst ex_in, int index1, int index2) {
1618  return ex_in.op(index1).op(index2);
1619  }
1620 
1629  ex get_op(const ex ex_in, int index1, int index2, int index3) {
1630  return ex_in.op(index1).op(index2).op(index3);
1631  }
1632 
1641  ex get_op(const lst ex_in, int index1, int index2, int index3) {
1642  return ex_in.op(index1).op(index2).op(index3);
1643  }
1644 
1645  //-----------------------------------------------------------
1646  // XIntegral Class
1647  //-----------------------------------------------------------
1648  //GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(XIntegral, basic,print_func<print_dflt>(&XIntegral::print))
1649  GiNaC::registered_class_info & XIntegral::get_class_info_static() { return reg_info; }
1651  XIntegral * XIntegral::duplicate() const { XIntegral * bp = new XIntegral(*this); bp->setflag(GiNaC::status_flags::dynallocated); return bp; }
1652  void XIntegral::accept(GiNaC::visitor & v) const { if (visitor *p = dynamic_cast<visitor *>(&v)) p->visit(*this); else inherited::accept(v); }
1653  const GiNaC::registered_class_info &XIntegral::get_class_info() const { return get_class_info_static(); }
1654  GiNaC::registered_class_info &XIntegral::get_class_info() { return get_class_info_static(); }
1655  const char *XIntegral::class_name() const { return get_class_info_static().options.get_name(); }
1656  //GINAC_IMPLEMENT_REGISTERED_CLASS END
1657 
1661 
1662  int XIntegral::compare_same_type(const basic &other) const {
1663  if(!is_a<XIntegral>(other)) throw Error("XIntegral::compare_same_type");
1664  const XIntegral &o = static_cast<const XIntegral &>(other);
1665  auto c = Function.compare(o.Function);
1666  if(c!=0) return c;
1667  c = Exponent.compare(o.Exponent);
1668  if(c!=0) return c;
1669  c = Deltas.compare(o.Deltas);
1670  if(c!=0) return c;
1671  return 0;
1672  }
1673 
1674  void XIntegral::print(const print_dflt &c, unsigned level) const {
1675  c.s << "[∬" << "(" << Function << ")^(" << Exponent << ")";
1676  if(Deltas.nops()>0) c.s << "*𝜹(" << Deltas << ")";
1677  c.s << "]";
1678  }
1679 
1680  size_t XIntegral::nops() const { return 3; }
1681  ex XIntegral::op(size_t i) const {
1682  if(i==0) return Function;
1683  else if(i==1) return Exponent;
1684  else if(i==2) return Deltas;
1685  else throw Error("XIntegral::op, It is required that i<3.");
1686  }
1687  ex & XIntegral::let_op(size_t i) {
1688  ensure_if_modifiable();
1689  if(i==0) return Function;
1690  else if(i==1) return Exponent;
1691  else if(i==2) return Deltas;
1692  else throw Error("XIntegral::let_op, It is required that i<3.");
1693  }
1694 
1695  void XIntegral::archive(archive_node & n) const {
1696  inherited::archive(n);
1697  n.add_ex("Function", Function);
1698  n.add_ex("Exponent", Exponent);
1699  n.add_ex("Deltas", Deltas);
1700  }
1701 
1702  void XIntegral::read_archive(const archive_node& n) {
1703  inherited::read_archive(n);
1704  n.find_ex("Function", Function);
1705  n.find_ex("Exponent", Exponent);
1706  n.find_ex("Deltas", Deltas);
1707  }
1708 
1710  Function = fed.op(0);
1711  Exponent = fed.op(1);
1712  if(fed.nops()>2) Deltas = fed.op(2);
1713  }
1714 
1715  XIntegral::XIntegral(ex loops, ex ps, ex ns) {
1716  ex a = 0;
1717  ex pre = 1;
1718  ex rem = 0;
1719  lst delta;
1720  int xni = 0;
1721  lst funs, exps;
1722  for(int i=0; i<ps.nops(); i++) {
1723  if(is_zero(ns.op(i))) continue;
1724  auto cpi = Symbol::set_all(ps.op(i).expand());
1725  bool ltQ = false;
1726  for(auto li : loops) {
1727  if(cpi.has(li)) {
1728  ltQ = true;
1729  break;
1730  }
1731  }
1732 
1733  ex sgn = 0;
1734  if(!ltQ) {
1735  pre *= pow(Symbol::set_all(ps.op(i).expand()), ex(0)-ns.op(i));
1736  ns.let_op(i) = 0;
1737  ps.let_op(i) = 1;
1738  continue;
1739  } else if(is_a<numeric>(ns.op(i)) && (ns.op(i)<0)) {
1740  throw Error("XIntegral, negative powers not supported yet.");
1741  }
1742 
1743  // check loop^2
1744  for(auto li : loops) {
1745  if(!is_a<numeric>(cpi.coeff(li,2))) continue; // TODO: make sure positive
1746  numeric nm = ex_to<numeric>(cpi.coeff(li,2));
1747  if(nm.is_zero()) continue;
1748  sgn = nm>0 ? -1 : 1;
1749  break;
1750  }
1751  // check iEpsilon
1752  if(sgn.is_zero()) {
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;
1756  }
1757  // others
1758  if(sgn.is_zero()) {
1759  sgn = 1;
1760  if(is_a<numeric>(cpi) && ex_to<numeric>(cpi)>0) sgn = -1;
1761  else sgn = 1;
1762  }
1763 
1764  cpi = (ps.op(i)*sgn).subs(iEpsilon==0);
1765  if(sgn==-1) pre *= exp(I * Pi * ns.op(i));
1766  a += ns.op(i);
1767  auto xi = x(xni);
1768  rem += xi * cpi;
1769  delta.append(xi);
1770  if(!is_zero(ns.op(i)-1)) {
1771  funs.append(xi);
1772  exps.append(ns.op(i)-1);
1773  }
1774  pre /= tgamma(ns.op(i));
1775  xni++;
1776  }
1777  rem = Symbol::set_all(rem.expand());
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));
1780 
1781  // Loop
1782  ex u=1;
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);
1787  u *= (-t2);
1788  if(t2==0) {
1789  Function = lst{0};
1790  Exponent = lst{1};
1791  return;
1792  }
1793  rem = (t0 - pow(t1,2)/(4*t2)).expand();
1794  }
1795  rem = Symbol::set_all(rem).normal();
1796  u = Symbol::set_all(u).normal();
1797  rem = (rem * u).normal();
1798 
1799  funs.prepend(rem);
1800  exps.prepend(dl2-a);
1801  funs.prepend(u);
1802  exps.prepend(a-dl2-(4-2*ep)/2);
1803  auto num_den = numer_denom(pre);
1804  funs.append(num_den.op(0));
1805  exps.append(1);
1806  if(!is_zero(num_den.op(1))) {
1807  funs.append(num_den.op(1));
1808  exps.append(-1);
1809  }
1810 
1811  Function = funs;
1812  Exponent = exps;
1813  Deltas = lst{delta};
1814  }
1815 
1820  int CpuCores() {
1821  return omp_get_num_procs();
1822  }
1823 
1824 
1831  ex normal_fermat(const ex & expr, bool dfactor) {
1832  auto nd = fermat_numer_denom(expr, dfactor);
1833  return nd.op(0)/nd.op(1);
1834  }
1835 
1841  ex collect_factors(const ex & expr) {
1842  try {
1843  return collect_common_factors(expr);
1844  } catch(...) { }
1845  return expr;
1846  }
1847 
1854  ex exfactor(const ex & expr, int opt) {
1855  if(opt==o_none) return expr;
1856  else if(opt==o_form) return factor_form(expr);
1857  else if(opt==o_flint || opt==o_flintf) return factor_flint(expr);
1858  else return expr;
1859  }
1860 
1867  ex exexpand(const ex & expr, int opt) {
1868  if(opt==1) return expand(expr);
1869  else return expand(expr);
1870  }
1871 
1872  namespace {
1873  exvector ev_fermat_pern(const exvector & ev, int np) {
1874  if(ev.size()<np) return ev;
1875  int tn = ev.size();
1876  exvector evs;
1877  int ci = 0;
1878  while(ci<tn) {
1879  ex sum = 0;
1880  for(int i=0; i<np; i++) {
1881  sum += ev[ci];
1882  ci++;
1883  if(ci==tn) break;
1884  }
1885  ex tt = normal_fermat(sum);
1886  evs.push_back(tt);
1887  }
1888  return exvector(std::move(evs));
1889  }
1890 
1891  ex normal_fermat_pern(const ex & e, int np) {
1892  if(!is_a<add>(e)) return normal_fermat(e);
1893  exvector evs(e.begin(), e.end());
1894  while(evs.size()>np) evs = ev_fermat_pern(evs,np);
1895  ex res = normal_fermat(add(evs));
1896  return res;
1897  }
1898  }
1899 
1906  ex exnormal(const ex & expr, int opt) {
1907  if(opt<0) return normal_fermat_pern(expr,-opt);
1908  else if(opt==o_normal) return normal(expr);
1909  else if(opt==o_fermat) return normal_fermat(expr);
1910  else if(opt==o_flint || opt==o_flintf || opt==o_flintfD) return normal_flint(expr, opt);
1911  else if(opt==o_fermatfD) return normal_fermat(expr,true);
1912  else if(opt==o_fermatN) return numer_fermat(expr);
1913  return expr;
1914  }
1915 
1922  ex exnd(const ex & expr, int opt) {
1923  if(opt==1) return numer_denom(factor_form(expr));
1924  else if(opt==2) return numer_denom_fermat(expr);
1925  else return numer_denom(expr);
1926  }
1927 
1928  ex form_eval(const ex & expr) {
1929  static map<pid_t, Form> form_map;
1930  auto pid = getpid();
1931  if((form_map.find(pid)==form_map.end())) { // init section
1932  ostringstream ss;
1933  ss << "AutoDeclare Symbols fv;" << endl;
1934  form_map[pid].Init("form");
1935  form_map[pid].Execute(ss.str());
1936  }
1937  Form &fprc = form_map[pid];
1938 
1939  auto expr_in = expr;
1940  exmap map_rat;
1941  expr_in = expr_in.to_polynomial(map_rat);
1942 
1943  exset vs;
1944  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
1945  auto e = (*i);
1946  if(is_a<symbol>(e)) vs.insert(e);
1947  }
1948 
1949  exmap s2v;
1950  symtab st;
1951  int cid = 0;
1952  for(auto ss : vs) {
1953  cid++;
1954  string cvv = "fv"+to_string(cid);
1955  s2v[ss] = Symbol(cvv);
1956  st[cvv] = ss.subs(map_rat);
1957  }
1958  ostringstream oss;
1959  expr_in = expr_in.subs(s2v);
1960  oss << "Local ff = " << expr_in << ";" << endl;
1961  oss << ".sort" << endl;
1962  try {
1963  auto ostr = fprc.Execute(oss.str(), "ff");
1964  Parser fp(st);
1965  ex ret = fp.Read(ostr);
1966  return ret;
1967  } catch(Error& err) {
1968  cout << oss.str() << endl;
1969  form_map.erase(pid);
1970  throw;
1971  }
1972  }
1973 
1974  ex inner_factor_form(const ex & expr) {
1975  if(is_a<mul>(expr)) {
1976  ex res = 1;
1977  for(auto item : expr) res *= inner_factor_form(item);
1978  return res;
1979  }
1980  static map<pid_t, Form> form_map;
1981  auto pid = getpid();
1982  if((form_map.find(pid)==form_map.end())) { // init section
1983  ostringstream ss;
1984  ss << "AutoDeclare Symbols fv;" << endl;
1985  form_map[pid].Init("form");
1986  form_map[pid].Execute(ss.str());
1987  }
1988  Form &fprc = form_map[pid];
1989 
1990  auto expr_in = expr;
1991  exmap map_rat;
1992  expr_in = expr_in.to_polynomial(map_rat);
1993 
1994  exset vs;
1995  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
1996  auto e = (*i);
1997  if(is_a<symbol>(e)) vs.insert(e);
1998  }
1999 
2000  exmap s2v;
2001  symtab st;
2002  int cid = 0;
2003  for(auto ss : vs) {
2004  cid++;
2005  string cvv = "fv"+to_string(cid);
2006  s2v[ss] = Symbol(cvv);
2007  st[cvv] = ss.subs(map_rat);
2008  }
2009  ostringstream oss;
2010  expr_in = expr_in.subs(s2v);
2011  oss << "Local ff = " << expr_in << ";" << endl;
2012  oss << "Factorize ff;" << endl;
2013  oss << ".sort" << endl;
2014  try {
2015  auto ostr = fprc.Execute(oss.str(), "ff");
2016  string_replace_all(ostr, "[", "(");
2017  string_replace_all(ostr, "]", ")");
2018  string_replace_all(ostr, "\\\n", "");
2019  string_replace_all(ostr, " ","");
2020  Parser fp(st);
2021  ex ret = fp.Read(ostr);
2022  ex res = 1;
2023  Symbol sfac("factor_");
2024  for(auto item : add2lst(ret)) res *= item.subs(sfac==1);
2025  return res;
2026  } catch(Error& err) {
2027  cout << oss.str() << endl;
2028  form_map.erase(pid);
2029  throw;
2030  }
2031  }
2032 
2038  ex factor_form(const ex & expr, bool nd) {
2039  if(!nd) return inner_factor_form(expr);
2040  auto num_den = fermat_numer_denom(expr);
2041  if(is_zero(num_den.op(1)-1)) return inner_factor_form(num_den.op(0));
2042  return inner_factor_form(num_den.op(0))/inner_factor_form(num_den.op(1));
2043  }
2044 
2045  //-----------------------------------------------------------
2046  // MMAFormat Output
2047  //-----------------------------------------------------------
2048  MMAFormat::MMAFormat(ostream &os, unsigned opt) : print_dflt(os, opt) {}
2049  MMAFormat::MMAFormat() : print_dflt(std::cout) {}
2050  GINAC_IMPLEMENT_PRINT_CONTEXT(MMAFormat, print_dflt)
2051 
2052  const MMAFormat & MMAFormat::operator << (const basic & v) const {
2053  v.print(*this);
2054  return *this;
2055  }
2056  const MMAFormat & MMAFormat::operator << (const ex & v) const {
2057  v.print(*this);
2058  return *this;
2059  }
2060  const MMAFormat & MMAFormat::operator << (const lst & v) const {
2061  v.print(*this);
2062  return *this;
2063  }
2064  const MMAFormat & MMAFormat::operator<<(std::ostream& (*v)(std::ostream&)) const {
2065  s << v;
2066  return *this;
2067  }
2068 
2069  const MMAFormat & MMAFormat::operator << (const matrix & mat) const {
2070  s << "{";
2071  int nr = mat.rows();
2072  int nc = mat.cols();
2073  for(int r=0; r<nr; r++) {
2074  s << "{";
2075  for(int c=0; c<nc; c++) {
2076  mat(r,c).print(*this);
2077  if(c+1!=nc) s << ",";
2078  }
2079  s << "}";
2080  if(r+1!=nr) s << ",";
2081  }
2082  s << "}";
2083  return *this;
2084  }
2085 
2086  const MMAFormat & MMAFormat::operator << (const exvector & e) const {
2087  auto i = e.begin();
2088  auto vend = e.end();
2089  if (i==vend) { s << "{}"; return *this; }
2090  s << "{";
2091  while (true) {
2092  i->print(*this);
2093  ++i;
2094  if(i==vend) break;
2095  s << ",";
2096  }
2097  s << "}";
2098  return *this;
2099  }
2100 
2101  const MMAFormat & MMAFormat::operator << (const exset & e) const {
2102  auto i = e.begin();
2103  auto send = e.end();
2104  if (i==send) { s << "{}"; return *this; }
2105  s << "{";
2106  while (true) {
2107  i->print(*this);
2108  ++i;
2109  if(i==send) break;
2110  s << ",";
2111  }
2112  s << "}";
2113  return *this;
2114  }
2115 
2116  const MMAFormat & MMAFormat::operator << (const exmap & e) const {
2117  auto i = e.begin();
2118  auto mend = e.end();
2119  if (i==mend) { s << "{}"; return *this; }
2120  s << "{";
2121  while (true) {
2122  i->first.print(*this);
2123  s << "->";
2124  i->second.print(*this);
2125  ++i;
2126  if(i==mend) break;
2127  s << ",";
2128  }
2129  s << "}";
2130  return *this;
2131  }
2132 
2133  void garWrite(const exvector &exv, string garfn) {
2134  archive ar;
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());
2138  }
2139  ofstream out(garfn);
2140  out << ar;
2141  out.close();
2142  }
2143 
2144  void garRead(exvector &exv, string garfn) {
2145  archive ar;
2146  ifstream in(garfn);
2147  in >> ar;
2148  in.close();
2149 
2150  map<string, ex> dict;
2151  for(int i=0; i<ar.num_expressions(); i++) {
2152  string name;
2153  ex res = ar.unarchive_ex(name, i);
2154  dict[name] = res;
2155  }
2156 
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)];
2161  }
2162 
2163  ex add_collect_normal(const exvector &exv, ex const &pats, int opt) {
2164  auto cvs_vec = GiNaC_Parallel(exv.size(), [&exv,pats](int idx)->ex {
2165  return collect_lst(exv[idx], pats);
2166  }, "ColEx");
2167  exmap res_map;
2168  for(auto cvs : cvs_vec) for(auto cv : cvs) res_map[cv.op(1)] += cv.op(0);
2169  exvector res_vec;
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);
2173  }, "NorEx");
2174  return add(res_vec);
2175  }
2176 
2177  ex add_collect_normal(const exvector &exv, init_list const &pats, int opt) {
2178  auto cvs_vec = GiNaC_Parallel(exv.size(), [&exv,pats](int idx)->ex {
2179  return collect_lst(exv[idx], pats);
2180  }, "ColEx");
2181  exmap res_map;
2182  for(auto cvs : cvs_vec) for(auto cv : cvs) res_map[cv.op(1)] += cv.op(0);
2183  exvector res_vec;
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);
2187  }, "NorEx");
2188  return add(res_vec);
2189  }
2190 
2191  ex add_collect_normal(const exvector &exv, lst const &pats, int opt) {
2192  auto cvs_vec = GiNaC_Parallel(exv.size(), [&exv,pats](int idx)->ex {
2193  return collect_lst(exv[idx], pats);
2194  }, "ColEx");
2195  exmap res_map;
2196  for(auto cvs : cvs_vec) for(auto cv : cvs) res_map[cv.op(1)] += cv.op(0);
2197  exvector res_vec;
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);
2201  }, "NorEx");
2202  return add(res_vec);
2203  }
2204 
2205  ex add_collect_normal(const ex & e, ex const &pats, int opt) {
2206  if(!is_a<add>(e)) throw Error("add_collect_normal: input is NOT a add class.");
2207  exvector exv(e.begin(), e.end());
2208  return add_collect_normal(exv, pats, opt);
2209  }
2210 
2211  ex add_collect_normal(const ex & e, lst const &pats, int opt) {
2212  if(!is_a<add>(e)) throw Error("add_collect_normal: input is NOT a add class.");
2213  exvector exv(e.begin(), e.end());
2214  return add_collect_normal(exv, pats, opt);
2215  }
2216 
2217  ex add_collect_normal(const ex & e, init_list const &pats, int opt) {
2218  if(!is_a<add>(e)) throw Error("add_collect_normal: input is NOT a add class.");
2219  exvector exv(e.begin(), e.end());
2220  return add_collect_normal(exv, pats, opt);
2221  }
2222 
2223  bool has_w(const ex & e) {
2224  for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i) if(is_a<wildcard>(*i)) return true;
2225  return false;
2226  }
2227 
2228  void subs_w(exmap & repl) {
2229  auto r = repl;
2230  for(auto kv : r) if(is_a<mul>(kv.first) && !has_w(kv.first)) repl[w*kv.first]=w*kv.second;
2231  }
2232 
2233  void subs_w(lst & repl) {
2234  auto r = repl;
2235  for(auto item : r) if(is_a<mul>(item) && !has_w(item)) repl.append(w*item.op(0)==w*item.op(1));
2236  repl.sort();
2237  repl.unique();
2238  }
2239 
2240  void ReShare(const ex & e) {
2241  archive ar;
2242  ar.archive_ex(e, "e");
2243  }
2244 
2245  void ReShare(const lst & es) {
2246  archive ar;
2247  for(auto const & e : es) ar.archive_ex(e, "e");
2248  }
2249 
2250  void ReShare(const ex & e1, const ex & e2) {
2251  archive ar;
2252  ar.archive_ex(e1, "e");
2253  ar.archive_ex(e2, "e");
2254  }
2255 
2256  void ReShare(const ex & e1, const ex & e2, const ex & e3) {
2257  archive ar;
2258  ar.archive_ex(e1, "e");
2259  ar.archive_ex(e2, "e");
2260  ar.archive_ex(e3, "e");
2261  }
2262 
2263  void ReShare(const exvector & ev) {
2264  archive ar;
2265  for(auto & e : ev) ar.archive_ex(e, "e");
2266  }
2267 
2268  void ReShare(const exvector & ev1, const exvector & ev2) {
2269  archive ar;
2270  for(auto & e : ev1) ar.archive_ex(e, "e");
2271  for(auto & e : ev2) ar.archive_ex(e, "e");
2272  }
2273 
2274  ex nextprime(const ex & n) {
2275  auto v = cln::the<cln::cl_I>(ex_to<numeric>(n).to_cl_N());
2276  return numeric(cln::nextprobprime(v));
2277  }
2278  ex nextprime(int n) {
2279  auto v = cln::the<cln::cl_I>(numeric(n).to_cl_N());
2280  return numeric(cln::nextprobprime(v));
2281  }
2282 
2283  ex Rationalize(const ex & e, int dn) {
2284 
2285  static MapFunction R([](const ex & e, MapFunction & self)->ex{
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);
2294  });
2295 
2296  if(dn>0) set_precision(dn);
2297  ex res = R(e);
2298  if(dn>0) reset_precision();
2299  return res;
2300  }
2301 
2302  extern std::stack<cln::float_format_t> cln_prec_stack;
2303  extern std::stack<long> digits_stack;
2304  void set_precision(long prec, bool push) {
2305  if(push) {
2306  cln_prec_stack.push(cln::default_float_format);
2307  digits_stack.push(Digits);
2308  }
2309  Digits = prec;
2310  cln::default_float_format = cln::float_format(prec);
2311  }
2312 
2314  if(cln_prec_stack.empty()) return;
2315  auto digits = digits_stack.top();
2316  auto prec = cln_prec_stack.top();
2317  Digits = digits;
2318  cln::default_float_format = prec;
2319  cln_prec_stack.pop();
2320  digits_stack.pop();
2321  }
2322 
2323  long get_precision() {
2324  return cln::default_float_format;
2325  }
2326 
2327  void arg2map(int argc, char** argv, const char *optstring, std::map<char,std::string> & kv) {
2328  auto o_opterr = opterr;
2329  opterr = 0;
2330  while(true) {
2331  int opt = getopt(argc, argv, optstring);
2332  if(opt==-1) break;
2333  else if(opt!='?') {
2334  if(optarg) kv[opt] = optarg;
2335  else kv[opt] = "";
2336  }
2337  }
2338  argc -= optind;
2339  argv += optind;
2340  opterr = o_opterr;
2341  }
2342 
2343  bool has_symbol(const ex & e) {
2344  for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i) if(is_a<symbol>(*i)) return true;
2345  return false;
2346  }
2347 
2348 
2349 }
2350 
int * a
Definition: Functions.cpp:234
unsigned golden_ratio_hash(uintptr_t n)
Definition: BASIC.cpp:9
Basic header file.
#define IMPLEMENT_HAS(classname)
Definition: BASIC.h:24
#define MAGENTA
Definition: BASIC.h:85
#define WHITE
Definition: BASIC.h:87
#define RED
Definition: BASIC.h:81
#define DEFAULT_CTOR(classname)
Definition: BASIC.h:21
#define RESET
Definition: BASIC.h:79
#define IMPLEMENT_ALL(classname)
Definition: BASIC.h:30
class used to wrap error message
Definition: BASIC.h:242
string msg
Definition: BASIC.h:244
const char * what() const
Definition: BASIC.cpp:24
Error(string _msg)
Error constructor,.
Definition: BASIC.cpp:23
interface to communicate with Form program
Definition: BASIC.h:820
string Execute(const string &script, const string &out_var="[o]")
Definition: Process.cpp:271
class for Mathematica Format Output
Definition: BASIC.h:843
const MMAFormat & operator<<(const T &v) const
Definition: BASIC.h:848
MMAFormat(ostream &os, unsigned opt=0)
Definition: BASIC.cpp:2048
class to wrap map_function of GiNaC
Definition: BASIC.h:632
class to parse for string or file, helpful with Symbol class
Definition: BASIC.h:645
ex Read(const string &instr, bool s2S=true)
Definition: Functions.cpp:111
virtual ~visitor()
Definition: BASIC.cpp:60
class extended to GiNaC symbol class, represent a positive symbol
Definition: BASIC.h:113
void unset() const
Definition: BASIC.cpp:140
static void unset_all()
Definition: BASIC.cpp:146
ex series(const relational &s, int order, unsigned options=0) const override
Definition: BASIC.cpp:125
void archive(archive_node &n) const override
Symbol archive.
Definition: BASIC.cpp:88
ex imag_part() const override
Definition: BASIC.cpp:104
Symbol * duplicate() const override
Definition: BASIC.cpp:61
bool is_equal_same_type(const basic &other) const override
Definition: BASIC.cpp:113
ex conjugate() const override
Definition: BASIC.cpp:102
static ex set_all(const ex &expr)
Definition: BASIC.cpp:147
const char * class_name() const override
Definition: BASIC.cpp:65
void read_archive(const archive_node &n) override
Symbol read_archive.
Definition: BASIC.cpp:96
static GiNaC::registered_class_info & get_class_info_static()
Definition: BASIC.cpp:59
ex derivative(const symbol &s) const override
Definition: BASIC.cpp:119
unsigned calchash() const override
Definition: BASIC.cpp:105
ex evalf() const override
Definition: BASIC.cpp:101
void set(const ex &v) const
Definition: BASIC.cpp:139
static std::map< std::string, ex > Table
Definition: BASIC.h:160
ex eval() const override
Definition: BASIC.cpp:100
void set_name(string n)
Definition: BASIC.cpp:137
void accept(GiNaC::visitor &v) const override
Definition: BASIC.cpp:62
static exmap vmap
Definition: BASIC.h:167
ex real_part() const override
Definition: BASIC.cpp:103
int compare_same_type(const GiNaC::basic &other) const override
Definition: BASIC.cpp:75
const GiNaC::registered_class_info & get_class_info() const override
Definition: BASIC.cpp:63
XIntegral Class, preface to SecDec.
Definition: BASIC.h:740
ex op(size_t i) const override
Definition: BASIC.cpp:1681
int compare_same_type(const GiNaC::basic &other) const override
Definition: BASIC.cpp:1662
void accept(GiNaC::visitor &v) const override
Definition: BASIC.cpp:1652
XIntegral * duplicate() const override
Definition: BASIC.cpp:1651
void archive(archive_node &n) const override
Definition: BASIC.cpp:1695
ex & let_op(size_t i) override
Definition: BASIC.cpp:1687
void print(const print_dflt &c, unsigned level=0) const
Definition: BASIC.cpp:1674
void read_archive(const archive_node &n) override
Definition: BASIC.cpp:1702
const char * class_name() const override
Definition: BASIC.cpp:1655
static GiNaC::registered_class_info & get_class_info_static()
Definition: BASIC.cpp:1649
const GiNaC::registered_class_info & get_class_info() const override
Definition: BASIC.cpp:1653
size_t nops() const override
Definition: BASIC.cpp:1680
class extended to GiNaC symbol class, pure imaginary symbol
Definition: BASIC.h:183
bool is_equal_same_type(const basic &other) const override
Definition: BASIC.cpp:188
ex eval() const override
Definition: BASIC.cpp:236
ex evalf() const override
Definition: BASIC.cpp:237
static std::map< std::string, ex > Table
Definition: BASIC.h:230
void archive(archive_node &n) const override
Symbol archive.
Definition: BASIC.cpp:224
void set_name(string n)
Definition: BASIC.cpp:212
iSymbol * duplicate() const override
Definition: BASIC.cpp:166
ex real_part() const override
Definition: BASIC.cpp:239
ex derivative(const symbol &s) const override
Definition: BASIC.cpp:194
static GiNaC::registered_class_info & get_class_info_static()
Definition: BASIC.cpp:164
ex conjugate() const override
Definition: BASIC.cpp:238
const char * class_name() const override
Definition: BASIC.cpp:170
ex series(const relational &s, int order, unsigned options=0) const override
Definition: BASIC.cpp:200
void accept(GiNaC::visitor &v) const override
Definition: BASIC.cpp:167
ex imag_part() const override
Definition: BASIC.cpp:240
const GiNaC::registered_class_info & get_class_info() const override
Definition: BASIC.cpp:168
int compare_same_type(const GiNaC::basic &other) const override
Definition: BASIC.cpp:179
void read_archive(const archive_node &n) override
Symbol read_archive.
Definition: BASIC.cpp:232
unsigned calchash() const override
Definition: BASIC.cpp:214
static ex sum(lst m)
sum all elements in the list
Definition: BASIC.cpp:514
HepLib namespace.
Definition: BASIC.cpp:17
ex expand_ex(ex const &expr_in, std::function< bool(const ex &)> has_func)
the expand like Mathematica
Definition: BASIC.cpp:1188
const char * Color_HighLight
Definition: BASIC.cpp:248
ex collect_factors(const ex &expr)
a wrapper for collect_common_factors, catch errors
Definition: BASIC.cpp:1841
void arg2map(int argc, char **argv, const char *optstring, std::map< char, std::string > &kv)
Definition: BASIC.cpp:2327
__float128 ex2q(ex num)
ex of numeric to __float128
Definition: BASIC.cpp:879
const iSymbol iEpsilon
void reset_precision()
Definition: BASIC.cpp:2313
ex EvalF(ex expr)
the nuerical evaluation, Digits=100 will be used
Definition: BASIC.cpp:1245
lst str2lst(const string &expr, symtab stab)
convert string to the lst, using Parser internally
Definition: BASIC.cpp:693
ex exnd(const ex &expr, int opt)
num_den a expression
Definition: BASIC.cpp:1922
ex form_factor(const ex &expr, bool nd=true)
Definition: BASIC.h:472
lst gather_symbols(const ex &e)
get all symbols from input expression
Definition: BASIC.cpp:539
void let_op_prepend(ex &ex_in, const ex item)
preppend item into expression
Definition: BASIC.cpp:1335
ex EvalL(ex expr)
Definition: BASIC.cpp:1256
lst xlst(int bi, int ei)
return a lst: x(bi), x(bi+1), ..., x(ei)
Definition: BASIC.cpp:946
int GiNaC_Parallel_Process
Definition: Init.cpp:143
void let_op(ex &ex_in, int index1, int index2, const ex item)
update index1-th.index2-th of expression with item
Definition: BASIC.cpp:1560
std::initializer_list< ex > init_list
Definition: BASIC.h:47
ex q2ex(__float128 num)
__float128 to ex
Definition: BASIC.cpp:867
bool xPositive(ex const expr)
check the expr is xPositive, i.e., each x-monomial item is postive
Definition: BASIC.cpp:1290
const int o_form
Definition: Init.cpp:109
const int o_fermat
Definition: Init.cpp:106
void ex2file(const ex &expr, string filename)
export expression file
Definition: BASIC.cpp:845
const int o_none
Definition: Init.cpp:104
ex normal_fermat(const ex &expr, bool dfactor)
return the normalizatied expression, using fermat_numer_denom
Definition: BASIC.cpp:1831
ex NN(ex expr, int digits)
the nuerical evaluation
Definition: BASIC.cpp:1278
const int o_flint
Definition: Init.cpp:110
map< string, bool > GiNaC_Parallel_RM
Definition: Init.cpp:148
ex factor_flint(const ex &e, bool nd=true)
const char * ErrColor
Definition: BASIC.cpp:246
const Symbol ep
const int o_fermatfD
Definition: Init.cpp:107
ex collect_ex(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica
Definition: BASIC.cpp:1202
ex EvalQ(ex expr)
Definition: BASIC.cpp:1261
void garRead(const string &garfn, map< string, ex > &resMap)
garRead from file, and output in a map
Definition: BASIC.cpp:591
map< string, int > GiNaC_Parallel_NB
Definition: Init.cpp:147
void subs_w(exmap &repl)
Definition: BASIC.cpp:2228
const int o_fermat_form
Definition: Init.cpp:115
void ReShare(const ex &e)
Definition: BASIC.cpp:2240
bool file_exists(string fn)
Definition: BASIC.h:289
bool has_w(const ex &e)
Definition: BASIC.cpp:2223
ex factor_form(const ex &expr, bool nd)
factorize a expression using FORM
Definition: BASIC.cpp:2038
ex series_ex(ex const &expr_in, ex const &s0, int sn0)
the series like Mathematica, include s^n
Definition: BASIC.cpp:968
const char * WarnColor
Definition: BASIC.cpp:247
ex str2ex(const string &expr, symtab stab)
convert string to ex expression, using Parser internally
Definition: BASIC.cpp:670
ex form_eval(const ex &expr)
Definition: BASIC.cpp:1928
string now(bool use_date)
date/time string
Definition: BASIC.cpp:525
vector< pair< ex, ex > > epvec_t
Definition: BASIC.cpp:1075
ex numer_denom_fermat(const ex &expr, bool dfactor=false)
return the numerator and denominator after normalization
Definition: Fermat.cpp:23
string RunOS(const string &cmd)
run symtem command and return the output as string
Definition: BASIC.cpp:571
lst add2lst(const ex &expr)
convert add to lst
Definition: BASIC.cpp:921
ex w
Definition: Init.cpp:90
ex fermat_normal(const ex &expr, bool dfactor=false)
Definition: BASIC.h:465
const Symbol vs
int CpuCores()
return the cpu cores using OpenMP
Definition: BASIC.cpp:1820
map< string, int > GiNaC_Parallel_NP
Definition: Init.cpp:144
map< string, int > GiNaC_Parallel_Verb
Definition: Init.cpp:145
map< string, bool > GiNaC_Parallel_ReWR
Definition: Init.cpp:150
exvector lst2vec(const lst &alst)
convert lst to exvector
Definition: BASIC.cpp:911
const int o_normal
Definition: Init.cpp:105
ex get_op(const ex ex_in, int index1, int index2)
return index1-th.index2-th of expression
Definition: BASIC.cpp:1606
vector< string > file2strvec(string filename, bool skip_empty)
read file content to string vector
Definition: BASIC.cpp:809
ex file2ex(string filename)
read file content to ex
Definition: BASIC.cpp:825
ex exnormal(const ex &expr, int opt)
normalize a expression
Definition: BASIC.cpp:1906
const int o_flintf
Definition: Init.cpp:111
const int o_flintfD
Definition: Init.cpp:112
ex EvalMP(ex expr)
Definition: BASIC.cpp:1266
const int o_normal_fermat
Definition: Init.cpp:113
ex nextprime(const ex &n)
Definition: BASIC.cpp:2274
ex normal_flint(const ex &expr, int opt=o_flint)
ex inner_factor_form(const ex &expr)
Definition: BASIC.cpp:1974
bool In_GiNaC_Parallel
Definition: Init.cpp:142
lst mul2lst(const ex &expr)
convert mul to lst
Definition: BASIC.cpp:933
string file2str(string filename)
read file content to string
Definition: BASIC.cpp:775
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}, ... }
Definition: BASIC.cpp:1222
ex diff_ex(ex const expr, ex const xp, unsigned nth, bool expand)
the differential like Mathematica
Definition: BASIC.cpp:1063
int Verbose
Definition: Init.cpp:139
ex Rationalize(const ex &e, int dn)
Definition: BASIC.cpp:2283
void string_replace_all(string &str, const string &from, const string &to)
Definition: Functions.cpp:148
ex exexpand(const ex &expr, int opt)
factorize a expression
Definition: BASIC.cpp:1867
void garWrite(const string &garfn, const map< string, ex > &resMap)
garWrite to write the string-key map to the archive
Definition: BASIC.cpp:639
void let_op_append(ex &ex_in, const ex item)
append item into expression
Definition: BASIC.cpp:1324
lst vec2lst(const exvector &ev)
convert exvector to lst
Definition: BASIC.cpp:902
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition: BASIC.cpp:715
void let_op_remove_first(ex &ex_in)
remove first from expression
Definition: BASIC.cpp:1355
ex exfactor(const ex &expr, int opt)
factorize a expression
Definition: BASIC.cpp:1854
matrix lst2mat(const lst &ls)
convert 2Dim lst to matrix
Definition: BASIC.cpp:737
bool dir_exists(string dir)
Definition: BASIC.h:297
pair< ex, epvec_t > inner_expand_collect(ex const &expr_in, std::function< bool(const ex &)> has_func, int depth=0)
Definition: BASIC.cpp:1121
void str2file(const string &ostr, string filename)
export string to a file
Definition: BASIC.cpp:787
const int o_normal_form
Definition: Init.cpp:114
co_epvec_t power_expand_2(const co_epvec_t &co_epv, int n)
Definition: BASIC.cpp:1077
void set_precision(long prec, bool push)
Definition: BASIC.cpp:2304
ex add_collect_normal(const exvector &exv, ex const &pats, int opt)
Definition: BASIC.cpp:2163
map< string, string > GiNaC_Parallel_PRE
Definition: Init.cpp:149
int xSign(ex const expr)
the always sign for expr
Definition: BASIC.cpp:1312
pair< ex, epvec_t > co_epvec_t
Definition: BASIC.cpp:1076
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
Definition: BASIC.cpp:260
void let_op_remove_last(ex &ex_in)
remove last from expression
Definition: BASIC.cpp:1345
std::vector< std::string > split(const std::string &s, char delimiter)
split the string into serveral part, separated by the delimiter
Definition: BASIC.cpp:499
const int o_fermatN
Definition: Init.cpp:108
bool has_symbol(const ex &e)
Definition: BASIC.cpp:2343
std::stack< cln::float_format_t > cln_prec_stack
Definition: Init.cpp:350
std::stack< long > digits_stack
Definition: Init.cpp:351
long get_precision()
Definition: BASIC.cpp:2323
ex numer_fermat(const ex &expr)
Definition: Fermat.cpp:170
ex fermat_numer_denom(const ex &expr, bool dfactor=false)
Definition: BASIC.h:461
string PRE
Definition: Init.cpp:140
int ex2int(ex num)
ex to integer
Definition: BASIC.cpp:893
ex subs(const ex &e, init_list sl)
Definition: BASIC.h:51
int GiNaC_Parallel_Batch
Definition: Init.cpp:146