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 Isabella Bierenbaum and Dirk Kreimer.
13 * For details, please see the diploma theses of Isabella Bierenbaum.
17 * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
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.
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.
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
38 // whether to run this beast or not:
39 static const bool do_test = true;
41 typedef pair<unsigned, unsigned> ijpair;
42 typedef pair<class node, bool> child;
44 const constant TrOne("Tr[One]", numeric(4));
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)
49 unsigned order = grad;
51 // maybe we have to generate more terms on the series (obnoxious):
53 exser = exarg.series(x==0, order);
55 } while (exser.degree(x) < 0);
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.
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)
66 if ((i==0 && a<=0) || (j==0 && b<=0))
69 return (tgamma(2-a-(i+1)*x)*
71 tgamma(a+b-2+(1+i+j)*x)/
73 tgamma(b+j*x)/tgamma(4-a-b-(2+i+j)*x));
76 /* Abstract base class (ABC) for all types of vertices. */
79 vertex(ijpair ij = ijpair(0,0)) : indices(ij) { }
80 void increment_indices(const ijpair &ind) { indices.first += ind.first; indices.second += ind.second; }
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;
91 * Class of vertices of type Sigma.
93 class Sigma : public vertex {
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;
103 ijpair Sigma::get_increment(void) const
106 return ijpair(indices.first+1, indices.second);
108 return ijpair(indices.first, indices.second+1);
111 ex Sigma::evaluate(const symbol &x) const
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;
120 *Class of vertices of type Gamma.
122 class Gamma : public vertex {
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;
131 ex Gamma::evaluate(const symbol &x) const
133 int i = indices.first;
134 int j = indices.second;
135 return F_ab(1,i,1,j,x);
140 * Class of vertices of type Vacuum.
142 class Vacuum : public vertex {
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;
151 ex Vacuum::evaluate(const symbol &x) const
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;
160 * Class of nodes (or trees or subtrees), including list of children.
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;
173 list<child> children;
176 const node & node::operator=(const node &n)
179 vert = (n.vert)->copy();
180 children = n.children;
184 void node::add_child(const node &childnode, bool cut)
186 children.push_back(child(childnode, cut));
188 vert->increment_indices(childnode.vert->get_increment());
191 ex node::evaluate(const symbol &x, unsigned grad) const
194 for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
196 product *= i->first.evaluate(x,grad);
198 product *= -div_part(i->first.evaluate(x,grad),x,grad);
200 return (product * vert->evaluate(x));
203 unsigned node::total_edges(void) const
206 for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
207 accu += i->first.total_edges();
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.
221 const node operator+(const node &n, const child &c)
224 nn.add_child(c.first, c.second);
227 void operator+=(node &n, const child &c)
229 n.add_child(c.first, c.second);
233 * Build this sample rooted tree characterized by a certain combination of
234 * cut or uncut edges as specified by the unsigned parameter:
239 * Sigma Sigma Sigma0 Sigma
241 static node mytree(unsigned cuts=0)
245 + child(Sigma(), bool(cuts & 1))
246 + child(Sigma(), bool(cuts & 2)),
249 + child(Sigma(ijpair(0,0),false), bool(cuts & 8))
250 + child(Sigma(), bool(cuts & 16)),
255 static unsigned test(void)
259 const unsigned edges = mytree().total_edges();
260 const unsigned vertices = edges+1;
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));
267 // ...the sum, when evaluated, is the antipode...
269 for (vector<node>::iterator i=counter.begin(); i!=counter.end(); ++i)
270 accu += i->evaluate(x,vertices);
272 // ...which is only interesting term-wise in the series expansion...
273 ex result = accu.series(x==0,vertices).expand().normal();
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;
284 unsigned time_antipode(void)
288 timer jaeger_le_coultre;
291 cout << "timing computation of an antipode in Yukawa theory" << flush;
292 clog << "-------computation of an antipode in Yukawa theory" << endl;
295 jaeger_le_coultre.start();
296 // correct for very small times:
300 } while ((time=jaeger_le_coultre.read())<0.1 && !result);
301 cout << '.' << flush;
305 clog << "(no output)" << endl;
309 cout << int(1000*(time/count))*0.001 << 's' << endl;
311 cout << " disabled" << endl;
312 clog << "(no output)" << endl;