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 see <http://www.arXiv.org/abs/hep-th/0111192>.
19 * GiNaC Copyright (C) 1999-2019 Johannes Gutenberg University Mainz, Germany
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.
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.
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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
38 using namespace GiNaC;
48 // whether to run this beast or not:
49 constexpr bool do_test = true;
51 // regularization parameter:
52 static const symbol x("x");
54 typedef pair<unsigned, unsigned> ijpair;
55 typedef pair<class node, bool> child;
57 const constant TrOne("Tr[One]", numeric(4));
59 /* Extract only the divergent part of a series and discard the rest. */
60 static ex div_part(const ex &exarg, const symbol &x, unsigned grad)
62 const ex exser = exarg.series(x==0, grad);
63 if (exser.degree(x)<0)
64 throw runtime_error("divergent part truncation disaster");
66 for (int i=exser.ldegree(x); i<0; ++i)
67 exser_trunc += exser.coeff(x,i)*pow(x,i);
68 // NB: exser_trunc is by construction collected in x.
72 /* F_ab(a, i, b, j, "x") is a common pattern in all vertex evaluators. */
73 static ex F_ab(int a, int i, int b, int j, const symbol &x)
76 if ((i==0 && a<=0) || (j==0 && b<=0))
79 return (tgamma(2-a-(i+1)*x)*
81 tgamma(a+b-2+(1+i+j)*x)/
83 tgamma(b+j*x)/tgamma(4-a-b-(2+i+j)*x));
86 /* Abstract base class (ABC) for all types of vertices. */
89 vertex(ijpair ij = ijpair(0,0)) : indices(ij) { }
90 void increment_indices(const ijpair &ind) { indices.first += ind.first; indices.second += ind.second; }
92 virtual vertex* copy() const = 0;
93 virtual ijpair get_increment() const { return indices; }
94 virtual const ex evaluate(const symbol &x, const unsigned grad) const = 0;
95 bool operator==(const vertex &v) const { return (indices==v.indices); }
96 bool operator<(const vertex &v) const { return (indices<v.indices); }
103 * Class of vertices of type Sigma.
105 class Sigma : public vertex {
107 Sigma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
108 vertex* copy() const override { return new Sigma(*this); }
109 ijpair get_increment() const override { return ijpair(indices.first+indices.second+1, 0); }
110 const ex evaluate(const symbol &x, const unsigned grad) const override;
114 const ex Sigma::evaluate(const symbol &x, const unsigned grad) const
116 // c.f. comments in node::evaluate()
117 static map<Sigma,ex> catalog;
118 static unsigned prev_grad = 0;
119 if (grad>prev_grad) {
123 map<Sigma,ex>::iterator pos = catalog.find(*this);
124 if (pos==catalog.end()) {
125 int i = indices.first;
126 int j = indices.second;
127 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();
128 pos = catalog.insert(map<Sigma,ex>::value_type(*this,result)).first;
136 /** Class of vertices of type Sigma_flipped, sitting in the upper fermion line of Vacuum; no consequences for Gamma. */
137 class Sigma_flipped : public Sigma {
139 Sigma_flipped(ijpair ij = ijpair(0,0)) : Sigma(ij) { }
140 vertex* copy() const override { return new Sigma_flipped(*this); }
141 ijpair get_increment() const override { return ijpair(0, indices.first+indices.second+1); }
142 const ex evaluate(const symbol &x, const unsigned grad) const override { return Sigma::evaluate(x, grad); }
148 *Class of vertices of type Gamma.
150 class Gamma : public vertex {
152 Gamma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
153 vertex* copy() const override { return new Gamma(*this); }
154 ijpair get_increment() const override { return ijpair(indices.first+indices.second+1, 0); }
155 const ex evaluate(const symbol &x, const unsigned grad) const override;
159 const ex Gamma::evaluate(const symbol &x, const unsigned grad) const
161 // c.f. comments in node::evaluate()
162 static map<Gamma,ex> catalog;
163 static unsigned prev_grad = 0;
164 if (prev_grad<grad) {
168 map<Gamma,ex>::iterator pos = catalog.find(*this);
169 if (pos==catalog.end()) {
170 int i = indices.first;
171 int j = indices.second;
172 const ex result = (F_ab(1,i,1,j,x)).series(x==0,grad).expand();
173 pos = catalog.insert(map<Gamma,ex>::value_type(*this,result)).first;
182 * Class of vertices of type Vacuum.
184 class Vacuum : public vertex {
186 Vacuum(ijpair ij = ijpair(0,0)) : vertex(ij) { }
187 vertex* copy() const override { return new Vacuum(*this); }
188 ijpair get_increment() const override { return ijpair(0, indices.first+indices.second+1); }
189 const ex evaluate(const symbol &x, const unsigned grad) const override;
193 const ex Vacuum::evaluate(const symbol &x, const unsigned grad) const
195 // c.f. comments in node::evaluate()
196 static map<Vacuum,ex> catalog;
197 static unsigned prev_grad = 0;
198 if (prev_grad<grad) {
202 map<Vacuum,ex>::iterator pos = catalog.find(*this);
203 if (pos==catalog.end()) {
204 int i = indices.first;
205 int j = indices.second;
206 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();
207 pos = catalog.insert(map<Vacuum,ex>::value_type(*this,result)).first;
216 * Class of nodes (or trees or subtrees), including list of children.
220 node(const vertex &v) { vert = v.copy(); }
221 node(const node &n) { vert = (n.vert)->copy(); children = n.children; }
222 const node & operator=(const node &);
223 ~node() { delete vert; }
224 void add_child(const node &, bool = false);
225 const ex evaluate(const symbol &x, unsigned grad) const;
226 unsigned total_edges() const;
227 bool operator==(const node &) const;
228 bool operator<(const node &) const;
231 multiset<child> children;
234 const node & node::operator=(const node &n)
238 vert = (n.vert)->copy();
239 children = n.children;
244 void node::add_child(const node &childnode, bool cut)
246 children.insert(child(childnode, cut));
248 vert->increment_indices(childnode.vert->get_increment());
251 const ex node::evaluate(const symbol &x, unsigned grad) const
253 static map<node,ex> catalog; // lookup table for already evaluated subnodes
254 static unsigned prev_grad = 0; // grad argument that the catalog has been build for
255 if (grad>prev_grad) {
256 // We haven't computed far enough last time. Our catalog cannot cope with
257 // the demands for this value of grad so let's clear it.
261 ex product = 1; // accumulator for all the children
262 for (multiset<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
263 map<node,ex>::iterator pos = catalog.find(i->first);
264 if (pos==catalog.end()) {
265 pos = catalog.insert(map<node,ex>::value_type(i->first,i->first.evaluate(x,grad).series(x==0,grad).expand())).first;
266 if (grad<prev_grad) {
267 // We have just spoiled the catalog by inserting a series computed
268 // to lower grad than the others in it. So let's make sure next time
269 // we don't use one of the newly inserted ones by bumping prev_grad
270 // down to the current value of grad.
275 product *= pos->second;
277 product *= -div_part(pos->second,x,grad);
279 return (product * vert->evaluate(x,grad));
282 unsigned node::total_edges() const
285 for (multiset<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
286 accu += i->first.total_edges();
292 bool node::operator==(const node &n) const
294 // Are the types of the top-level node vertices equal?
295 if (typeid(*vert)!=typeid(*n.vert))
297 // Are the indices of the top-level nodes equal?
298 if (!(*vert==*n.vert))
300 // Are the sets of children equal, one by one?
301 return (children==n.children);
304 bool node::operator<(const node &n) const
306 // Are the types of the top-level node vertices different?
307 if (typeid(*vert)!=typeid(*n.vert))
308 return typeid(*vert).before(typeid(*n.vert));
309 // Are the indices of the top-level nodes different?
310 if (!(*vert==*n.vert))
311 return (vert<n.vert);
312 // Are the sets of children different, one by one?
313 return (children<n.children);
317 * These operators let us write down trees in an intuitive way, by adding
318 * arbitrarily complex children to a given vertex. The eye candy that can be
319 * produced with it makes detection of errors much simpler than with code
320 * written using calls to node's method add_child() because it allows for
321 * editor-assisted indentation.
323 const node operator+(const node &n, const child &c)
326 nn.add_child(c.first, c.second);
330 void operator+=(node &n, const child &c)
332 n.add_child(c.first, c.second);
340 static const node tree1(unsigned cuts=0)
349 * Vacuum Gamma Vacuum
353 static const node tree2(unsigned cuts=0)
357 + child(Sigma(), bool(cuts & 1))
358 + child(Sigma(), bool(cuts & 2))
359 + child(Sigma_flipped(), bool(cuts & 4)),
361 + child(Gamma(), bool(cuts & 16))
362 + child(Gamma(), bool(cuts & 32)));
375 static const node tree3(unsigned cuts=0)
381 + child(Sigma(), bool(cuts & 1)),
384 + child(Sigma(), bool(cuts & 4))
385 + child(Sigma(), bool(cuts & 8)),
395 * Sigma Sigma Sigma0 Sigma
397 static const node tree4(unsigned cuts=0)
401 + child(Sigma(), bool(cuts & 1))
402 + child(Sigma(), bool(cuts & 2)),
405 + child(Sigma_flipped(), bool(cuts & 8))
406 + child(Sigma(), bool(cuts & 16)),
412 * Sigma Vacuum Vacuum
416 static const node tree5(unsigned cuts=0)
419 + child(Sigma(), bool(cuts & 1))
421 + child(Sigma(), bool(cuts & 2))
422 + child(Sigma_flipped(), bool(cuts & 4)),
425 + child(Sigma(), bool(cuts & 16)),
437 static const node tree6(unsigned cuts=0)
441 + child(Sigma(), bool(cuts & 1)),
443 + child(Sigma_flipped()
445 + child(Vacuum(), bool(cuts & 4)),
450 static unsigned test_tree(const node tree_generator(unsigned))
452 const int edges = tree_generator(0).total_edges();
453 const int vertices = edges+1;
455 // fill a vector of all possible 2^edges combinations of cuts...
456 vector<node> counter;
457 for (unsigned i=0; i<(1U<<edges); ++i)
458 counter.push_back(tree_generator(i));
460 // ...the sum, when evaluated and reexpanded, is the antipode...
462 for (vector<node>::iterator i=counter.begin(); i!=counter.end(); ++i)
463 result = (result+i->evaluate(x,vertices-1)).series(x==0,vertices-1).expand();
465 // ...and has the nice property that in each term all the Eulers cancel:
466 if (result.has(Euler)) {
467 clog << "The antipode was miscalculated\nAntipode==" << result
468 << "\nshould not have any occurrence of Euler" << endl;
470 } else if (result.ldegree(x)!=-vertices || result.degree(x)!=0) {
471 clog << "The antipode was miscalculated\nAntipode==" << result
472 << "\nshould run from " << x << "^(" << -vertices << ") to "
473 << x << "^(0)" << "but it runs from " << x << "^("
474 << result.ldegree(x) << ")" << "to " << x << "^("
475 << result.degree(x) << ")" << endl;
481 unsigned time_antipode()
484 timer jaeger_le_coultre;
486 cout << "timing computation of antipodes in Yukawa theory" << flush;
489 jaeger_le_coultre.start();
490 result += test_tree(tree1); cout << '.' << flush;
491 result += test_tree(tree2); cout << '.' << flush;
492 result += test_tree(tree3); cout << '.' << flush;
493 result += test_tree(tree4); cout << '.' << flush;
494 result += test_tree(tree5); cout << '.' << flush;
495 result += test_tree(tree6); cout << '.' << flush;
497 cout << jaeger_le_coultre.read() << "s (total)" << endl;
499 cout << " disabled" << endl;
504 extern void randomify_symbol_serials();
506 int main(int argc, char** argv)
508 randomify_symbol_serials();
509 cout << setprecision(2) << showpoint;
510 return time_antipode();