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