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 see <http://www.arXiv.org/abs/hep-th/0111192>.
  16  */
  17
  18 /*
  19  *  GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
  20  *
  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.
  25  *
  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.
  30  *
  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
  34  */
  35
  36 #include "ginac.h"
  37 #include "timer.h"
  38 using namespace GiNaC;
  39
  40 #include <map>
  41 #include <set>
  42 #include <stdexcept>
  43 #include <typeinfo>
  44 #include <utility>
  45 #include <vector>
  46 using namespace std;
  47
  48 // whether to run this beast or not:
  49 constexpr bool do_test = true;
  50
  51 // regularization parameter:
  52 static const symbol x("x");
  53
  54 typedef pair<unsigned, unsigned> ijpair;
  55 typedef pair<class node, bool> child;
  56
  57 const constant TrOne("Tr[One]", numeric(4));
  58
  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)
  61 {
  62         const ex exser = exarg.series(x==0, grad);
  63         if (exser.degree(x)<0)
  64                 throw runtime_error("divergent part truncation disaster");
  65         ex exser_trunc;
  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.
  69         return exser_trunc;
  70 }
  71
  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)
  74 {
  75         using GiNaC::tgamma;
  76         if ((i==0 && a<=0) || (j==0 && b<=0))
  77                 return 0;
  78         else
  79                 return (tgamma(2-a-(i+1)*x)*
  80                         tgamma(2-b-(1+j)*x)*
  81                         tgamma(a+b-2+(1+i+j)*x)/
  82                         tgamma(a+i*x)/
  83                         tgamma(b+j*x)/tgamma(4-a-b-(2+i+j)*x));
  84 }
  85
  86 /* Abstract base class (ABC) for all types of vertices. */
  87 class vertex {
  88 public:
  89         vertex(ijpair ij = ijpair(0,0)) : indices(ij) { }
  90         void increment_indices(const ijpair &ind) { indices.first += ind.first; indices.second += ind.second; }
  91         virtual ~vertex() { }
  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); }
  97 protected:
  98         ijpair indices;
  99 };
 100
 101
 102 /*
 103  * Class of vertices of type Sigma.
 104  */
 105 class Sigma : public vertex {
 106 public:
 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;
 111 private:
 112 };
 113
 114 const ex Sigma::evaluate(const symbol &x, const unsigned grad) const
 115 {
 116         // c.f. comments in node::evaluate()
 117         static map<Sigma,ex> catalog;
 118         static unsigned prev_grad = 0;
 119         if (grad>prev_grad) {
 120                 catalog.clear();
 121                 prev_grad = grad;
 122         }
 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;
 129                 if (grad<prev_grad)
 130                         prev_grad = grad;
 131         }
 132         return pos->second;
 133 }
 134
 135
 136 /** Class of vertices of type Sigma_flipped, sitting in the upper fermionline of Vacuum; no consequences for Gamma. */
 137 class Sigma_flipped : public Sigma {
 138 public:
 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); }
 143 private:
 144 };
 145
 146
 147 /*
 148  *Class of vertices of type Gamma.
 149  */
 150 class Gamma : public vertex {
 151 public:
 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;
 156 private:
 157 };
 158
 159 const ex Gamma::evaluate(const symbol &x, const unsigned grad) const
 160 {
 161         // c.f. comments in node::evaluate()
 162         static map<Gamma,ex> catalog;
 163         static unsigned prev_grad = 0;
 164         if (prev_grad<grad) {
 165                 catalog.clear();
 166                 prev_grad = grad;
 167         }
 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;
 174                 if (grad<prev_grad)
 175                         prev_grad = grad;
 176         }
 177         return pos->second;
 178 }
 179
 180
 181 /*
 182  * Class of vertices of type Vacuum.
 183  */
 184 class Vacuum : public vertex {
 185 public:
 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;
 190 private:
 191 };
 192
 193 const ex Vacuum::evaluate(const symbol &x, const unsigned grad) const
 194 {
 195         // c.f. comments in node::evaluate()
 196         static map<Vacuum,ex> catalog;
 197         static unsigned prev_grad = 0;
 198         if (prev_grad<grad) {
 199                 catalog.clear();
 200                 prev_grad = grad;
 201         }
 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;
 208                 if (grad<prev_grad)
 209                         prev_grad = grad;
 210         }
 211         return pos->second;
 212 }
 213
 214
 215 /*
 216  * Class of nodes (or trees or subtrees), including list of children.
 217  */
 218 class node {
 219 public:
 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;
 229 private:
 230         vertex *vert;
 231         multiset<child> children;
 232 };
 233
 234 const node & node::operator=(const node &n)
 235 {
 236         if (this!=&n) {
 237                 delete vert;
 238                 vert = (n.vert)->copy();
 239                 children = n.children;
 240         }
 241         return *this;
 242 }
 243
 244 void node::add_child(const node &childnode, bool cut)
 245 {
 246         children.insert(child(childnode, cut));
 247         if(!cut)
 248                 vert->increment_indices(childnode.vert->get_increment());
 249 }
 250
 251 const ex node::evaluate(const symbol &x, unsigned grad) const
 252 {
 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.
 258                 catalog.clear();
 259                 prev_grad = grad;
 260         }
 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.
 271                                 prev_grad = grad;
 272                         }
 273                 }
 274                 if (!i->second)
 275                         product *= pos->second;
 276                 else
 277                         product *= -div_part(pos->second,x,grad);
 278         }
 279         return (product * vert->evaluate(x,grad));
 280 }
 281
 282 unsigned node::total_edges() const
 283 {
 284         unsigned accu = 0;
 285         for (multiset<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
 286                 accu += i->first.total_edges();
 287                 ++accu;
 288         }
 289         return accu;
 290 }
 291
 292 bool node::operator==(const node &n) const
 293 {
 294         // Are the types of the top-level node vertices equal?
 295         if (typeid(*vert)!=typeid(*n.vert))
 296                 return false;
 297         // Are the indices of the top-level nodes equal?
 298         if (!(*vert==*n.vert))
 299                 return false;
 300         // Are the sets of children equal, one by one?
 301         return (children==n.children);
 302 }
 303
 304 bool node::operator<(const node &n) const
 305 {
 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);
 314 }
 315
 316 /*
 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.
 322  */
 323 const node operator+(const node &n, const child &c)
 324 {
 325         node nn(n);
 326         nn.add_child(c.first, c.second);
 327         return nn;
 328 }
 329
 330 void operator+=(node &n, const child &c)
 331 {
 332         n.add_child(c.first, c.second);
 333 }
 334
 335
 336 /*                Gamma
 337  *                  |
 338  *                Gamma
 339  */
 340 static const node tree1(unsigned cuts=0)
 341 {
 342         return (Gamma()
 343                 + child(Gamma(),
 344                         bool(cuts & 1)));
 345 }
 346
 347 /*                Gamma
 348  *               /  |  \
 349  *        Vacuum  Gamma  Vacuum
 350  *       /   |  \
 351  *  Sigma Sigma Sigma0
 352  */
 353 static const node tree2(unsigned cuts=0)
 354 {
 355         return (Gamma()
 356                 + child(Vacuum()
 357                         + child(Sigma(), bool(cuts & 1))
 358                         + child(Sigma(), bool(cuts & 2))
 359                         + child(Sigma_flipped(), bool(cuts & 4)),
 360                         bool(cuts & 8))
 361                 + child(Gamma(), bool(cuts & 16))
 362                 + child(Gamma(), bool(cuts & 32)));
 363 }
 364
 365 /*                Gamma
 366  *                  |
 367  *                Gamma
 368  *                  |
 369  *                Gamma
 370  *                /   \
 371  *           Vacuum    Gamma
 372  *           /    \       \
 373  *        Sigma  Sigma   Sigma
 374  */
 375 static const node tree3(unsigned cuts=0)
 376 {
 377         return (Gamma()
 378                 + child(Gamma()
 379                         + child(Gamma()
 380                                 + child(Gamma()
 381                                         + child(Sigma(), bool(cuts & 1)),
 382                                         bool(cuts & 2))
 383                                 + child(Vacuum()
 384                                         + child(Sigma(), bool(cuts & 4))
 385                                         + child(Sigma(), bool(cuts & 8)),
 386                                 bool(cuts & 16)),
 387                         bool(cuts & 32)),
 388                 bool(cuts & 64)));
 389 }
 390
 391 /*                Gamma
 392  *                /   \
 393  *           Sigma     Vacuum
 394  *          /   \       /   \
 395  *      Sigma Sigma  Sigma0 Sigma
 396  */
 397 static const node tree4(unsigned cuts=0)
 398 {
 399         return (Gamma()
 400                 + child(Sigma()
 401                         + child(Sigma(), bool(cuts & 1))
 402                         + child(Sigma(), bool(cuts & 2)),
 403                         bool(cuts & 4))
 404                 + child(Vacuum()
 405                         + child(Sigma_flipped(), bool(cuts & 8))
 406                         + child(Sigma(), bool(cuts & 16)),
 407                         bool(cuts & 32)));
 408 }
 409
 410 /*                Sigma
 411  *               /  |  \
 412  *         Sigma Vacuum  Vacuum
 413  *               /    \       \
 414  *            Sigma  Sigma0   Sigma
 415  */
 416 static const node tree5(unsigned cuts=0)
 417 {
 418         return (Sigma()
 419                 + child(Sigma(), bool(cuts & 1))
 420                 + child(Vacuum()
 421                         + child(Sigma(), bool(cuts & 2))
 422                         + child(Sigma_flipped(), bool(cuts & 4)),
 423                         bool(cuts & 8))
 424                 + child(Vacuum()
 425                         + child(Sigma(), bool(cuts & 16)),
 426                         bool(cuts & 32)));
 427 }
 428
 429 /*               Vacuum
 430  *               /    \
 431  *           Sigma    Sigma0
 432  *             |        |
 433  *           Sigma    Sigma
 434  *                      |
 435  *                    Vacuum
 436  */
 437 static const node tree6(unsigned cuts=0)
 438 {
 439         return (Vacuum()
 440                 + child(Sigma()
 441                         + child(Sigma(), bool(cuts & 1)),
 442                         bool(cuts & 2))
 443                 + child(Sigma_flipped()
 444                         + child(Sigma()
 445                                 + child(Vacuum(), bool(cuts & 4)),
 446                                 bool(cuts & 8)),
 447                         bool(cuts & 16)));
 448 }
 449
 450 static unsigned test_tree(const node tree_generator(unsigned))
 451 {
 452         const int edges = tree_generator(0).total_edges();
 453         const int vertices = edges+1;
 454
 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));
 459
 460         // ...the sum, when evaluated and reexpanded, is the antipode...
 461         ex result = 0;
 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();
 464
 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;
 469                 return 1;
 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;
 476                 return 1;
 477         }
 478         return 0;
 479 }
 480
 481 unsigned time_antipode()
 482 {
 483         unsigned result = 0;
 484         timer jaeger_le_coultre;
 485
 486         cout << "timing computation of antipodes in Yukawa theory" << flush;
 487
 488         if (do_test) {
 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;
 496
 497                 cout << jaeger_le_coultre.read() << "s (total)" << endl;
 498         } else {
 499                 cout << " disabled" << endl;
 500         }
 501         return result;
 502 }
 503
 504 extern void randomify_symbol_serials();
 505
 506 int main(int argc, char** argv)
 507 {
 508         randomify_symbol_serials();
 509         cout << setprecision(2) << showpoint;
 510         return time_antipode();
 511 }