- Banned exZERO(), exONE(), exMINUSHALF() and all this from the interface.
[ginac.git] / ginac / simp_lor.cpp
1 /** @file simp_lor.cpp
2  *
3  *  Implementation of GiNaC's simp_lor objects.
4  *  No real implementation yet, to be done.     */
5
6 /*
7  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
8  *
9  *  This program is free software; you can redistribute it and/or modify
10  *  it under the terms of the GNU General Public License as published by
11  *  the Free Software Foundation; either version 2 of the License, or
12  *  (at your option) any later version.
13  *
14  *  This program is distributed in the hope that it will be useful,
15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17  *  GNU General Public License for more details.
18  *
19  *  You should have received a copy of the GNU General Public License
20  *  along with this program; if not, write to the Free Software
21  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
22  */
23
24 #include <string>
25 #include <list>
26 #include <algorithm>
27 #include <iostream>
28 #include <stdexcept>
29 #include <map>
30
31 #include "simp_lor.h"
32 #include "ex.h"
33 #include "mul.h"
34 #include "symbol.h"
35 #include "debugmsg.h"
36 #include "utils.h"
37
38 #ifndef NO_GINAC_NAMESPACE
39 namespace GiNaC {
40 #endif // ndef NO_GINAC_NAMESPACE
41
42 //////////
43 // default constructor, destructor, copy constructor assignment operator and helpers
44 //////////
45
46 // public
47
48 simp_lor::simp_lor() : type(invalid)
49 {
50     debugmsg("simp_lor default constructor",LOGLEVEL_CONSTRUCT);
51     tinfo_key=TINFO_simp_lor;
52 }
53
54 simp_lor::~simp_lor()
55 {
56     debugmsg("simp_lor destructor",LOGLEVEL_DESTRUCT);
57     destroy(0);
58 }
59
60 simp_lor::simp_lor(simp_lor const & other)
61 {
62     debugmsg("simp_lor copy constructor",LOGLEVEL_CONSTRUCT);
63     copy (other);
64 }
65
66 simp_lor const & simp_lor::operator=(simp_lor const & other)
67 {
68     debugmsg("simp_lor operator=",LOGLEVEL_ASSIGNMENT);
69     if (this != &other) {
70         destroy(1);
71         copy(other);
72     }
73     return *this;
74 }
75
76 // protected
77
78 void simp_lor::copy(simp_lor const & other)
79 {
80     indexed::copy(other);
81     type=other.type;
82     name=other.name;
83 }
84
85 void simp_lor::destroy(bool call_parent)
86 {
87     if (call_parent) {
88         indexed::destroy(call_parent);
89     }
90 }
91
92 //////////
93 // other constructors
94 //////////
95
96 // protected
97
98 simp_lor::simp_lor(simp_lor_types const t) : type(t)
99 {
100     debugmsg("simp_lor constructor from simp_lor_types",LOGLEVEL_CONSTRUCT);
101     tinfo_key=TINFO_simp_lor;
102 }
103
104 simp_lor::simp_lor(simp_lor_types const t, ex const & i1, ex const & i2) :
105     indexed(i1,i2), type(t)
106 {
107     debugmsg("simp_lor constructor from simp_lor_types,ex,ex",LOGLEVEL_CONSTRUCT);
108     tinfo_key=TINFO_simp_lor;
109     GINAC_ASSERT(all_of_type_lorentzidx());
110 }
111
112 simp_lor::simp_lor(simp_lor_types const t, string const & n, ex const & i1) :
113     indexed(i1), type(t), name(n)
114 {
115     debugmsg("simp_lor constructor from simp_lor_types,string,ex",LOGLEVEL_CONSTRUCT);
116     tinfo_key=TINFO_simp_lor;
117     GINAC_ASSERT(all_of_type_lorentzidx());
118 }
119
120 simp_lor::simp_lor(simp_lor_types const t, string const & n, exvector const & iv) :
121     indexed(iv), type(t), name(n)
122 {
123     debugmsg("simp_lor constructor from simp_lor_types,string,exvector",LOGLEVEL_CONSTRUCT);
124     tinfo_key=TINFO_simp_lor;
125     GINAC_ASSERT(all_of_type_lorentzidx());
126 }
127
128 simp_lor::simp_lor(simp_lor_types const t, string const & n, exvector * ivp) :
129     indexed(ivp), type(t), name(n)
130 {
131     debugmsg("simp_lor constructor from simp_lor_types,string,exvector*",LOGLEVEL_CONSTRUCT);
132     tinfo_key=TINFO_simp_lor;
133     GINAC_ASSERT(all_of_type_lorentzidx());
134 }
135
136 //////////
137 // functions overriding virtual functions from bases classes
138 //////////
139
140 // public
141
142 basic * simp_lor::duplicate() const
143 {
144     debugmsg("simp_lor duplicate",LOGLEVEL_DUPLICATE);
145     return new simp_lor(*this);
146 }
147
148 void simp_lor::printraw(ostream & os) const
149 {
150     debugmsg("simp_lor printraw",LOGLEVEL_PRINT);
151     os << "simp_lor(type=" << (unsigned)type
152        << ",name=" << name << ",indices=";
153     printrawindices(os);
154     os << ",hash=" << hashvalue << ",flags=" << flags << ")";
155 }
156
157 void simp_lor::printtree(ostream & os, unsigned indent) const
158 {
159     debugmsg("simp_lor printtree",LOGLEVEL_PRINT);
160     os << string(indent,' ') << "simp_lor object: "
161        << "type=" << (unsigned)type
162        << ", name=" << name << ", ";
163     os << seq.size() << " indices" << endl;
164     printtreeindices(os,indent);
165     os << string(indent,' ') << "hash=" << hashvalue
166        << " (0x" << hex << hashvalue << dec << ")"
167        << ", flags=" << flags << endl;
168 }
169
170 void simp_lor::print(ostream & os, unsigned upper_precedence) const
171 {
172     debugmsg("simp_lor print",LOGLEVEL_PRINT);
173     switch (type) {
174     case simp_lor_g:
175         os << "g";
176         break;
177     case simp_lor_vec:
178         os << name;
179         break;
180     case invalid:
181     default:
182         os << "INVALID_SIMP_LOR_OBJECT";
183         break;
184     }
185     printindices(os);
186 }
187
188 void simp_lor::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
189 {
190     debugmsg("simp_lor print csrc",LOGLEVEL_PRINT);
191     print(os,upper_precedence);
192 }
193
194 bool simp_lor::info(unsigned inf) const
195 {
196     return indexed::info(inf);
197 }
198
199 ex simp_lor::eval(int level) const
200 {
201     if (type==simp_lor_g) {
202         // canonicalize indices
203         exvector iv=seq;
204         int sig=canonicalize_indices(iv,false); // symmetric
205         if (sig!=INT_MAX) {
206             // something has changed while sorting indices, more evaluations later
207             if (sig==0) return _ex0();
208             return ex(sig)*simp_lor(type,name,iv);
209         }
210         lorentzidx const & idx1=ex_to_lorentzidx(seq[0]);
211         lorentzidx const & idx2=ex_to_lorentzidx(seq[1]);
212         if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
213             // both indices are numeric
214             if ((idx1.get_value()==idx2.get_value())) {
215                 // both on diagonal
216                 if (idx1.get_value()==0) {
217                     // (0,0)
218                     return _ex1();
219                 } else {
220                     if (idx1.is_covariant()!=idx2.is_covariant()) {
221                         // (_i,~i) or (~i,_i), i=1..3
222                         return _ex1();
223                     } else {
224                         // (_i,_i) or (~i,~i), i=1..3
225                         return _ex_1();
226                     }
227                 }
228             } else {
229                 // at least one off-diagonal
230                 return _ex0();
231             }
232         } else if (idx1.is_symbolic() &&
233                    idx1.is_co_contra_pair(idx2)) {
234             return Dim()-idx1.get_dim_parallel_space();
235         }
236     }
237
238     return this->hold();
239 }
240     
241 // protected
242
243 int simp_lor::compare_same_type(basic const & other) const
244 {
245     GINAC_ASSERT(other.tinfo() == TINFO_simp_lor);
246     const simp_lor *o = static_cast<const simp_lor *>(&other);
247     if (type==o->type) {
248         if (name==o->name) {
249             return indexed::compare_same_type(other);
250         }
251         return name.compare(o->name);
252     }
253     return type < o->type ? -1 : 1;
254 }
255
256 bool simp_lor::is_equal_same_type(basic const & other) const
257 {
258     GINAC_ASSERT(other.tinfo() == TINFO_simp_lor);
259     const simp_lor *o = static_cast<const simp_lor *>(&other);
260     if (type!=o->type) return false;
261     if (name!=o->name) return false;
262     return indexed::is_equal_same_type(other);
263 }
264
265 unsigned simp_lor::return_type(void) const
266 {
267     return return_types::commutative;
268 }
269    
270 unsigned simp_lor::return_type_tinfo(void) const
271 {
272     return tinfo_key;
273 }
274
275 ex simp_lor::thisexprseq(exvector const & v) const
276 {
277     return simp_lor(type,name,v);
278 }
279
280 ex simp_lor::thisexprseq(exvector * vp) const
281 {
282     return simp_lor(type,name,vp);
283 }
284
285 //////////
286 // virtual functions which can be overridden by derived classes
287 //////////
288
289 // none
290
291 //////////
292 // non-virtual functions in this class
293 //////////
294
295 // protected
296
297 bool simp_lor::all_of_type_lorentzidx(void) const
298 {
299     // used only inside of ASSERTs
300     for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
301         if (!is_ex_of_type(*cit,lorentzidx)) return false;
302     }
303     return true;
304 }
305
306 //////////
307 // static member variables
308 //////////
309
310 // none
311
312 //////////
313 // global constants
314 //////////
315
316 const simp_lor some_simp_lor;
317 type_info const & typeid_simp_lor=typeid(some_simp_lor);
318
319 //////////
320 // friend functions
321 //////////
322
323 simp_lor lor_g(ex const & mu, ex const & nu)
324 {
325     return simp_lor(simp_lor::simp_lor_g,mu,nu);
326 }
327
328 simp_lor lor_vec(string const & n, ex const & mu)
329 {
330     return simp_lor(simp_lor::simp_lor_vec,n,mu);
331 }
332
333 ex simplify_simp_lor_mul(ex const & m, scalar_products const & sp)
334 {
335     GINAC_ASSERT(is_ex_exactly_of_type(m,mul));
336     exvector v_contracted;
337
338     // collect factors in an exvector, store squares twice
339     int n=m.nops();
340     v_contracted.reserve(2*n);
341     for (int i=0; i<n; ++i) {
342         ex f=m.op(i);
343         if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
344             v_contracted.push_back(f.op(0));
345             v_contracted.push_back(f.op(0));
346         } else {
347             v_contracted.push_back(f);
348         }
349     }
350
351     unsigned replacements;
352     bool something_changed=false;
353
354     exvector::iterator it=v_contracted.begin();
355     while (it!=v_contracted.end()) {
356         // process only lor_g objects
357         if (is_ex_exactly_of_type(*it,simp_lor) &&
358             (ex_to_simp_lor(*it).type==simp_lor::simp_lor_g)) {
359             simp_lor const & g=ex_to_simp_lor(*it);
360             GINAC_ASSERT(g.seq.size()==2);
361             idx const & first_idx=ex_to_lorentzidx(g.seq[0]);
362             idx const & second_idx=ex_to_lorentzidx(g.seq[1]);
363             // g_{mu,mu} should have been contracted in simp_lor::eval()
364             GINAC_ASSERT(!first_idx.is_equal(second_idx));
365             ex saved_g=*it; // save to restore it later
366
367             // try to contract first index
368             replacements=0;
369             if (first_idx.is_symbolic()) {
370                 replacements = subs_index_in_exvector(v_contracted,
371                                    first_idx.toggle_covariant(),second_idx);
372                 if (replacements==0) {
373                     // not contracted, restore g object
374                     *it=saved_g;
375                 } else {
376                     // a contracted index should occur exactly once
377                     GINAC_ASSERT(replacements==1);
378                     *it=_ex1();
379                     something_changed=true;
380                 }
381             }
382
383             // try second index only if first was not contracted
384             if ((replacements==0)&&(second_idx.is_symbolic())) {
385                 // first index not contracted, *it is again the original g object
386                 replacements = subs_index_in_exvector(v_contracted,
387                                    second_idx.toggle_covariant(),first_idx);
388                 if (replacements==0) {
389                     // not contracted except in itself, restore g object
390                     *it=saved_g;
391                 } else {
392                     // a contracted index should occur exactly once
393                     GINAC_ASSERT(replacements==1);
394                     *it=_ex1();
395                     something_changed=true;
396                 }
397             }
398         }
399         ++it;
400     }
401
402     // process only lor_vec objects
403     bool jump_to_next=false;
404     exvector::iterator it1=v_contracted.begin();
405     while (it1!=v_contracted.end()-1) {
406         if (is_ex_exactly_of_type(*it1,simp_lor) && 
407             (ex_to_simp_lor(*it1).type==simp_lor::simp_lor_vec)) {
408             exvector::iterator it2=it1+1;
409             while ((it2!=v_contracted.end())&&!jump_to_next) {
410                 if (is_ex_exactly_of_type(*it2,simp_lor) && 
411                     (ex_to_simp_lor(*it2).type==simp_lor::simp_lor_vec)) {
412                     simp_lor const & vec1=ex_to_simp_lor(*it1);
413                     simp_lor const & vec2=ex_to_simp_lor(*it2);
414                     GINAC_ASSERT(vec1.seq.size()==1);
415                     GINAC_ASSERT(vec2.seq.size()==1);
416                     lorentzidx const & idx1=ex_to_lorentzidx(vec1.seq[0]);
417                     lorentzidx const & idx2=ex_to_lorentzidx(vec2.seq[0]);
418                     if (idx1.is_symbolic() &&
419                         idx1.is_co_contra_pair(idx2) &&
420                         sp.is_defined(vec1,vec2)) {
421                         *it1=sp.evaluate(vec1,vec2);
422                         *it2=_ex1();
423                         something_changed=true;
424                         jump_to_next=true;
425                     }
426                 }
427                 ++it2;
428             }
429             jump_to_next=false;
430         }
431         ++it1;
432     }
433     if (something_changed) {
434         return mul(v_contracted);
435     }
436     return m;
437 }
438
439 ex simplify_simp_lor(ex const & e, scalar_products const & sp)
440 {
441     // all simplification is done on expanded objects
442     ex e_expanded=e.expand();
443
444     // simplification of sum=sum of simplifications
445     if (is_ex_exactly_of_type(e_expanded,add)) {
446         ex sum=_ex0();
447         for (int i=0; i<e_expanded.nops(); ++i) {
448             sum += simplify_simp_lor(e_expanded.op(i),sp);
449         }
450         return sum;
451     }
452
453     // simplification of commutative product=commutative product of simplifications
454     if (is_ex_exactly_of_type(e_expanded,mul)) {
455         return simplify_simp_lor_mul(e,sp);
456     }
457
458     // cannot do anything
459     return e_expanded;
460 }
461
462 //ex Dim(void)   // FIXME: what's going on here?
463 //{
464 //    static symbol * d=new symbol("dim");
465 //    return *d;
466 //}
467
468 //////////
469 // helper classes
470 //////////
471
472 void scalar_products::reg(simp_lor const & v1, simp_lor const & v2,
473                           ex const & sp)
474 {
475     if (v1.compare_same_type(v2)>0) {
476         reg(v2,v1,sp);
477         return;
478     }
479     spm[make_key(v1,v2)]=sp;
480 }
481
482 bool scalar_products::is_defined(simp_lor const & v1, simp_lor const & v2) const
483 {
484     if (v1.compare_same_type(v2)>0) {
485         return is_defined(v2,v1);
486     }
487     return spm.find(make_key(v1,v2))!=spm.end();
488 }
489
490 ex scalar_products::evaluate(simp_lor const & v1, simp_lor const & v2) const
491 {
492     if (v1.compare_same_type(v2)>0) {
493         return evaluate(v2,v1);
494     }
495     return spm.find(make_key(v1,v2))->second;
496 }
497
498 void scalar_products::debugprint(void) const
499 {
500     cerr << "map size=" << spm.size() << endl;
501     for (spmap::const_iterator cit=spm.begin(); cit!=spm.end(); ++cit) {
502         spmapkey const & k=(*cit).first;
503         cerr << "item key=((" << k.first.first
504              << "," << k.first.second << "),";
505         k.second.printraw(cerr);
506         cerr << ") value=" << (*cit).second << endl;
507     }
508 }
509
510 spmapkey scalar_products::make_key(simp_lor const & v1, simp_lor const & v2)
511 {
512     GINAC_ASSERT(v1.type==simp_lor::simp_lor_vec);
513     GINAC_ASSERT(v2.type==simp_lor::simp_lor_vec);
514     lorentzidx anon=ex_to_lorentzidx(v1.seq[0]).create_anonymous_representative();
515     GINAC_ASSERT(anon.is_equal_same_type(ex_to_lorentzidx(v2.seq[0]).create_anonymous_representative()));
516     return spmapkey(strstrpair(v1.name,v2.name),anon);
517 }
518
519 #ifndef NO_GINAC_NAMESPACE
520 } // namespace GiNaC
521 #endif // ndef NO_GINAC_NAMESPACE