HepLib
Loading...
Searching...
No Matches
Fermat.cpp
Go to the documentation of this file.
1// Functions related to fermat program
2
3#include "BASIC.h"
4
5namespace 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
HepLib namespace.
Definition BASIC.cpp:17
HepLib::fermat_Redrowech_Sparse
matrix fermat_Redrowech_Sparse(const matrix &mat)
Definition Fermat.cpp:516
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:155
HepLib::factor_form
ex factor_form(const ex &expr, bool nd)
factorize a expression using FORM
Definition BASIC.cpp:2048
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:154
HepLib::fermat_Redrowech
matrix fermat_Redrowech(const matrix &mat)
Definition Fermat.cpp:409
HepLib::nextprime
ex nextprime(const ex &n)
Definition BASIC.cpp:2284
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::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.8