1 /** @file time_antipode.cpp
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.
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.
21 * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
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.
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.
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
43 // whether to run this beast or not:
44 static const bool do_test = true;
46 typedef pair<unsigned, unsigned> ijpair;
47 typedef pair<class node, bool> child;
49 const constant TrOne("Tr[One]", numeric(4));
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)
54 const ex exser = exarg.series(x==0, grad);
55 if (exser.degree(x)<0)
56 throw runtime_error("divergent part truncation disaster");
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.
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)
67 if ((i==0 && a<=0) || (j==0 && b<=0))
70 return (tgamma(2-a-(i+1)*x)*
72 tgamma(a+b-2+(1+i+j)*x)/
74 tgamma(b+j*x)/tgamma(4-a-b-(2+i+j)*x));
77 /* Abstract base class (ABC) for all types of vertices. */
80 vertex(ijpair ij = ijpair(0,0)) : indices(ij) { }
81 void increment_indices(const ijpair &ind) { indices.first += ind.first; indices.second += ind.second; }
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;
92 * Class of vertices of type Sigma.
94 class Sigma : public vertex {
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;
104 ijpair Sigma::get_increment(void) const
107 return ijpair(indices.first+1, indices.second);
109 return ijpair(indices.first, indices.second+1);
112 ex Sigma::evaluate(const symbol &x) const
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;
121 *Class of vertices of type Gamma.
123 class Gamma : public vertex {
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;
132 ex Gamma::evaluate(const symbol &x) const
134 int i = indices.first;
135 int j = indices.second;
136 return F_ab(1,i,1,j,x);
141 * Class of vertices of type Vacuum.
143 class Vacuum : public vertex {
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;
152 ex Vacuum::evaluate(const symbol &x) const
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;
161 * Class of nodes (or trees or subtrees), including list of children.
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;
174 vector<child> children;
177 const node & node::operator=(const node &n)
180 vert = (n.vert)->copy();
181 children = n.children;
185 void node::add_child(const node &childnode, bool cut)
187 children.push_back(child(childnode, cut));
189 vert->increment_indices(childnode.vert->get_increment());
192 ex node::evaluate(const symbol &x, unsigned grad) const
195 for (vector<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
197 product *= i->first.evaluate(x,grad);
199 product *= -div_part(i->first.evaluate(x,grad),x,grad);
201 return (product * vert->evaluate(x));
204 unsigned node::total_edges(void) const
207 for (vector<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
208 accu += i->first.total_edges();
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.
222 const node operator+(const node &n, const child &c)
225 nn.add_child(c.first, c.second);
228 void operator+=(node &n, const child &c)
230 n.add_child(c.first, c.second);
234 * Build this sample rooted tree characterized by a certain combination of
235 * cut or uncut edges as specified by the unsigned parameter:
240 * Sigma Sigma Sigma0 Sigma
242 static node mytree(unsigned cuts=0)
246 + child(Sigma(), bool(cuts & 1))
247 + child(Sigma(), bool(cuts & 2)),
250 + child(Sigma(ijpair(0,0),false), bool(cuts & 8))
251 + child(Sigma(), bool(cuts & 16)),
256 static unsigned test(void)
260 const unsigned edges = mytree().total_edges();
261 const unsigned vertices = edges+1;
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));
268 // ...the sum, when evaluated, is the antipode...
270 for (vector<node>::iterator i=counter.begin(); i!=counter.end(); ++i)
271 accu += i->evaluate(x,vertices);
273 // ...which is only interesting term-wise in the series expansion...
274 ex result = accu.series(x==0,vertices).expand().normal();
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;
285 unsigned time_antipode(void)
289 timer jaeger_le_coultre;
292 cout << "timing computation of an antipode in Yukawa theory" << flush;
293 clog << "-------computation of an antipode in Yukawa theory" << endl;
296 jaeger_le_coultre.start();
297 // correct for very small times:
301 } while ((time=jaeger_le_coultre.read())<0.1 && !result);
302 cout << '.' << flush;
306 clog << "(no output)" << endl;
310 cout << int(1000*(time/count))*0.001 << 's' << endl;
312 cout << " disabled" << endl;
313 clog << "(no output)" << endl;