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