HepLib
Fermat.cpp
Go to the documentation of this file.
1 // Functions related to fermat program
2 
3 #include "BASIC.h"
4 
5 namespace HepLib {
6 
7  Fermat& Fermat::get() {
8  static map<pid_t, Fermat> fermat_map;
9  auto pid = getpid();
10  if((fermat_map.find(pid)==fermat_map.end())) { // init section
11  fermat_map[pid].Init();
12  fermat_map[pid].vmax = 0;
13  }
14  return fermat_map[pid];
15  }
16 
23  ex numer_denom_fermat(const ex & expr, bool dfactor) {
24  int _fermat_using_array = fermat_using_array;
25  bool use_ncheck = false;
26 
27  Fermat &fermat = Fermat::get();
28  int &v_max = fermat.vmax;
29 
30  auto expr_in = expr;
31  exmap map_rat;
32  expr_in = expr_in.to_rational(map_rat);
33 
34  lst rep_vs;
35  if(true) {
36  map<ex,long long,ex_is_less> s2c;
37  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
38  if(is_a<symbol>(*i)) s2c[*i]++;
39  }
40  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
41  for(auto kv : s2c) {
42  auto fw = fermat_weight.find(kv.first);
43  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
44  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
45  else sv1.push_back(lst{kv.second, kv.first});
46  }
47  sort_vec(sv1);
48  sort_vec(sv2);
49  sort_vec(sv3);
50  for(auto sv : sv1) rep_vs.append(sv.op(1));
51  for(auto sv : sv2) rep_vs.append(sv.op(1));
52  for(auto sv : sv3) rep_vs.append(sv.op(1));
53  }
54 
55  exmap v2f, f2v;
56  exmap nn_map;
57  auto nn_pi1 = nextprime(3);
58  auto nn_pi2 = nextprime(3);
59  int fvi = 0;
60  for(auto vi : rep_vs) {
61  auto name = "v" + to_string(fvi);
62  Symbol s(name);
63  v2f[vi] = s;
64  f2v[s] = vi;
65  fvi++;
66  nn_pi1 = nextprime(nn_pi2+1);
67  nn_pi2 = nextprime(nn_pi1+1);
68  nn_map[s] = nn_pi1/nn_pi2;
69  }
70 
71  stringstream ss;
72  if(fvi>111) {
73  cout << rep_vs << endl;
74  throw Error("Fermat: Too many variables.");
75  }
76  if(fvi>v_max) {
77  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
78  fermat.Execute(ss.str());
79  ss.clear();
80  ss.str("");
81  v_max = fvi;
82  }
83 
84  ex nn_chk=0, num, den;
85  lst item;
86  if(!is_a<add>(expr_in)) item = lst{ expr_in };
87  else for(auto ii : expr_in) item.append(ii);
88  //sort_lst(item); // no need
89  if(item.nops()>999999) _fermat_using_array = 0;
90  if(_fermat_using_array) ss << "Array m[" << item.nops() << "];" << endl;
91  else ss << "res:=0;" << endl;
92  fermat.Execute(ss.str());
93  ss.clear();
94  ss.str("");
95 
96  for(int i=0; i<item.nops(); i++) {
97  ex tt = item.op(i).subs(v2f, subs_options::no_pattern);
98  if(use_ncheck) nn_chk += tt.subs(nn_map, subs_options::no_pattern);
99  if(_fermat_using_array) ss << "m[" << (i+1) << "]:=";
100  else ss << "item:=";
101  ss << tt << ";" << endl;
102  if(!_fermat_using_array) ss << "res:=*res+*item;" << endl;
103  fermat.Execute(ss.str());
104  ss.clear();
105  ss.str("");
106  }
107  if(_fermat_using_array) {
108  if(_fermat_using_array==1) ss << "res:=Sumup([m]);" << endl;
109  else if(_fermat_using_array==2) ss << "res:=Sigma<i=1,"<<item.nops()<<">(*m[i]);" << endl;
110  else {
111  ss << "n:=" << item.nops() << ";" << endl;
112  ss << "while n>1 do n2:=n\\2; for i=1,n2 do m[i]:=*m[i]+*m[n+1-i] od; &_G; if (n|2)=0 then n:=n2 else n:=n2+1 fi od;" << endl;
113  ss << "res:=*m[1];" << endl;
114  }
115  fermat.Execute(ss.str());
116  ss.clear();
117  ss.str("");
118  }
119 
120  static string bstr("[-begin-]"), estr("[-end-]");
121  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
122  ss << "!('" <<bstr<< "','{',Numer(^res),',',Denom(^res),'}','" <<estr<< "')" << endl;
123  auto ostr = fermat.Execute(ss.str());
124  ss.clear();
125  ss.str("");
126 
127  // note the order,(normal_fermat will be called again in factor_form)
128  ss << "&(U=0);" << endl; // disable ugly printing
129  if(_fermat_using_array) ss << "@(res,[m]);" << endl;
130  else ss << "@(res,item);" << endl;
131  ss << "&_G;" << endl;
132  fermat.Execute(ss.str());
133  ss.clear();
134  ss.str("");
135 
136  // make sure last char is 0
137  if(ostr[ostr.length()-1]!='0') throw Error("fermat_together: last char is NOT 0.");
138  ostr = ostr.substr(0, ostr.length()-1);
139  auto cpos = ostr.find(bstr);
140  if(cpos==string::npos) throw Error(bstr+" NOT Found.");
141  ostr = ostr.substr(cpos+bstr.length(),string::npos);
142  cpos = ostr.find(estr);
143  if(cpos==string::npos) throw Error(estr+" NOT Found.");
144  ostr = ostr.substr(0,cpos);
145  string_trim(ostr);
146 
147  symtab st;
148  Parser fp(st);
149  auto ret = fp.Read(ostr);
150  //ReShare(ret,expr);
151  num = ret.op(0);
152  if(dfactor) den = factor_form(ret.op(1));
153  else den = ret.op(1);
154  //fermat.Exit();
155 
156  if(use_ncheck) {
157  auto nn_ret = subs(num/den,nn_map);
158  if(nn_chk-nn_ret!=0) {
159  cout << nn_chk << " : " << nn_ret << endl;
160  throw Error("fermat_together: N Check Failed.");
161  }
162  }
163 
164  num = num.subs(f2v,subs_options::no_pattern).subs(map_rat,nopat);
165  den = den.subs(f2v,subs_options::no_pattern).subs(map_rat,nopat);
166  return lst{num, den};
167 
168  }
169 
170  ex numer_fermat(const ex & expr) {
171  int _fermat_using_array = fermat_using_array;
172  bool use_ncheck = false;
173 
174  Fermat &fermat = Fermat::get();
175  int &v_max = fermat.vmax;
176 
177  auto expr_in = expr;
178  exmap map_rat;
179  expr_in = expr_in.to_polynomial(map_rat);
180 
181  lst rep_vs;
182  if(true) {
183  map<ex,long long,ex_is_less> s2c;
184  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
185  if(is_a<symbol>(*i)) s2c[*i]++;
186  }
187  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
188  for(auto kv : s2c) {
189  auto fw = fermat_weight.find(kv.first);
190  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
191  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
192  else sv1.push_back(lst{kv.second, kv.first});
193  }
194  sort_vec(sv1);
195  sort_vec(sv2);
196  sort_vec(sv3);
197  for(auto sv : sv1) rep_vs.append(sv.op(1));
198  for(auto sv : sv2) rep_vs.append(sv.op(1));
199  for(auto sv : sv3) rep_vs.append(sv.op(1));
200  }
201 
202  exmap v2f, f2v;
203  exmap nn_map;
204  auto nn_pi1 = nextprime(3);
205  auto nn_pi2 = nextprime(3);
206  int fvi = 0;
207  for(auto vi : rep_vs) {
208  auto name = "v" + to_string(fvi);
209  Symbol s(name);
210  v2f[vi] = s;
211  f2v[s] = vi;
212  fvi++;
213  nn_pi1 = nextprime(nn_pi2+1);
214  nn_pi2 = nextprime(nn_pi1+1);
215  nn_map[s] = nn_pi1/nn_pi2;
216  }
217 
218  stringstream ss;
219  if(fvi>111) {
220  cout << rep_vs << endl;
221  throw Error("Fermat: Too many variables.");
222  }
223  if(fvi>v_max) {
224  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
225  fermat.Execute(ss.str());
226  ss.clear();
227  ss.str("");
228  v_max = fvi;
229  }
230 
231  ex nn_chk=0, nres;
232  lst item;
233  if(!is_a<add>(expr_in)) item = lst{ expr_in };
234  else for(auto ii : expr_in) item.append(ii);
235  //sort_lst(item); // no need
236  if(item.nops()>999999) _fermat_using_array = 0;
237  if(_fermat_using_array) ss << "Array m[" << item.nops() << "];" << endl;
238  else ss << "res:=0;" << endl;
239  fermat.Execute(ss.str());
240  ss.clear();
241  ss.str("");
242 
243  for(int i=0; i<item.nops(); i++) {
244  ex tt = item.op(i).subs(v2f, subs_options::no_pattern);
245  if(use_ncheck) nn_chk += tt.subs(nn_map, subs_options::no_pattern);
246  if(_fermat_using_array) ss << "m[" << (i+1) << "]:=";
247  else ss << "item:=";
248  ss << tt << ";" << endl;
249  if(!_fermat_using_array) ss << "res:=*res+*item;" << endl;
250  fermat.Execute(ss.str());
251  ss.clear();
252  ss.str("");
253  }
254  if(_fermat_using_array) {
255  if(_fermat_using_array==1) ss << "res:=Sumup([m]);" << endl;
256  else if(_fermat_using_array==2) ss << "res:=Sigma<i=1,"<<item.nops()<<">(*m[i]);" << endl;
257  else {
258  ss << "n:=" << item.nops() << ";" << endl;
259  ss << "while n>1 do n2:=n\\2; for i=1,n2 do m[i]:=*m[i]+*m[n+1-i] od; &_G; if (n|2)=0 then n:=n2 else n:=n2+1 fi od;" << endl;
260  ss << "res:=*m[1];" << endl;
261  }
262  fermat.Execute(ss.str());
263  ss.clear();
264  ss.str("");
265  }
266 
267  static string bstr("[-begin-]"), estr("[-end-]");
268  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
269  ss << "!('" <<bstr<< "',res,'" <<estr<< "')" << endl;
270  auto ostr = fermat.Execute(ss.str());
271  ss.clear();
272  ss.str("");
273 
274  // note the order,(normal_fermat will be called again in factor_form)
275  ss << "&(U=0);" << endl; // disable ugly printing
276  if(_fermat_using_array) ss << "@(res,[m]);" << endl;
277  else ss << "@(res,item);" << endl;
278  ss << "&_G;" << endl;
279  fermat.Execute(ss.str());
280  ss.clear();
281  ss.str("");
282 
283  // make sure last char is 0
284  if(ostr[ostr.length()-1]!='0') throw Error("fermat_together: last char is NOT 0.");
285  ostr = ostr.substr(0, ostr.length()-1);
286  auto cpos = ostr.find(bstr);
287  if(cpos==string::npos) throw Error(bstr+" NOT Found.");
288  ostr = ostr.substr(cpos+bstr.length(),string::npos);
289  cpos = ostr.find(estr);
290  if(cpos==string::npos) throw Error(estr+" NOT Found.");
291  ostr = ostr.substr(0,cpos);
292  string_trim(ostr);
293 
294  symtab st;
295  Parser fp(st);
296  auto ret = fp.Read(ostr);
297  //fermat.Exit();
298 
299  if(use_ncheck) {
300  auto nn_ret = subs(ret,nn_map);
301  if(nn_chk-nn_ret!=0) {
302  cout << nn_chk << " : " << nn_ret << endl;
303  throw Error("fermat_together: N Check Failed.");
304  }
305  }
306 
307  ret = ret.subs(f2v,subs_options::no_pattern).subs(map_rat,subs_options::no_pattern);
308  return ret;
309 
310  }
311 
317  ex fermat_eval(const ex & expr) {
318  Fermat &fermat = Fermat::get();
319  int &v_max = fermat.vmax;
320 
321  auto expr_in = expr;
322  exmap map_rat;
323  expr_in = expr_in.to_rational(map_rat);
324 
325  lst rep_vs;
326  if(true) {
327  map<ex,long long,ex_is_less> s2c;
328  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
329  if(is_a<symbol>(*i)) s2c[*i]++;
330  }
331  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
332  for(auto kv : s2c) {
333  auto fw = fermat_weight.find(kv.first);
334  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
335  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
336  else sv1.push_back(lst{kv.second, kv.first});
337  }
338  sort_vec(sv1);
339  sort_vec(sv2);
340  sort_vec(sv3);
341  for(auto sv : sv1) rep_vs.append(sv.op(1));
342  for(auto sv : sv2) rep_vs.append(sv.op(1));
343  for(auto sv : sv3) rep_vs.append(sv.op(1));
344  }
345 
346  exmap v2f, f2v;
347  int fvi = 0;
348  for(auto vi : rep_vs) {
349  auto name = "v" + to_string(fvi);
350  Symbol s(name);
351  v2f[vi] = s;
352  f2v[s] = vi;
353  fvi++;
354  }
355  expr_in = expr_in.subs(v2f);
356 
357  stringstream ss;
358  if(fvi>111) {
359  cout << rep_vs << endl;
360  throw Error("Fermat: Too many variables.");
361  }
362  if(fvi>v_max) {
363  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
364  fermat.Execute(ss.str());
365  ss.clear();
366  ss.str("");
367  v_max = fvi;
368  }
369 
370  ex res;
371  ss << "res:=" << expr_in << ";" << endl;
372  fermat.Execute(ss.str());
373  ss.clear();
374  ss.str("");
375 
376  static string bstr("[-begin-]"), estr("[-end-]");
377  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
378  ss << "!('" <<bstr<< "',res,'" <<estr<< "')" << endl;
379  auto ostr = fermat.Execute(ss.str());
380  ss.clear();
381  ss.str("");
382 
383  // note the order,(normal_fermat will be called again in factor_form)
384  ss << "&(U=0);" << endl; // disable ugly printing
385  ss << "@(res);" << endl;
386  ss << "&_G;" << endl;
387  fermat.Execute(ss.str());
388  ss.clear();
389  ss.str("");
390 
391  // make sure last char is 0
392  if(ostr[ostr.length()-1]!='0') throw Error("fermat_together: last char is NOT 0.");
393  ostr = ostr.substr(0, ostr.length()-1);
394  auto cpos = ostr.find(bstr);
395  if(cpos==string::npos) throw Error(bstr+" NOT Found.");
396  ostr = ostr.substr(cpos+bstr.length(),string::npos);
397  cpos = ostr.find(estr);
398  if(cpos==string::npos) throw Error(estr+" NOT Found.");
399  ostr = ostr.substr(0,cpos);
400  string_trim(ostr);
401 
402  symtab st;
403  Parser fp(st);
404  res = fp.Read(ostr);
405  res = res.subs(f2v,nopat);
406  return res.subs(map_rat,nopat);
407  }
408 
409  matrix fermat_Redrowech(const matrix & mat_in) {
410  Fermat &fermat = Fermat::get();
411  int &v_max = fermat.vmax;
412 
413  exmap map_rat;
414  int nrow = mat_in.rows();
415  int ncol = mat_in.cols();
416  matrix mat(nrow, ncol);
417  for(int r=0; r<nrow; r++) for(int c=0; c<ncol; c++) mat(r,c) = mat_in(r,c).to_rational(map_rat);
418 
419  lst rep_vs;
420  if(true) {
421  map<ex,long long,ex_is_less> s2c;
422  ex expr_in = mat;
423  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
424  if(is_a<symbol>(*i)) s2c[*i]++;
425  }
426  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
427  for(auto kv : s2c) {
428  auto fw = fermat_weight.find(kv.first);
429  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
430  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
431  else sv1.push_back(lst{kv.second, kv.first});
432  }
433  sort_vec(sv1);
434  sort_vec(sv2);
435  sort_vec(sv3);
436  for(auto sv : sv1) rep_vs.append(sv.op(1));
437  for(auto sv : sv2) rep_vs.append(sv.op(1));
438  for(auto sv : sv3) rep_vs.append(sv.op(1));
439  }
440 
441  exmap v2f;
442  symtab st;
443  int fvi = 0;
444  for(auto vi : rep_vs) {
445  auto name = "v" + to_string(fvi);
446  v2f[vi] = Symbol(name);
447  st[name] = vi;
448  fvi++;
449  }
450 
451  stringstream ss;
452  if(fvi>111) {
453  cout << rep_vs << endl;
454  throw Error("Fermat: Too many variables.");
455  }
456  if(fvi>v_max) {
457  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
458  fermat.Execute(ss.str());
459  ss.clear();
460  ss.str("");
461  v_max = fvi;
462  }
463 
464  ss << "Array m[" << nrow << "," << ncol << "];" << endl;
465  fermat.Execute(ss.str());
466  ss.clear();
467  ss.str("");
468 
469  ss << "[m]:=[(";
470  for(int c=0; c<ncol; c++) {
471  for(int r=0; r<nrow; r++) {
472  ss << mat(r,c).subs(iEpsilon==0,nopat).subs(v2f,nopat) << ",";
473  }
474  }
475  ss << ")];" << endl;
476  ss << "Redrowech([m]);" << endl;
477  auto tmp = ss.str();
478  string_replace_all(tmp,",)]",")]");
479  fermat.Execute(tmp);
480  ss.clear();
481  ss.str("");
482 
483  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
484  ss << "![m" << endl;
485  auto ostr = fermat.Execute(ss.str());
486  ss.clear();
487  ss.str("");
488  //fermat.Exit();
489 
490  // note the order, before exfactor (normal_fermat will be called again here)
491  ss << "&(U=0);" << endl; // disable ugly printing
492  ss << "@([m]);" << endl;
493  ss << "&_G;" << endl;
494  fermat.Execute(ss.str());
495  ss.clear();
496  ss.str("");
497 
498  // make sure last char is 0
499  if(ostr[ostr.length()-1]!='0') throw Error("Direc::Export, last char is NOT 0.");
500  ostr = ostr.substr(0, ostr.length()-1);
501  string_trim(ostr);
502 
503  ostr.erase(0, ostr.find(":=")+2);
504  string_replace_all(ostr, "[", "{");
505  string_replace_all(ostr, "]", "}");
506  Parser fp(st);
507  matrix mr(nrow, ncol);
508  auto res = fp.Read(ostr);
509  for(int r=0; r<nrow; r++) {
510  auto cur = res.op(r);
511  for(int c=0; c<ncol; c++) mr(r,c) = cur.op(c).subs(map_rat);
512  }
513  return mr;
514  }
515 
516  matrix fermat_Redrowech_Sparse(const matrix & mat_in) {
517  Fermat &fermat = Fermat::get();
518  int &v_max = fermat.vmax;
519 
520  exmap map_rat;
521  int nrow = mat_in.rows();
522  int ncol = mat_in.cols();
523  matrix mat(nrow, ncol);
524  matrix ret_mat(nrow, ncol);
525  for(int r=0; r<nrow; r++) for(int c=0; c<ncol; c++) mat(r,c) = mat_in(r,c).to_rational(map_rat);
526 
527  lst rep_vs;
528  if(true) {
529  map<ex,long long,ex_is_less> s2c;
530  ex expr_in = mat;
531  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
532  if(is_a<symbol>(*i)) s2c[*i]++;
533  }
534  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
535  for(auto kv : s2c) {
536  auto fw = fermat_weight.find(kv.first);
537  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
538  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
539  else sv1.push_back(lst{kv.second, kv.first});
540  }
541  sort_vec(sv1);
542  sort_vec(sv2);
543  sort_vec(sv3);
544  for(auto sv : sv1) rep_vs.append(sv.op(1));
545  for(auto sv : sv2) rep_vs.append(sv.op(1));
546  for(auto sv : sv3) rep_vs.append(sv.op(1));
547  }
548 
549  exmap v2f;
550  symtab st;
551  int fvi = 0;
552  for(auto vi : rep_vs) {
553  auto name = "v" + to_string(fvi);
554  v2f[vi] = Symbol(name);
555  st[name] = vi;
556  fvi++;
557  }
558 
559  stringstream ss;
560  if(fvi>111) {
561  cout << rep_vs << endl;
562  throw Error("Fermat: Too many variables.");
563  }
564  if(fvi>v_max) {
565  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
566  fermat.Execute(ss.str());
567  ss.clear();
568  ss.str("");
569  v_max = fvi;
570  }
571 
572  ostringstream oss;
573  oss << "[";
574  bool r1 = true;
575  for(int r=0; r<nrow; r++) {
576  bool c1 = true;
577  for(int c=0; c<ncol; c++) {
578  if(mat(r,c).is_zero()) continue;
579  if(c1) {
580  if(!r1) oss << "`" << endl;
581  oss << "[" << r+1 << ",";
582  r1 = c1 = false;
583  } else oss << ",";
584  oss << "[" << c+1 << "," << mat(r,c).subs(v2f) << "]";
585  }
586  if(!c1) oss << "]";
587  }
588  oss << "];" << endl;
589 
590  ss << "Array m[" << nrow << "," << ncol+1 << "] Sparse;" << endl;
591  fermat.Execute(ss.str());
592  ss.clear();
593  ss.str("");
594 
595  ss << "[m]:=" << oss.str();
596  fermat.Execute(ss.str());
597  ss.clear();
598  ss.str("");
599 
600  fermat.Execute("Redrowech([m]);");
601  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
602 
603  if(true) { // read matrix m
604  ss << "![m]" << endl;
605  auto ostr = fermat.Execute(ss.str());
606  ss.clear();
607  ss.str("");
608 
609  // make sure last char is 0
610  if(ostr[ostr.length()-1]!='0') throw Error("RowReduce, last char is NOT 0.");
611  ostr = ostr.substr(0, ostr.length()-1);
612  string_trim(ostr);
613 
614  ostr.erase(0, ostr.find(":=")+2);
615  size_t sn = ostr.length();
616  char lc;
617  for(size_t i=0; i<sn; i++) {
618  char & c = ostr[i];
619  if(c=='[') {
620  c = '{';
621  if(i>0 && lc=='}') ostr[i-1] = ',';
622  } else if(c==']') c = '}';
623  else if(c==' '||c=='\t'||c=='\n'||c=='\r') continue;
624  lc = c;
625  }
626 
627  Parser fp(st);
628  auto res = fp.Read(ostr);
629  for(auto const & item : res) {
630  int r = -1;
631  for(auto const & it : item) {
632  if(r==-1) r = ex2int(it)-1;
633  else {
634  int c = ex2int(it.op(0))-1;
635  ret_mat(r,c) = it.op(1).subs(map_rat);
636  }
637  }
638  }
639  }
640  // note the order, before exfactor (normal_fermat will be called again here)
641  ss << "&(U=0);" << endl; // disable ugly printing
642  ss << "@([m]);" << endl;
643  ss << "&_G;" << endl;
644  fermat.Execute(ss.str());
645  ss.clear();
646  ss.str("");
647  return ret_mat;
648  }
649 
650  ex fermat_Det(const matrix & mat_in) {
651  Fermat &fermat = Fermat::get();
652  int &v_max = fermat.vmax;
653 
654  exmap map_rat;
655  int nrow = mat_in.rows();
656  int ncol = mat_in.cols();
657  matrix mat(nrow, ncol);
658  for(int r=0; r<nrow; r++) for(int c=0; c<ncol; c++) mat(r,c) = mat_in(r,c).to_rational(map_rat);
659 
660  lst rep_vs;
661  if(true) {
662  map<ex,long long,ex_is_less> s2c;
663  ex expr_in = mat;
664  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
665  if(is_a<symbol>(*i)) s2c[*i]++;
666  }
667  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
668  for(auto kv : s2c) {
669  auto fw = fermat_weight.find(kv.first);
670  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
671  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
672  else sv1.push_back(lst{kv.second, kv.first});
673  }
674  sort_vec(sv1);
675  sort_vec(sv2);
676  sort_vec(sv3);
677  for(auto sv : sv1) rep_vs.append(sv.op(1));
678  for(auto sv : sv2) rep_vs.append(sv.op(1));
679  for(auto sv : sv3) rep_vs.append(sv.op(1));
680  }
681 
682  exmap v2f;
683  symtab st;
684  int fvi = 0;
685  for(auto vi : rep_vs) {
686  auto name = "v" + to_string(fvi);
687  v2f[vi] = Symbol(name);
688  st[name] = vi;
689  fvi++;
690  }
691 
692  stringstream ss;
693  if(fvi>111) {
694  cout << rep_vs << endl;
695  throw Error("Fermat: Too many variables.");
696  }
697  if(fvi>v_max) {
698  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
699  fermat.Execute(ss.str());
700  ss.clear();
701  ss.str("");
702  v_max = fvi;
703  }
704 
705  ss << "Array m[" << nrow << "," << ncol << "];" << endl;
706  fermat.Execute(ss.str());
707  ss.clear();
708  ss.str("");
709 
710  ss << "[m]:=[(";
711  for(int c=0; c<ncol; c++) {
712  for(int r=0; r<nrow; r++) {
713  ss << mat(r,c).subs(iEpsilon==0,nopat).subs(v2f,nopat) << ",";
714  }
715  }
716  ss << ")];" << endl;
717  ss << "res:=Det[m];" << endl;
718  auto tmp = ss.str();
719  string_replace_all(tmp,",)]",")]");
720  fermat.Execute(tmp);
721  ss.clear();
722  ss.str("");
723 
724  static string bstr("[-begin-]"), estr("[-end-]");
725  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
726  ss << "!('" <<bstr<< "',res,'" <<estr<< "')" << endl;
727  auto ostr = fermat.Execute(ss.str());
728  ss.clear();
729  ss.str("");
730 
731  // note the order,(normal_fermat will be called again in factor_form)
732  ss << "&(U=0);" << endl; // disable ugly printing
733  ss << "@res;" << endl;
734  ss << "@[**];" << endl;
735  ss << "&_G;" << endl;
736  fermat.Execute(ss.str());
737  ss.clear();
738  ss.str("");
739 
740  // make sure last char is 0
741  if(ostr[ostr.length()-1]!='0') throw Error("fermat_together: last char is NOT 0.");
742  ostr = ostr.substr(0, ostr.length()-1);
743  auto cpos = ostr.find(bstr);
744  if(cpos==string::npos) throw Error(bstr+" NOT Found.");
745  ostr = ostr.substr(cpos+bstr.length(),string::npos);
746  cpos = ostr.find(estr);
747  if(cpos==string::npos) throw Error(estr+" NOT Found.");
748  ostr = ostr.substr(0,cpos);
749  string_trim(ostr);
750 
751  Parser fp(st);
752  auto res = fp.Read(ostr);
753  res = res.subs(map_rat);
754  return res;
755  }
756 
757  ex fermat_Det_Sparse(const matrix & mat_in) {
758  Fermat &fermat = Fermat::get();
759  int &v_max = fermat.vmax;
760 
761  exmap map_rat;
762  int nrow = mat_in.rows();
763  int ncol = mat_in.cols();
764  matrix mat(nrow, ncol);
765  for(int r=0; r<nrow; r++) for(int c=0; c<ncol; c++) mat(r,c) = mat_in(r,c).to_rational(map_rat);
766 
767  lst rep_vs;
768  if(true) {
769  map<ex,long long,ex_is_less> s2c;
770  ex expr_in = mat;
771  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
772  if(is_a<symbol>(*i)) s2c[*i]++;
773  }
774  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
775  for(auto kv : s2c) {
776  auto fw = fermat_weight.find(kv.first);
777  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
778  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
779  else sv1.push_back(lst{kv.second, kv.first});
780  }
781  sort_vec(sv1);
782  sort_vec(sv2);
783  sort_vec(sv3);
784  for(auto sv : sv1) rep_vs.append(sv.op(1));
785  for(auto sv : sv2) rep_vs.append(sv.op(1));
786  for(auto sv : sv3) rep_vs.append(sv.op(1));
787  }
788 
789  exmap v2f;
790  symtab st;
791  int fvi = 0;
792  for(auto vi : rep_vs) {
793  auto name = "v" + to_string(fvi);
794  v2f[vi] = Symbol(name);
795  st[name] = vi;
796  fvi++;
797  }
798 
799  stringstream ss;
800  if(fvi>111) {
801  cout << rep_vs << endl;
802  throw Error("Fermat: Too many variables.");
803  }
804  if(fvi>v_max) {
805  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
806  fermat.Execute(ss.str());
807  ss.clear();
808  ss.str("");
809  v_max = fvi;
810  }
811 
812  ostringstream oss;
813  oss << "[";
814  bool r1 = true;
815  for(int r=0; r<nrow; r++) {
816  bool c1 = true;
817  for(int c=0; c<ncol; c++) {
818  if(mat(r,c).is_zero()) continue;
819  if(c1) {
820  if(!r1) oss << "`" << endl;
821  oss << "[" << r+1 << ",";
822  r1 = c1 = false;
823  } else oss << ",";
824  oss << "[" << c+1 << "," << mat(r,c).subs(v2f) << "]";
825  }
826  if(!c1) oss << "]";
827  }
828  oss << "];" << endl;
829 
830  ss << "Array m[" << nrow << "," << ncol << "] Sparse;" << endl;
831  fermat.Execute(ss.str());
832  ss.clear();
833  ss.str("");
834 
835  ss << "[m]:=" << oss.str();
836  fermat.Execute(ss.str());
837  ss.clear();
838  ss.str("");
839 
840  ss << "res:=Det[m];" << endl;
841  auto tmp = ss.str();
842  string_replace_all(tmp,",)]",")]");
843  fermat.Execute(tmp);
844  ss.clear();
845  ss.str("");
846 
847  static string bstr("[-begin-]"), estr("[-end-]");
848  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
849  ss << "!('" <<bstr<< "',res,'" <<estr<< "')" << endl;
850  auto ostr = fermat.Execute(ss.str());
851  ss.clear();
852  ss.str("");
853 
854  // note the order,(normal_fermat will be called again in factor_form)
855  ss << "&(U=0);" << endl; // disable ugly printing
856  ss << "@res;" << endl;
857  ss << "@[**];" << endl;
858  ss << "&_G;" << endl;
859  fermat.Execute(ss.str());
860  ss.clear();
861  ss.str("");
862 
863  // make sure last char is 0
864  if(ostr[ostr.length()-1]!='0') throw Error("fermat_together: last char is NOT 0.");
865  ostr = ostr.substr(0, ostr.length()-1);
866  auto cpos = ostr.find(bstr);
867  if(cpos==string::npos) throw Error(bstr+" NOT Found.");
868  ostr = ostr.substr(cpos+bstr.length(),string::npos);
869  cpos = ostr.find(estr);
870  if(cpos==string::npos) throw Error(estr+" NOT Found.");
871  ostr = ostr.substr(0,cpos);
872  string_trim(ostr);
873 
874  Parser fp(st);
875  auto res = fp.Read(ostr);
876  res = res.subs(map_rat);
877  return res;
878  }
879 
880  matrix fermat_inv(const matrix & mat_in) {
881  Fermat &fermat = Fermat::get();
882  int &v_max = fermat.vmax;
883 
884  exmap map_rat;
885  int nrow = mat_in.rows();
886  int ncol = mat_in.cols();
887  matrix mat(nrow, ncol);
888  for(int r=0; r<nrow; r++) for(int c=0; c<ncol; c++) mat(r,c) = mat_in(r,c).to_rational(map_rat);
889 
890  lst rep_vs;
891  if(true) {
892  map<ex,long long,ex_is_less> s2c;
893  ex expr_in = mat;
894  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
895  if(is_a<symbol>(*i)) s2c[*i]++;
896  }
897  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
898  for(auto kv : s2c) {
899  auto fw = fermat_weight.find(kv.first);
900  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
901  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
902  else sv1.push_back(lst{kv.second, kv.first});
903  }
904  sort_vec(sv1);
905  sort_vec(sv2);
906  sort_vec(sv3);
907  for(auto sv : sv1) rep_vs.append(sv.op(1));
908  for(auto sv : sv2) rep_vs.append(sv.op(1));
909  for(auto sv : sv3) rep_vs.append(sv.op(1));
910  }
911 
912  exmap v2f;
913  symtab st;
914  int fvi = 0;
915  for(auto vi : rep_vs) {
916  auto name = "v" + to_string(fvi);
917  v2f[vi] = Symbol(name);
918  st[name] = vi;
919  fvi++;
920  }
921 
922  stringstream ss;
923  if(fvi>111) {
924  cout << rep_vs << endl;
925  throw Error("Fermat: Too many variables.");
926  }
927  if(fvi>v_max) {
928  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
929  fermat.Execute(ss.str());
930  ss.clear();
931  ss.str("");
932  v_max = fvi;
933  }
934 
935  ss << "Array m[" << nrow << "," << ncol << "];" << endl;
936  fermat.Execute(ss.str());
937  ss.clear();
938  ss.str("");
939 
940  ss << "[m]:=[(";
941  for(int c=0; c<ncol; c++) {
942  for(int r=0; r<nrow; r++) {
943  ss << mat(r,c).subs(iEpsilon==0,nopat).subs(v2f,nopat) << ",";
944  }
945  }
946  ss << ")];" << endl;
947  ss << "[m]:=1/[m];" << endl;
948  auto tmp = ss.str();
949  string_replace_all(tmp,",)]",")]");
950  fermat.Execute(tmp);
951  ss.clear();
952  ss.str("");
953 
954  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
955  ss << "![m" << endl;
956  auto ostr = fermat.Execute(ss.str());
957  ss.clear();
958  ss.str("");
959  //fermat.Exit();
960 
961  // note the order, before exfactor (normal_fermat will be called again here)
962  ss << "&(U=0);" << endl; // disable ugly printing
963  ss << "@([m]);" << endl;
964  ss << "&_G;" << endl;
965  fermat.Execute(ss.str());
966  ss.clear();
967  ss.str("");
968 
969  // make sure last char is 0
970  if(ostr[ostr.length()-1]!='0') throw Error("Direc::Export, last char is NOT 0.");
971  ostr = ostr.substr(0, ostr.length()-1);
972  string_trim(ostr);
973 
974  ostr.erase(0, ostr.find(":=")+2);
975  string_replace_all(ostr, "[", "{");
976  string_replace_all(ostr, "]", "}");
977  Parser fp(st);
978  matrix mr(nrow, ncol);
979  auto res = fp.Read(ostr);
980  for(int r=0; r<nrow; r++) {
981  auto cur = res.op(r);
982  for(int c=0; c<ncol; c++) mr(r,c) = cur.op(c).subs(map_rat);
983  }
984  return mr;
985  }
986 
987  matrix fermat_mul(const matrix & m1, const matrix & m2) {
988  Fermat &fermat = Fermat::get();
989  int &v_max = fermat.vmax;
990 
991  exmap map_rat;
992  int nr1 = m1.rows();
993  int nc1 = m1.cols();
994  matrix mat1(nr1, nc1);
995  for(int r=0; r<nr1; r++) for(int c=0; c<nc1; c++) mat1(r,c) = m1(r,c).to_rational(map_rat);
996  int nr2 = m2.rows();
997  int nc2 = m2.cols();
998  matrix mat2(nr2, nc2);
999  for(int r=0; r<nr2; r++) for(int c=0; c<nc2; c++) mat2(r,c) = m2(r,c).to_rational(map_rat);
1000 
1001  lst rep_vs;
1002  if(true) {
1003  map<ex,long long,ex_is_less> s2c;
1004  ex expr_in = mat1;
1005  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
1006  if(is_a<symbol>(*i)) s2c[*i]++;
1007  }
1008  expr_in = mat2;
1009  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
1010  if(is_a<symbol>(*i)) s2c[*i]++;
1011  }
1012  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
1013  for(auto kv : s2c) {
1014  auto fw = fermat_weight.find(kv.first);
1015  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
1016  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
1017  else sv1.push_back(lst{kv.second, kv.first});
1018  }
1019  sort_vec(sv1);
1020  sort_vec(sv2);
1021  sort_vec(sv3);
1022  for(auto sv : sv1) rep_vs.append(sv.op(1));
1023  for(auto sv : sv2) rep_vs.append(sv.op(1));
1024  for(auto sv : sv3) rep_vs.append(sv.op(1));
1025  }
1026 
1027  exmap v2f;
1028  symtab st;
1029  int fvi = 0;
1030  for(auto vi : rep_vs) {
1031  auto name = "v" + to_string(fvi);
1032  v2f[vi] = Symbol(name);
1033  st[name] = vi;
1034  fvi++;
1035  }
1036 
1037  stringstream ss;
1038  if(fvi>111) {
1039  cout << rep_vs << endl;
1040  throw Error("Fermat: Too many variables.");
1041  }
1042  if(fvi>v_max) {
1043  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
1044  fermat.Execute(ss.str());
1045  ss.clear();
1046  ss.str("");
1047  v_max = fvi;
1048  }
1049 
1050  ss << "Array m1[" << nr1 << "," << nc1 << "];" << endl;
1051  ss << "Array m2[" << nr2 << "," << nc2 << "];" << endl;
1052  fermat.Execute(ss.str());
1053  ss.clear();
1054  ss.str("");
1055 
1056  ss << "[m1]:=[(";
1057  for(int c=0; c<nc1; c++) {
1058  for(int r=0; r<nr1; r++) {
1059  ss << mat1(r,c).subs(iEpsilon==0,nopat).subs(v2f,nopat) << ",";
1060  }
1061  }
1062  ss << ")];" << endl;
1063 
1064  ss << "[m2]:=[(";
1065  for(int c=0; c<nc2; c++) {
1066  for(int r=0; r<nr2; r++) {
1067  ss << mat2(r,c).subs(iEpsilon==0,nopat).subs(v2f,nopat) << ",";
1068  }
1069  }
1070  ss << ")];" << endl;
1071 
1072  ss << "[m]:=[m1]*[m2];" << endl;
1073  auto tmp = ss.str();
1074  string_replace_all(tmp,",)]",")]");
1075  fermat.Execute(tmp);
1076  ss.clear();
1077  ss.str("");
1078 
1079  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
1080  ss << "![m" << endl;
1081  auto ostr = fermat.Execute(ss.str());
1082  ss.clear();
1083  ss.str("");
1084  //fermat.Exit();
1085 
1086  // note the order, before exfactor (normal_fermat will be called again here)
1087  ss << "&(U=0);" << endl; // disable ugly printing
1088  ss << "@([m1],[m2],[m]);" << endl;
1089  ss << "&_G;" << endl;
1090  fermat.Execute(ss.str());
1091  ss.clear();
1092  ss.str("");
1093 
1094  // make sure last char is 0
1095  if(ostr[ostr.length()-1]!='0') throw Error("Direc::Export, last char is NOT 0.");
1096  ostr = ostr.substr(0, ostr.length()-1);
1097  string_trim(ostr);
1098 
1099  ostr.erase(0, ostr.find(":=")+2);
1100  string_replace_all(ostr, "[", "{");
1101  string_replace_all(ostr, "]", "}");
1102  Parser fp(st);
1103  matrix mr(nr1, nc2);
1104  auto res = fp.Read(ostr);
1105  for(int r=0; r<nr1; r++) {
1106  auto cur = res.op(r);
1107  for(int c=0; c<nc2; c++) mr(r,c) = cur.op(c).subs(map_rat);
1108  }
1109  return mr;
1110  }
1111 
1112  matrix fermat_pow(const matrix & mat_in, int n) {
1113  Fermat &fermat = Fermat::get();
1114  int &v_max = fermat.vmax;
1115 
1116  exmap map_rat;
1117  int nrow = mat_in.rows();
1118  int ncol = mat_in.cols();
1119  matrix mat(nrow, ncol);
1120  for(int r=0; r<nrow; r++) for(int c=0; c<ncol; c++) mat(r,c) = mat_in(r,c).to_rational(map_rat);
1121 
1122  lst rep_vs;
1123  if(true) {
1124  map<ex,long long,ex_is_less> s2c;
1125  ex expr_in = mat;
1126  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
1127  if(is_a<symbol>(*i)) s2c[*i]++;
1128  }
1129  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
1130  for(auto kv : s2c) {
1131  auto fw = fermat_weight.find(kv.first);
1132  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
1133  else if(map_rat.find(kv.first)!=map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
1134  else sv1.push_back(lst{kv.second, kv.first});
1135  }
1136  sort_vec(sv1);
1137  sort_vec(sv2);
1138  sort_vec(sv3);
1139  for(auto sv : sv1) rep_vs.append(sv.op(1));
1140  for(auto sv : sv2) rep_vs.append(sv.op(1));
1141  for(auto sv : sv3) rep_vs.append(sv.op(1));
1142  }
1143 
1144  exmap v2f;
1145  symtab st;
1146  int fvi = 0;
1147  for(auto vi : rep_vs) {
1148  auto name = "v" + to_string(fvi);
1149  v2f[vi] = Symbol(name);
1150  st[name] = vi;
1151  fvi++;
1152  }
1153 
1154  stringstream ss;
1155  if(fvi>111) {
1156  cout << rep_vs << endl;
1157  throw Error("Fermat: Too many variables.");
1158  }
1159  if(fvi>v_max) {
1160  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
1161  fermat.Execute(ss.str());
1162  ss.clear();
1163  ss.str("");
1164  v_max = fvi;
1165  }
1166 
1167  ss << "Array m[" << nrow << "," << ncol << "];" << endl;
1168  fermat.Execute(ss.str());
1169  ss.clear();
1170  ss.str("");
1171 
1172  ss << "[m]:=[(";
1173  for(int c=0; c<ncol; c++) {
1174  for(int r=0; r<nrow; r++) {
1175  ss << mat(r,c).subs(iEpsilon==0,nopat).subs(v2f,nopat) << ",";
1176  }
1177  }
1178  ss << ")];" << endl;
1179  ss << "[m]:=[m]^(" << n << ");" << endl;
1180  auto tmp = ss.str();
1181  string_replace_all(tmp,",)]",")]");
1182  fermat.Execute(tmp);
1183  ss.clear();
1184  ss.str("");
1185 
1186  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
1187  ss << "![m" << endl;
1188  auto ostr = fermat.Execute(ss.str());
1189  ss.clear();
1190  ss.str("");
1191  //fermat.Exit();
1192 
1193  // note the order, before exfactor (normal_fermat will be called again here)
1194  ss << "&(U=0);" << endl; // disable ugly printing
1195  ss << "@([m]);" << endl;
1196  ss << "&_G;" << endl;
1197  fermat.Execute(ss.str());
1198  ss.clear();
1199  ss.str("");
1200 
1201  // make sure last char is 0
1202  if(ostr[ostr.length()-1]!='0') throw Error("Direc::Export, last char is NOT 0.");
1203  ostr = ostr.substr(0, ostr.length()-1);
1204  string_trim(ostr);
1205 
1206  ostr.erase(0, ostr.find(":=")+2);
1207  string_replace_all(ostr, "[", "{");
1208  string_replace_all(ostr, "]", "}");
1209  Parser fp(st);
1210  matrix mr(nrow, ncol);
1211  auto res = fp.Read(ostr);
1212  for(int r=0; r<nrow; r++) {
1213  auto cur = res.op(r);
1214  for(int c=0; c<ncol; c++) mr(r,c) = cur.op(c).subs(map_rat);
1215  }
1216  return mr;
1217  }
1218 
1219  namespace {
1220  exmap mat_map_rat;
1221  symtab mat_st;
1222  }
1223 
1224  void fermat_mat(const matrix & mat_in, const string & _name) {
1225  string name = _name;
1226  string_replace_all(name, "[", "");
1227  string_replace_all(name, "]", "");
1228 
1229  Fermat &fermat = Fermat::get();
1230  int &v_max = fermat.vmax;
1231 
1232  int nrow = mat_in.rows();
1233  int ncol = mat_in.cols();
1234  matrix mat(nrow, ncol);
1235  for(int r=0; r<nrow; r++) for(int c=0; c<ncol; c++) mat(r,c) = mat_in(r,c).to_rational(mat_map_rat);
1236 
1237  lst rep_vs;
1238  if(true) {
1239  map<ex,long long,ex_is_less> s2c;
1240  ex expr_in = mat;
1241  for(const_preorder_iterator i = expr_in.preorder_begin(); i != expr_in.preorder_end(); ++i) {
1242  if(is_a<symbol>(*i)) s2c[*i]++;
1243  }
1244  exvector sv1, sv2, sv3; // sv2:fermat_weight, sv3:map_rat
1245  for(auto kv : s2c) {
1246  auto fw = fermat_weight.find(kv.first);
1247  if(fw!=fermat_weight.end()) sv2.push_back(lst{fw->second, fw->first});
1248  else if(mat_map_rat.find(kv.first)!=mat_map_rat.end()) sv3.push_back(lst{kv.second, kv.first});
1249  else sv1.push_back(lst{kv.second, kv.first});
1250  }
1251  sort_vec(sv1);
1252  sort_vec(sv2);
1253  sort_vec(sv3);
1254  for(auto sv : sv1) rep_vs.append(sv.op(1));
1255  for(auto sv : sv2) rep_vs.append(sv.op(1));
1256  for(auto sv : sv3) rep_vs.append(sv.op(1));
1257  }
1258 
1259  exmap v2f;
1260  int fvi = 0;
1261  for(auto vi : rep_vs) {
1262  auto name = "v" + to_string(fvi);
1263  v2f[vi] = Symbol(name);
1264  mat_st[name] = vi;
1265  fvi++;
1266  }
1267 
1268  stringstream ss;
1269  if(fvi>111) {
1270  cout << rep_vs << endl;
1271  throw Error("Fermat: Too many variables.");
1272  }
1273  if(fvi>v_max) {
1274  for(int i=v_max; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
1275  fermat.Execute(ss.str());
1276  ss.clear();
1277  ss.str("");
1278  v_max = fvi;
1279  }
1280 
1281  ss << "Array " << name << "[" << nrow << "," << ncol << "];" << endl;
1282  fermat.Execute(ss.str());
1283  ss.clear();
1284  ss.str("");
1285 
1286  ss << "["<<name<<"]:=[(";
1287  for(int c=0; c<ncol; c++) {
1288  for(int r=0; r<nrow; r++) {
1289  ss << mat(r,c).subs(iEpsilon==0,nopat).subs(v2f,nopat) << ",";
1290  }
1291  }
1292  ss << ")];" << endl;
1293  //ss << "Redrowech([m]);" << endl;
1294  auto tmp = ss.str();
1295  string_replace_all(tmp,",)]",")]");
1296  fermat.Execute(tmp);
1297  ss.clear();
1298  ss.str("");
1299  }
1300 
1301  matrix fermat_mat(const string & _name) {
1302  string name = _name;
1303  string_replace_all(name, "[", "");
1304  string_replace_all(name, "]", "");
1305 
1306  Fermat &fermat = Fermat::get();
1307  int &v_max = fermat.vmax;
1308 
1309  stringstream ss;
1310  ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
1311  ss << "![" << name << endl;
1312  auto ostr = fermat.Execute(ss.str());
1313  ss.clear();
1314  ss.str("");
1315  //fermat.Exit();
1316 
1317  // note the order, before exfactor (normal_fermat will be called again here)
1318  ss << "&(U=0);" << endl; // disable ugly printing
1319  ss << "&_G;" << endl;
1320  fermat.Execute(ss.str());
1321  ss.clear();
1322  ss.str("");
1323 
1324  // make sure last char is 0
1325  if(ostr[ostr.length()-1]!='0') throw Error("Direc::Export, last char is NOT 0.");
1326  ostr = ostr.substr(0, ostr.length()-1);
1327  string_trim(ostr);
1328  ostr.erase(0, ostr.find(":=")+2);
1329  string_replace_all(ostr, "[", "{");
1330  string_replace_all(ostr, "]", "}");
1331  Parser fp(mat_st);
1332  auto res = fp.Read(ostr);
1333  int nrow = res.nops();
1334  int ncol = res.op(0).nops();
1335  matrix mr(nrow, ncol);
1336  for(int r=0; r<nrow; r++) {
1337  auto cur = res.op(r);
1338  for(int c=0; c<ncol; c++) mr(r,c) = cur.op(c).subs(mat_map_rat);
1339  }
1340  return mr;
1341  }
1342 
1343  void fermat_eval(const string & fcmd) {
1344  Fermat &fermat = Fermat::get();
1345  int &v_max = fermat.vmax;
1346 
1347  stringstream ss;
1348  ss << fcmd << endl; // ugly printing, the whitespace matters
1349  auto ostr = fermat.Execute(ss.str());
1350  ss.clear();
1351  ss.str("");
1352  //fermat.Exit();
1353  }
1354 
1355 }
BASIC.h
Basic header file.
HepLib::Error
class used to wrap error message
Definition: BASIC.h:242
HepLib::Fermat
interface to communicate with Fermat program
Definition: BASIC.h:797
HepLib::Fermat::get
static Fermat & get()
Definition: Fermat.cpp:7
HepLib::Fermat::Execute
string Execute(string)
Definition: Process.cpp:102
HepLib::Parser
class to parse for string or file, helpful with Symbol class
Definition: BASIC.h:645
HepLib::Parser::Read
ex Read(const string &instr, bool s2S=true)
Definition: Functions.cpp:111
HepLib::Symbol
class extended to GiNaC symbol class, represent a positive symbol
Definition: BASIC.h:113
HepLib::Symbol::subs
ex subs(const exmap &m, unsigned options=0) const override
Definition: BASIC.h:145
HepLib
HepLib namespace.
Definition: BASIC.cpp:17
HepLib::fermat_Redrowech_Sparse
matrix fermat_Redrowech_Sparse(const matrix &mat)
Definition: Fermat.cpp:516
HepLib::iEpsilon
const iSymbol iEpsilon
HepLib::fermat_mul
matrix fermat_mul(const matrix &m1, const matrix &m2)
Definition: Fermat.cpp:987
HepLib::fermat_weight
map< ex, long long, ex_is_less > fermat_weight
Definition: Init.cpp:152
HepLib::factor_form
ex factor_form(const ex &expr, bool nd)
factorize a expression using FORM
Definition: BASIC.cpp:2038
HepLib::fermat_mat
void fermat_mat(const matrix &mat_in, const string &name)
Definition: Fermat.cpp:1224
HepLib::numer_denom_fermat
ex numer_denom_fermat(const ex &expr, bool dfactor=false)
return the numerator and denominator after normalization
Definition: Fermat.cpp:23
HepLib::fermat_using_array
int fermat_using_array
Definition: Init.cpp:151
HepLib::fermat_Redrowech
matrix fermat_Redrowech(const matrix &mat)
Definition: Fermat.cpp:409
HepLib::nextprime
ex nextprime(const ex &n)
Definition: BASIC.cpp:2274
HepLib::fermat_eval
ex fermat_eval(const ex &expr)
return the numerator and denominator after normalization
Definition: Fermat.cpp:317
HepLib::fermat_Det
ex fermat_Det(const matrix &mat)
Definition: Fermat.cpp:650
HepLib::fermat_inv
matrix fermat_inv(const matrix &mat)
Definition: Fermat.cpp:880
HepLib::sort_vec
void sort_vec(exvector &ivec, bool less=true)
sort the list in less order, or the reverse
Definition: Sort.cpp:54
HepLib::string_replace_all
void string_replace_all(string &str, const string &from, const string &to)
Definition: Functions.cpp:148
HepLib::nopat
unsigned nopat
Definition: Init.cpp:88
HepLib::fermat_Det_Sparse
ex fermat_Det_Sparse(const matrix &mat)
Definition: Fermat.cpp:757
HepLib::string_trim
void string_trim(string &str)
Definition: Functions.cpp:156
HepLib::numer_fermat
ex numer_fermat(const ex &expr)
Definition: Fermat.cpp:170
HepLib::fermat_pow
matrix fermat_pow(const matrix &mat_in, int n)
Definition: Fermat.cpp:1112
HepLib::ex2int
int ex2int(ex num)
ex to integer
Definition: BASIC.cpp:893
HepLib::subs
ex subs(const ex &e, init_list sl)
Definition: BASIC.h:51
Generated by doxygen 1.9.1