d716ea24117031c417fbb6ae0a294255ad3ebdb0
[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 Isabella Bierenbaum and Dirk Kreimer.
13  *  For details, please see the diploma theses of Isabella Bierenbaum.
14  */
15
16 /*
17  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
18  *
19  *  This program is free software; you can redistribute it and/or modify
20  *  it under the terms of the GNU General Public License as published by
21  *  the Free Software Foundation; either version 2 of the License, or
22  *  (at your option) any later version.
23  *
24  *  This program is distributed in the hope that it will be useful,
25  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
26  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27  *  GNU General Public License for more details.
28  *
29  *  You should have received a copy of the GNU General Public License
30  *  along with this program; if not, write to the Free Software
31  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
32  */
33
34 #include "times.h"
35 #include <utility>
36 #include <list>
37
38 // whether to run this beast or not:
39 static const bool do_test = true;
40
41 typedef pair<unsigned, unsigned> ijpair;
42 typedef pair<class node, bool> child;
43
44 const constant TrOne("Tr[One]", numeric(4));
45
46 /* Extract only the divergent part of a series and discard the rest. */
47 static ex div_part(const ex &exarg, const symbol &x, unsigned grad)
48 {
49         unsigned order = grad;
50         ex exser;
51         // maybe we have to generate more terms on the series (obnoxious):
52         do {
53                 exser = exarg.series(x==0, order);
54                 ++order;
55         } while (exser.degree(x) < 0);
56         ex exser_trunc;
57         for (int i=exser.ldegree(x); i<0; ++i)
58                 exser_trunc += exser.coeff(x,i)*pow(x,i);
59         // NB: exser_trunc is by construction collected in x.
60         return exser_trunc;
61 }
62
63 /* F_ab(a, i, b, j, "x") is a common pattern in all vertex evaluators. */
64 static ex F_ab(int a, int i, int b, int j, const symbol &x)
65 {
66         if ((i==0 && a<=0) || (j==0 && b<=0))
67                 return 0;
68         else
69                 return (tgamma(2-a-(i+1)*x)*
70                         tgamma(2-b-(1+j)*x)*
71                         tgamma(a+b-2+(1+i+j)*x)/
72                         tgamma(a+i*x)/
73                         tgamma(b+j*x)/tgamma(4-a-b-(2+i+j)*x));
74 }
75
76 /* Abstract base class (ABC) for all types of vertices. */
77 class vertex {
78 public:
79         vertex(ijpair ij = ijpair(0,0)) : indices(ij) { }
80         void increment_indices(const ijpair &ind) { indices.first += ind.first; indices.second += ind.second; }
81         virtual ~vertex() { }
82         virtual vertex* copy(void) const = 0;
83         virtual ijpair get_increment(void) const { return indices; }
84         virtual ex evaluate(const symbol &x) const = 0; 
85 protected:
86         ijpair indices;
87 };
88
89
90 /*
91  * Class of vertices of type Sigma.
92  */
93 class Sigma : public vertex {
94 public:
95         Sigma(ijpair ij = ijpair(0,0), bool f = true) : vertex(ij), flag(f) { }
96         vertex* copy(void) const { return new Sigma(*this); }
97         ijpair get_increment(void) const;
98         ex evaluate(const symbol &x) const;
99 private:
100         bool flag;
101 };
102
103 ijpair Sigma::get_increment(void) const
104 {
105         if (flag == true)
106             return  ijpair(indices.first+1, indices.second);
107         else
108             return  ijpair(indices.first, indices.second+1);
109 }
110
111 ex Sigma::evaluate(const symbol &x) const
112 {
113         int i = indices.first;
114         int j = indices.second;
115         return (F_ab(0,i,1,j,x)+F_ab(1,i,1,j,x)-F_ab(1,i,0,j,x))/2;
116 }
117
118
119 /*
120  *Class of vertices of type Gamma.
121  */
122 class Gamma : public vertex {
123 public:
124         Gamma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
125         vertex* copy(void) const { return new Gamma(*this); }
126         ijpair get_increment(void) const { return ijpair(indices.first+indices.second+1, 0); }
127         ex evaluate(const symbol &x) const;
128 private:
129 };
130
131 ex Gamma::evaluate(const symbol &x) const
132 {
133         int i = indices.first;
134         int j = indices.second;
135         return F_ab(1,i,1,j,x);
136 }
137
138
139 /*
140  * Class of vertices of type Vacuum.
141  */
142 class Vacuum : public vertex {
143 public:
144         Vacuum(ijpair ij = ijpair(0,0)) : vertex(ij) { }
145         vertex* copy(void) const { return new Vacuum(*this); }
146         ijpair get_increment() const { return ijpair(0, indices.first+indices.second+1); }
147         ex evaluate(const symbol &x) const;
148 private:
149 };
150
151 ex Vacuum::evaluate(const symbol &x) const
152 {
153         int i = indices.first;
154         int j = indices.second;
155         return (-TrOne*(F_ab(0,i,1,j,x)-F_ab(1,i,1,j,x)+F_ab(1,i,0,j,x)))/2;
156 }
157
158
159 /*
160  * Class of nodes (or trees or subtrees), including list of children.
161  */
162 class node {
163 public:
164         node(const vertex &v) { vert = v.copy(); }
165         node(const node &n) { vert = (n.vert)->copy(); children = n.children; }
166         const node & operator=(const node &);
167         ~node() { delete vert; }
168         void add_child(const node &, bool = false);
169         ex evaluate(const symbol &x, unsigned grad) const;
170         unsigned total_edges(void) const;
171 private:
172         vertex *vert;
173         list<child> children;
174 };
175
176 const node & node::operator=(const node &n)
177 {
178         delete vert;
179         vert = (n.vert)->copy();
180         children = n.children;
181         return *this;
182 }
183
184 void node::add_child(const node &childnode, bool cut)
185 {
186         children.push_back(child(childnode, cut));
187         if(!cut)
188                 vert->increment_indices(childnode.vert->get_increment());
189 }
190
191 ex node::evaluate(const symbol &x, unsigned grad) const
192 {
193         ex product = 1;
194         for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
195                 if (!i->second)
196                         product *= i->first.evaluate(x,grad);
197                 else
198                         product *= -div_part(i->first.evaluate(x,grad),x,grad);
199         }
200         return (product * vert->evaluate(x));
201 }
202
203 unsigned node::total_edges(void) const
204 {
205         unsigned accu = 0;
206         for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
207                 accu += i->first.total_edges();
208                 ++accu;
209         }
210         return accu;
211 }
212
213
214 /*
215  * These operators let us write down trees in an intuitive way, by adding
216  * arbitrarily complex children to a given vertex.  The eye candy that can be
217  * produced with it makes detection of errors much simpler than with code
218  * written using calls to node's method add_child() because it allows for
219  * editor-assisted indentation.
220  */
221 const node operator+(const node &n, const child &c)
222 {
223         node nn(n);
224         nn.add_child(c.first, c.second);
225         return nn;
226 }
227 void operator+=(node &n, const child &c)
228 {
229         n.add_child(c.first, c.second);
230 }
231
232 /*
233  * Build this sample rooted tree characterized by a certain combination of
234  * cut or uncut edges as specified by the unsigned parameter:
235  *              Gamma
236  *              /   \
237  *         Sigma     Vacuum
238  *        /   \       /   \
239  *    Sigma Sigma  Sigma0 Sigma
240  */
241 static node mytree(unsigned cuts=0)
242 {
243         return (Gamma()
244                 + child(Sigma()
245                         + child(Sigma(), bool(cuts & 1))
246                         + child(Sigma(), bool(cuts & 2)),
247                         bool(cuts & 4))
248                 + child(Vacuum()
249                         + child(Sigma(ijpair(0,0),false), bool(cuts & 8))
250                         + child(Sigma(), bool(cuts & 16)),
251                         bool(cuts & 32)));
252 }
253
254
255 static unsigned test(void)
256 {
257         const symbol x("x");
258         
259         const unsigned edges = mytree().total_edges();
260         const unsigned vertices = edges+1;
261         
262         // fill a vector of all possible 2^edges combinations of cuts...
263         vector<node> counter;
264         for (unsigned i=0; i<(1U<<edges); ++i)
265                 counter.push_back(mytree(i));
266         
267         // ...the sum, when evaluated, is the antipode...
268         ex accu = 0;
269         for (vector<node>::iterator i=counter.begin(); i!=counter.end(); ++i)
270                 accu += i->evaluate(x,vertices);
271         
272         // ...which is only interesting term-wise in the series expansion...
273         ex result = accu.series(x==0,vertices).expand().normal();
274         
275         // ...and has the nice property that in each term all the Eulers cancel:
276         if (result.has(Euler)) {
277                 clog << "The antipode was miscalculated\nAntipode==" << result
278                      << "\nshould not have any occurrence of Euler" << endl;
279                 return 1;
280         }
281         return 0;
282 }
283
284 unsigned time_antipode(void)
285 {
286         unsigned result = 0;
287         unsigned count = 0;
288         timer jaeger_le_coultre;
289         double time = .0;
290         
291         cout << "timing computation of an antipode in Yukawa theory" << flush;
292         clog << "-------computation of an antipode in Yukawa theory" << endl;
293         
294         if (do_test) {
295                 jaeger_le_coultre.start();
296                 // correct for very small times:
297                 do {
298                         result = test();
299                         ++count;
300                 } while ((time=jaeger_le_coultre.read())<0.1 && !result);
301                 cout << '.' << flush;
302                 
303                 if (!result) {
304                         cout << " passed ";
305                         clog << "(no output)" << endl;
306                 } else {
307                         cout << " failed ";
308                 }
309                 cout << int(1000*(time/count))*0.001 << 's' << endl;
310         } else {
311                 cout << " disabled" << endl;
312                 clog << "(no output)" << endl;
313         }
314         
315         return result;
316 }