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