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