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