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