Merge git://www.ginac.de/ginac
[ginac.git] / check / time_antipode.cpp
1 /** @file time_antipode.cpp
2  *
3  *  This is a beautiful example that calculates the counterterm for the
4  *  overall divergence of some special sorts of Feynman diagrams in a massless
5  *  Yukawa theory.  For this end it computes the antipode of the corresponding
6  *  decorated rooted tree using dimensional regularization in the parameter
7  *  x==-(D-4)/2, which leads to a Laurent series in x.  The renormalization
8  *  scheme used is the minimal subtraction scheme (MS).  From an efficiency
9  *  point of view it boils down to power series expansion.  It also has quite
10  *  an intriguing check for consistency, which is why we include it here.
11  *
12  *  This program is based on work by
13  *      Isabella Bierenbaum <bierenbaum@thep.physik.uni-mainz.de> and
14  *      Dirk Kreimer <dkreimer@bu.edu>.
15  *  For details, please see <http://www.arXiv.org/abs/hep-th/0111192>.
16  */
17
18 /*
19  *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
20  *
21  *  This program is free software; you can redistribute it and/or modify
22  *  it under the terms of the GNU General Public License as published by
23  *  the Free Software Foundation; either version 2 of the License, or
24  *  (at your option) any later version.
25  *
26  *  This program is distributed in the hope that it will be useful,
27  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
28  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
29  *  GNU General Public License for more details.
30  *
31  *  You should have received a copy of the GNU General Public License
32  *  along with this program; if not, write to the Free Software
33  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
34  */
35
36 #include <utility>
37 #include <vector>
38 #include <set>
39 #include <map>
40 #include <typeinfo>
41 #include <stdexcept>
42 #include "timer.h"
43 #include "ginac.h"
44 using namespace std;
45 using namespace GiNaC;
46
47 // whether to run this beast or not:
48 static const bool do_test = true;
49
50 // regularization parameter:
51 static const symbol x("x");
52
53 typedef pair<unsigned, unsigned> ijpair;
54 typedef pair<class node, bool> child;
55
56 const constant TrOne("Tr[One]", numeric(4));
57
58 /* Extract only the divergent part of a series and discard the rest. */
59 static ex div_part(const ex &exarg, const symbol &x, unsigned grad)
60 {
61         const ex exser = exarg.series(x==0, grad);
62         if (exser.degree(x)<0)
63                 throw runtime_error("divergent part truncation disaster");
64         ex exser_trunc;
65         for (int i=exser.ldegree(x); i<0; ++i)
66                 exser_trunc += exser.coeff(x,i)*pow(x,i);
67         // NB: exser_trunc is by construction collected in x.
68         return exser_trunc;
69 }
70
71 /* F_ab(a, i, b, j, "x") is a common pattern in all vertex evaluators. */
72 static ex F_ab(int a, int i, int b, int j, const symbol &x)
73 {
74         if ((i==0 && a<=0) || (j==0 && b<=0))
75                 return 0;
76         else
77                 return (tgamma(2-a-(i+1)*x)*
78                         tgamma(2-b-(1+j)*x)*
79                         tgamma(a+b-2+(1+i+j)*x)/
80                         tgamma(a+i*x)/
81                         tgamma(b+j*x)/tgamma(4-a-b-(2+i+j)*x));
82 }
83
84 /* Abstract base class (ABC) for all types of vertices. */
85 class vertex {
86 public:
87         vertex(ijpair ij = ijpair(0,0)) : indices(ij) { }
88         void increment_indices(const ijpair &ind) { indices.first += ind.first; indices.second += ind.second; }
89         virtual ~vertex() { }
90         virtual vertex* copy() const = 0;
91         virtual ijpair get_increment() const { return indices; }
92         virtual const ex evaluate(const symbol &x, const unsigned grad) const = 0;
93         bool operator==(const vertex &v) const { return (indices==v.indices); }
94         bool operator<(const vertex &v) const { return (indices<v.indices); }
95 protected:
96         ijpair indices;
97 };
98
99
100 /*
101  * Class of vertices of type Sigma.
102  */
103 class Sigma : public vertex {
104 public:
105         Sigma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
106         vertex* copy() const { return new Sigma(*this); }
107         ijpair get_increment() const { return ijpair(indices.first+indices.second+1, 0); }
108         const ex evaluate(const symbol &x, const unsigned grad) const;
109 private:
110 };
111
112 const ex Sigma::evaluate(const symbol &x, const unsigned grad) const
113 {
114         // c.f. comments in node::evaluate()
115         static map<Sigma,ex> catalog;
116         static unsigned prev_grad = 0;
117         if (grad>prev_grad) {
118                 catalog.clear();
119                 prev_grad = grad;
120         }
121         map<Sigma,ex>::iterator pos = catalog.find(*this);
122         if (pos==catalog.end()) {
123                 int i = indices.first;
124                 int j = indices.second;
125                 const ex result = ((F_ab(0,i,1,j,x)+F_ab(1,i,1,j,x)-F_ab(1,i,0,j,x))/2).series(x==0, grad).expand();
126                 pos = catalog.insert(map<Sigma,ex>::value_type(*this,result)).first;
127                 if (grad<prev_grad)
128                         prev_grad = grad;
129         }
130         return pos->second;
131 }
132
133
134 /** Class of vertices of type Sigma_flipped, sitting in the upper fermionline of Vacuum; no consequences for Gamma. */
135 class Sigma_flipped : public Sigma {
136 public:
137         Sigma_flipped(ijpair ij = ijpair(0,0)) : Sigma(ij) { }
138         vertex* copy() const { return new Sigma_flipped(*this); }
139         ijpair get_increment() const { return ijpair(0, indices.first+indices.second+1); }
140         const ex evaluate(const symbol &x, const unsigned grad) const { return Sigma::evaluate(x, grad); }
141 private:
142 };
143
144
145 /*
146  *Class of vertices of type Gamma.
147  */
148 class Gamma : public vertex {
149 public:
150         Gamma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
151         vertex* copy() const { return new Gamma(*this); }
152         ijpair get_increment() const { return ijpair(indices.first+indices.second+1, 0); }
153         const ex evaluate(const symbol &x, const unsigned grad) const;
154 private:
155 };
156
157 const ex Gamma::evaluate(const symbol &x, const unsigned grad) const
158 {
159         // c.f. comments in node::evaluate()
160         static map<Gamma,ex> catalog;
161         static unsigned prev_grad = 0;
162         if (prev_grad<grad) {
163                 catalog.clear();
164                 prev_grad = grad;
165         }
166         map<Gamma,ex>::iterator pos = catalog.find(*this);
167         if (pos==catalog.end()) {
168                 int i = indices.first;
169                 int j = indices.second;
170                 const ex result = (F_ab(1,i,1,j,x)).series(x==0,grad).expand();
171                 pos = catalog.insert(map<Gamma,ex>::value_type(*this,result)).first;
172                 if (grad<prev_grad)
173                         prev_grad = grad;
174         }
175         return pos->second;
176 }
177
178
179 /*
180  * Class of vertices of type Vacuum.
181  */
182 class Vacuum : public vertex {
183 public:
184         Vacuum(ijpair ij = ijpair(0,0)) : vertex(ij) { }
185         vertex* copy() const { return new Vacuum(*this); }
186         ijpair get_increment() const { return ijpair(0, indices.first+indices.second+1); }
187         const ex evaluate(const symbol &x, const unsigned grad) const;
188 private:
189 };
190
191 const ex Vacuum::evaluate(const symbol &x, const unsigned grad) const
192 {
193         // c.f. comments in node::evaluate()
194         static map<Vacuum,ex> catalog;
195         static unsigned prev_grad = 0;
196         if (prev_grad<grad) {
197                 catalog.clear();
198                 prev_grad = grad;
199         }
200         map<Vacuum,ex>::iterator pos = catalog.find(*this);
201         if (pos==catalog.end()) {
202                 int i = indices.first;
203                 int j = indices.second;
204                 const ex result = ((-TrOne*(F_ab(0,i,1,j,x)-F_ab(1,i,1,j,x)+F_ab(1,i,0,j,x)))/2).series(x==0,grad).expand();
205                 pos = catalog.insert(map<Vacuum,ex>::value_type(*this,result)).first;
206                 if (grad<prev_grad)
207                         prev_grad = grad;
208         }
209         return pos->second;
210 }
211
212
213 /*
214  * Class of nodes (or trees or subtrees), including list of children.
215  */
216 class node {
217 public:
218         node(const vertex &v) { vert = v.copy(); }
219         node(const node &n) { vert = (n.vert)->copy(); children = n.children; }
220         const node & operator=(const node &);
221         ~node() { delete vert; }
222         void add_child(const node &, bool = false);
223         const ex evaluate(const symbol &x, unsigned grad) const;
224         unsigned total_edges() const;
225         bool operator==(const node &) const;
226         bool operator<(const node &) const;
227 private:
228         vertex *vert;
229         multiset<child> children;
230 };
231
232 const node & node::operator=(const node &n)
233 {
234         if (this!=&n) {
235                 delete vert;
236                 vert = (n.vert)->copy();
237                 children = n.children;
238         }
239         return *this;
240 }
241
242 void node::add_child(const node &childnode, bool cut)
243 {
244         children.insert(child(childnode, cut));
245         if(!cut)
246                 vert->increment_indices(childnode.vert->get_increment());
247 }
248
249 const ex node::evaluate(const symbol &x, unsigned grad) const
250 {
251         static map<node,ex> catalog;    // lookup table for already evaluated subnodes
252         static unsigned prev_grad = 0;  // grad argument that the catalog has been build for
253         if (grad>prev_grad) {
254                 // We haven't computed far enough last time.  Our catalog cannot cope with
255                 // the demands for this value of grad so let's clear it.
256                 catalog.clear();
257                 prev_grad = grad;
258         }
259         ex product = 1;   // accumulator for all the children
260         for (multiset<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
261                 map<node,ex>::iterator pos = catalog.find(i->first);
262                 if (pos==catalog.end()) {
263                         pos = catalog.insert(map<node,ex>::value_type(i->first,i->first.evaluate(x,grad).series(x==0,grad).expand())).first;
264                         if (grad<prev_grad) {
265                                 // We have just spoiled the catalog by inserting a series computed
266                                 // to lower grad than the others in it.  So let's make sure next time
267                                 // we don't use one of the newly inserted ones by bumping prev_grad
268                                 // down to the current value of grad.
269                                 prev_grad = grad;
270                         }
271                 }
272                 if (!i->second)
273                         product *= pos->second;
274                 else
275                         product *= -div_part(pos->second,x,grad);
276         }
277         return (product * vert->evaluate(x,grad));
278 }
279
280 unsigned node::total_edges() const
281 {
282         unsigned accu = 0;
283         for (multiset<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
284                 accu += i->first.total_edges();
285                 ++accu;
286         }
287         return accu;
288 }
289
290 bool node::operator==(const node &n) const
291 {
292         // Are the types of the top-level node vertices equal?
293         if (typeid(*vert)!=typeid(*n.vert))
294                 return false;
295         // Are the indices of the top-level nodes equal?
296         if (!(*vert==*n.vert))
297                 return false;
298         // Are the sets of children equal, one by one?
299         return (children==n.children);
300 }
301
302 bool node::operator<(const node &n) const
303 {
304         // Are the types of the top-level node vertices different?
305         if (typeid(*vert)!=typeid(*n.vert))
306                 return typeid(*vert).before(typeid(*n.vert));
307         // Are the indices of the top-level nodes different?
308         if (!(*vert==*n.vert))
309                 return (vert<n.vert);
310         // Are the sets of children different, one by one?
311         return (children<n.children);
312 }
313
314 /*
315  * These operators let us write down trees in an intuitive way, by adding
316  * arbitrarily complex children to a given vertex.  The eye candy that can be
317  * produced with it makes detection of errors much simpler than with code
318  * written using calls to node's method add_child() because it allows for
319  * editor-assisted indentation.
320  */
321 const node operator+(const node &n, const child &c)
322 {
323         node nn(n);
324         nn.add_child(c.first, c.second);
325         return nn;
326 }
327
328 void operator+=(node &n, const child &c)
329 {
330         n.add_child(c.first, c.second);
331 }
332
333
334 /*                Gamma
335  *                  |
336  *                Gamma
337  */
338 static const node tree1(unsigned cuts=0)
339 {
340         return (Gamma()
341                 + child(Gamma(),
342                         bool(cuts & 1)));
343 }
344
345 /*                Gamma
346  *               /  |  \
347  *        Vacuum  Gamma  Vacuum
348  *       /   |  \
349  *  Sigma Sigma Sigma0
350  */
351 static const node tree2(unsigned cuts=0)
352 {
353         return (Gamma()
354                 + child(Vacuum()
355                         + child(Sigma(), bool(cuts & 1))
356                         + child(Sigma(), bool(cuts & 2))
357                         + child(Sigma_flipped(), bool(cuts & 4)),
358                         bool(cuts & 8))
359                 + child(Gamma(), bool(cuts & 16))
360                 + child(Gamma(), bool(cuts & 32)));
361 }
362
363 /*                Gamma
364  *                  |
365  *                Gamma
366  *                  |
367  *                Gamma
368  *                /   \
369  *           Vacuum    Gamma
370  *           /    \       \
371  *        Sigma  Sigma   Sigma
372  */
373 static const node tree3(unsigned cuts=0)
374 {
375         return (Gamma()
376                 + child(Gamma()
377                         + child(Gamma()
378                                 + child(Gamma()
379                                         + child(Sigma(), bool(cuts & 1)),
380                                         bool(cuts & 2))
381                                 + child(Vacuum()
382                                         + child(Sigma(), bool(cuts & 4))
383                                         + child(Sigma(), bool(cuts & 8)),
384                                 bool(cuts & 16)),
385                         bool(cuts & 32)),
386                 bool(cuts & 64)));
387 }
388
389 /*                Gamma
390  *                /   \
391  *           Sigma     Vacuum
392  *          /   \       /   \
393  *      Sigma Sigma  Sigma0 Sigma
394  */
395 static const node tree4(unsigned cuts=0)
396 {
397         return (Gamma()
398                 + child(Sigma()
399                         + child(Sigma(), bool(cuts & 1))
400                         + child(Sigma(), bool(cuts & 2)),
401                         bool(cuts & 4))
402                 + child(Vacuum()
403                         + child(Sigma_flipped(), bool(cuts & 8))
404                         + child(Sigma(), bool(cuts & 16)),
405                         bool(cuts & 32)));
406 }
407
408 /*                Sigma
409  *               /  |  \
410  *         Sigma Vacuum  Vacuum
411  *               /    \       \
412  *            Sigma  Sigma0   Sigma
413  */
414 static const node tree5(unsigned cuts=0)
415 {
416         return (Sigma()
417                 + child(Sigma(), bool(cuts & 1))
418                 + child(Vacuum()
419                         + child(Sigma(), bool(cuts & 2))
420                         + child(Sigma_flipped(), bool(cuts & 4)),
421                         bool(cuts & 8))
422                 + child(Vacuum()
423                         + child(Sigma(), bool(cuts & 16)),
424                         bool(cuts & 32)));
425 }
426
427 /*               Vacuum
428  *               /    \
429  *           Sigma    Sigma0
430  *             |        |
431  *           Sigma    Sigma
432  *                      |
433  *                    Vacuum
434  */
435 static const node tree6(unsigned cuts=0)
436 {
437         return (Vacuum()
438                 + child(Sigma()
439                         + child(Sigma(), bool(cuts & 1)),
440                         bool(cuts & 2))
441                 + child(Sigma_flipped()
442                         + child(Sigma()
443                                 + child(Vacuum(), bool(cuts & 4)),
444                                 bool(cuts & 8)),
445                         bool(cuts & 16)));
446 }
447
448 static unsigned test_tree(const node tree_generator(unsigned))
449 {
450         const int edges = tree_generator(0).total_edges();
451         const int vertices = edges+1;
452         
453         // fill a vector of all possible 2^edges combinations of cuts...
454         vector<node> counter;
455         for (unsigned i=0; i<(1U<<edges); ++i)
456                 counter.push_back(tree_generator(i));
457         
458         // ...the sum, when evaluated and reexpanded, is the antipode...
459         ex result = 0;
460         for (vector<node>::iterator i=counter.begin(); i!=counter.end(); ++i)
461                 result = (result+i->evaluate(x,vertices-1)).series(x==0,vertices-1).expand();
462         
463         // ...and has the nice property that in each term all the Eulers cancel:
464         if (result.has(Euler)) {
465                 clog << "The antipode was miscalculated\nAntipode==" << result
466                      << "\nshould not have any occurrence of Euler" << endl;
467                 return 1;
468         } else if (result.ldegree(x)!=-vertices || result.degree(x)!=0) {
469                 clog << "The antipode was miscalculated\nAntipode==" << result
470                      << "\nshould run from " << x << "^(" << -vertices << ") to "
471                      << x << "^(0)" << "but it runs from " << x << "^("
472                      << result.ldegree(x) << ")" << "to " << x << "^("
473                      << result.degree(x) << ")" << endl;
474                 return 1;
475         }
476         return 0;
477 }
478
479 unsigned time_antipode()
480 {
481         unsigned result = 0;
482         timer jaeger_le_coultre;
483         
484         cout << "timing computation of antipodes in Yukawa theory" << flush;
485         
486         if (do_test) {
487                 jaeger_le_coultre.start();
488                 result += test_tree(tree1);  cout << '.' << flush;
489                 result += test_tree(tree2);  cout << '.' << flush;
490                 result += test_tree(tree3);  cout << '.' << flush;
491                 result += test_tree(tree4);  cout << '.' << flush;
492                 result += test_tree(tree5);  cout << '.' << flush;
493                 result += test_tree(tree6);  cout << '.' << flush;
494                 
495                 cout << jaeger_le_coultre.read() << "s (total)" << endl;
496         } else {
497                 cout << " disabled" << endl;
498         }
499         return result;
500 }
501
502 extern void randomify_symbol_serials();
503
504 int main(int argc, char** argv)
505 {
506         randomify_symbol_serials();
507         cout << setprecision(2) << showpoint;
508         return time_antipode();
509 }