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  */
  17
  18 /*
  19  *  GiNaC Copyright (C) 1999-2001 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  34  */
  35
  36 #include "times.h"
  37 #include <utility>
  38 #include <vector>
  39 #include <set>
  40 #include <map>
  41 #include <typeinfo>
  42 #include <stdexcept>
  43
  44 // whether to run this beast or not:
  45 static const bool do_test = true;
  46
  47 // regularization parameter:
  48 static const symbol x("x");
  49
  50 typedef pair<unsigned, unsigned> ijpair;
  51 typedef pair<class node, bool> child;
  52
  53 const constant TrOne("Tr[One]", numeric(4));
  54
  55 /* Extract only the divergent part of a series and discard the rest. */
  56 static ex div_part(const ex &exarg, const symbol &x, unsigned grad)
  57 {
  58         const ex exser = exarg.series(x==0, grad);
  59         if (exser.degree(x)<0)
  60                 throw runtime_error("divergent part truncation disaster");
  61         ex exser_trunc;
  62         for (int i=exser.ldegree(x); i<0; ++i)
  63                 exser_trunc += exser.coeff(x,i)*pow(x,i);
  64         // NB: exser_trunc is by construction collected in x.
  65         return exser_trunc;
  66 }
  67
  68 /* F_ab(a, i, b, j, "x") is a common pattern in all vertex evaluators. */
  69 static ex F_ab(int a, int i, int b, int j, const symbol &x)
  70 {
  71         if ((i==0 && a<=0) || (j==0 && b<=0))
  72                 return 0;
  73         else
  74                 return (tgamma(2-a-(i+1)*x)*
  75                         tgamma(2-b-(1+j)*x)*
  76                         tgamma(a+b-2+(1+i+j)*x)/
  77                         tgamma(a+i*x)/
  78                         tgamma(b+j*x)/tgamma(4-a-b-(2+i+j)*x));
  79 }
  80
  81 /* Abstract base class (ABC) for all types of vertices. */
  82 class vertex {
  83 public:
  84         vertex(ijpair ij = ijpair(0,0)) : indices(ij) { }
  85         void increment_indices(const ijpair &ind) { indices.first += ind.first; indices.second += ind.second; }
  86         virtual ~vertex() { }
  87         virtual vertex* copy(void) const = 0;
  88         virtual ijpair get_increment(void) const { return indices; }
  89         virtual const ex evaluate(const symbol &x, const unsigned grad) const = 0;
  90         bool operator==(const vertex &v) const { return (indices==v.indices); }
  91         bool operator<(const vertex &v) const { return (indices<v.indices); }
  92 protected:
  93         ijpair indices;
  94 };
  95
  96
  97 /*
  98  * Class of vertices of type Sigma.
  99  */
 100 class Sigma : public vertex {
 101 public:
 102         Sigma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
 103         vertex* copy(void) const { return new Sigma(*this); }
 104         ijpair get_increment(void) const { return ijpair(indices.first+indices.second+1, 0); }
 105         const ex evaluate(const symbol &x, const unsigned grad) const;
 106 private:
 107 };
 108
 109 const ex Sigma::evaluate(const symbol &x, const unsigned grad) const
 110 {
 111         // c.f. comments in node::evaluate()
 112         static map<Sigma,ex> catalog;
 113         static unsigned prev_grad = 0;
 114         if (grad>prev_grad) {
 115                 catalog.clear();
 116                 prev_grad = grad;
 117         }
 118         map<Sigma,ex>::iterator pos = catalog.find(*this);
 119         if (pos==catalog.end()) {
 120                 int i = indices.first;
 121                 int j = indices.second;
 122                 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();
 123                 pos = catalog.insert(map<Sigma,ex>::value_type(*this,result)).first;
 124                 if (grad<prev_grad)
 125                         prev_grad = grad;
 126         }
 127         return pos->second;
 128 }
 129
 130
 131 /** Class of vertices of type Sigma_flipped, sitting in the upper fermionline of Vacuum; no consequences for Gamma. */
 132 class Sigma_flipped : public Sigma {
 133 public:
 134         Sigma_flipped(ijpair ij = ijpair(0,0)) : Sigma(ij) { }
 135         vertex* copy(void) const { return new Sigma_flipped(*this); }
 136         ijpair get_increment(void) const { return ijpair(0, indices.first+indices.second+1); }
 137         const ex evaluate(const symbol &x, const unsigned grad) const { return Sigma::evaluate(x, grad); }
 138 private:
 139 };
 140
 141
 142 /*
 143  *Class of vertices of type Gamma.
 144  */
 145 class Gamma : public vertex {
 146 public:
 147         Gamma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
 148         vertex* copy(void) const { return new Gamma(*this); }
 149         ijpair get_increment(void) const { return ijpair(indices.first+indices.second+1, 0); }
 150         const ex evaluate(const symbol &x, const unsigned grad) const;
 151 private:
 152 };
 153
 154 const ex Gamma::evaluate(const symbol &x, const unsigned grad) const
 155 {
 156         // c.f. comments in node::evaluate()
 157         static map<Gamma,ex> catalog;
 158         static unsigned prev_grad = 0;
 159         if (prev_grad<grad) {
 160                 catalog.clear();
 161                 prev_grad = grad;
 162         }
 163         map<Gamma,ex>::iterator pos = catalog.find(*this);
 164         if (pos==catalog.end()) {
 165                 int i = indices.first;
 166                 int j = indices.second;
 167                 const ex result = (F_ab(1,i,1,j,x)).series(x==0,grad).expand();
 168                 pos = catalog.insert(map<Gamma,ex>::value_type(*this,result)).first;
 169                 if (grad<prev_grad)
 170                         prev_grad = grad;
 171         }
 172         return pos->second;
 173 }
 174
 175
 176 /*
 177  * Class of vertices of type Vacuum.
 178  */
 179 class Vacuum : public vertex {
 180 public:
 181         Vacuum(ijpair ij = ijpair(0,0)) : vertex(ij) { }
 182         vertex* copy(void) const { return new Vacuum(*this); }
 183         ijpair get_increment() const { return ijpair(0, indices.first+indices.second+1); }
 184         const ex evaluate(const symbol &x, const unsigned grad) const;
 185 private:
 186 };
 187
 188 const ex Vacuum::evaluate(const symbol &x, const unsigned grad) const
 189 {
 190         // c.f. comments in node::evaluate()
 191         static map<Vacuum,ex> catalog;
 192         static unsigned prev_grad = 0;
 193         if (prev_grad<grad) {
 194                 catalog.clear();
 195                 prev_grad = grad;
 196         }
 197         map<Vacuum,ex>::iterator pos = catalog.find(*this);
 198         if (pos==catalog.end()) {
 199                 int i = indices.first;
 200                 int j = indices.second;
 201                 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();
 202                 pos = catalog.insert(map<Vacuum,ex>::value_type(*this,result)).first;
 203                 if (grad<prev_grad)
 204                         prev_grad = grad;
 205         }
 206         return pos->second;
 207 }
 208
 209
 210 /*
 211  * Class of nodes (or trees or subtrees), including list of children.
 212  */
 213 class node {
 214 public:
 215         node(const vertex &v) { vert = v.copy(); }
 216         node(const node &n) { vert = (n.vert)->copy(); children = n.children; }
 217         const node & operator=(const node &);
 218         ~node() { delete vert; }
 219         void add_child(const node &, bool = false);
 220         const ex evaluate(const symbol &x, unsigned grad) const;
 221         unsigned total_edges(void) const;
 222         bool operator==(const node &) const;
 223         bool operator<(const node &) const;
 224 private:
 225         vertex *vert;
 226         multiset<child> children;
 227 };
 228
 229 const node & node::operator=(const node &n)
 230 {
 231         if (this!=&n) {
 232                 delete vert;
 233                 vert = (n.vert)->copy();
 234                 children = n.children;
 235         }
 236         return *this;
 237 }
 238
 239 void node::add_child(const node &childnode, bool cut)
 240 {
 241         children.insert(child(childnode, cut));
 242         if(!cut)
 243                 vert->increment_indices(childnode.vert->get_increment());
 244 }
 245
 246 const ex node::evaluate(const symbol &x, unsigned grad) const
 247 {
 248         static map<node,ex> catalog;    // lookup table for already evaluated subnodes
 249         static unsigned prev_grad = 0;  // grad argument that the catalog has been build for
 250         if (grad>prev_grad) {
 251                 // We haven't computed far enough last time.  Our catalog cannot cope with
 252                 // the demands for this value of grad so let's clear it.
 253                 catalog.clear();
 254                 prev_grad = grad;
 255         }
 256         ex product = 1;   // accumulator for all the children
 257         for (multiset<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
 258                 map<node,ex>::iterator pos = catalog.find(i->first);
 259                 if (pos==catalog.end()) {
 260                         pos = catalog.insert(map<node,ex>::value_type(i->first,i->first.evaluate(x,grad).series(x==0,grad).expand())).first;
 261                         if (grad<prev_grad) {
 262                                 // We have just spoiled the catalog by inserting a series computed
 263                                 // to lower grad than the others in it.  So let's make sure next time
 264                                 // we don't use one of the newly inserted ones by bumping prev_grad
 265                                 // down to the current value of grad.
 266                                 prev_grad = grad;
 267                         }
 268                 }
 269                 if (!i->second)
 270                         product *= pos->second;
 271                 else
 272                         product *= -div_part(pos->second,x,grad);
 273         }
 274         return (product * vert->evaluate(x,grad));
 275 }
 276
 277 unsigned node::total_edges(void) const
 278 {
 279         unsigned accu = 0;
 280         for (multiset<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
 281                 accu += i->first.total_edges();
 282                 ++accu;
 283         }
 284         return accu;
 285 }
 286
 287 bool node::operator==(const node &n) const
 288 {
 289         // Are the types of the top-level node vertices equal?
 290         if (typeid(*vert)!=typeid(*n.vert))
 291                 return false;
 292         // Are the indices of the top-level nodes equal?
 293         if (!(*vert==*n.vert))
 294                 return false;
 295         // Are the sets of children equal, one by one?
 296         return (children==n.children);
 297 }
 298
 299 bool node::operator<(const node &n) const
 300 {
 301         // Are the types of the top-level node vertices different?
 302         if (typeid(*vert)!=typeid(*n.vert))
 303                 return typeid(*vert).before(typeid(*n.vert));
 304         // Are the indices of the top-level nodes different?
 305         if (!(*vert==*n.vert))
 306                 return (vert<n.vert);
 307         // Are the sets of children different, one by one?
 308         return (children<n.children);
 309 }
 310
 311 /*
 312  * These operators let us write down trees in an intuitive way, by adding
 313  * arbitrarily complex children to a given vertex.  The eye candy that can be
 314  * produced with it makes detection of errors much simpler than with code
 315  * written using calls to node's method add_child() because it allows for
 316  * editor-assisted indentation.
 317  */
 318 const node operator+(const node &n, const child &c)
 319 {
 320         node nn(n);
 321         nn.add_child(c.first, c.second);
 322         return nn;
 323 }
 324
 325 void operator+=(node &n, const child &c)
 326 {
 327         n.add_child(c.first, c.second);
 328 }
 329
 330
 331 /*                Gamma
 332  *                  |
 333  *                Gamma
 334  */
 335 static const node tree1(unsigned cuts=0)
 336 {
 337         return (Gamma()
 338                 + child(Gamma(),
 339                         bool(cuts & 1)));
 340 }
 341
 342 /*                Gamma
 343  *               /  |  \
 344  *        Vacuum  Gamma  Vacuum
 345  *       /   |  \
 346  *  Sigma Sigma Sigma0
 347  */
 348 static const node tree2(unsigned cuts=0)
 349 {
 350         return (Gamma()
 351                 + child(Vacuum()
 352                         + child(Sigma(), bool(cuts & 1))
 353                         + child(Sigma(), bool(cuts & 2))
 354                         + child(Sigma_flipped(), bool(cuts & 4)),
 355                         bool(cuts & 8))
 356                 + child(Gamma(), bool(cuts & 16))
 357                 + child(Gamma(), bool(cuts & 32)));
 358 }
 359
 360 /*                Gamma
 361  *                  |
 362  *                Gamma
 363  *                  |
 364  *                Gamma
 365  *                /   \
 366  *           Vacuum    Gamma
 367  *           /    \       \
 368  *        Sigma  Sigma   Sigma
 369  */
 370 static const node tree3(unsigned cuts=0)
 371 {
 372         return (Gamma()
 373                 + child(Gamma()
 374                         + child(Gamma()
 375                                 + child(Gamma()
 376                                         + child(Sigma(), bool(cuts & 1)),
 377                                         bool(cuts & 2))
 378                                 + child(Vacuum()
 379                                         + child(Sigma(), bool(cuts & 4))
 380                                         + child(Sigma(), bool(cuts & 8)),
 381                                 bool(cuts & 16)),
 382                         bool(cuts & 32)),
 383                 bool(cuts & 64)));
 384 }
 385
 386 /*                Gamma
 387  *                /   \
 388  *           Sigma     Vacuum
 389  *          /   \       /   \
 390  *      Sigma Sigma  Sigma0 Sigma
 391  */
 392 static const node tree4(unsigned cuts=0)
 393 {
 394         return (Gamma()
 395                 + child(Sigma()
 396                         + child(Sigma(), bool(cuts & 1))
 397                         + child(Sigma(), bool(cuts & 2)),
 398                         bool(cuts & 4))
 399                 + child(Vacuum()
 400                         + child(Sigma_flipped(), bool(cuts & 8))
 401                         + child(Sigma(), bool(cuts & 16)),
 402                         bool(cuts & 32)));
 403 }
 404
 405 /*                Sigma
 406  *               /  |  \
 407  *         Sigma Vacuum  Vacuum
 408  *               /    \       \
 409  *            Sigma  Sigma0   Sigma
 410  */
 411 static const node tree5(unsigned cuts=0)
 412 {
 413         return (Sigma()
 414                 + child(Sigma(), bool(cuts & 1))
 415                 + child(Vacuum()
 416                         + child(Sigma(), bool(cuts & 2))
 417                         + child(Sigma_flipped(), bool(cuts & 4)),
 418                         bool(cuts & 8))
 419                 + child(Vacuum()
 420                         + child(Sigma(), bool(cuts & 16)),
 421                         bool(cuts & 32)));
 422 }
 423
 424 /*               Vacuum
 425  *               /    \
 426  *           Sigma    Sigma0
 427  *             |        |
 428  *           Sigma    Sigma
 429  *                      |
 430  *                    Vacuum
 431  */
 432 static const node tree6(unsigned cuts=0)
 433 {
 434         return (Vacuum()
 435                 + child(Sigma()
 436                         + child(Sigma(), bool(cuts & 1)),
 437                         bool(cuts & 2))
 438                 + child(Sigma_flipped()
 439                         + child(Sigma()
 440                                 + child(Vacuum(), bool(cuts & 4)),
 441                                 bool(cuts & 8)),
 442                         bool(cuts & 16)));
 443 }
 444
 445 static unsigned test_tree(const node (*tree_generator)(unsigned=0))
 446 {
 447         const int edges = tree_generator().total_edges();
 448         const int vertices = edges+1;
 449
 450         // fill a vector of all possible 2^edges combinations of cuts...
 451         vector<node> counter;
 452         for (unsigned i=0; i<(1U<<edges); ++i)
 453                 counter.push_back(tree_generator(i));
 454
 455         // ...the sum, when evaluated and reexpanded, is the antipode...
 456         ex result = 0;
 457         for (vector<node>::iterator i=counter.begin(); i!=counter.end(); ++i)
 458                 result = (result+i->evaluate(x,vertices)).series(x==0,vertices).expand();
 459
 460         // ...and has the nice property that in each term all the Eulers cancel:
 461         if (result.has(Euler)) {
 462                 clog << "The antipode was miscalculated\nAntipode==" << result
 463                      << "\nshould not have any occurrence of Euler" << endl;
 464                 return 1;
 465         } else if (result.ldegree(x)!=-vertices || result.degree(x)!=0) {
 466                 clog << "The antipode was miscalculated\nAntipode==" << result
 467                      << "\nshould run from x^(" << -vertices << ") to x^(0) but it runs from x^("
 468                      << result.ldegree(x) << ") to x^(" << result.degree(x) << ")" << endl;
 469                 return 1;
 470         }
 471         return 0;
 472 }
 473
 474 unsigned time_antipode(void)
 475 {
 476         unsigned result = 0;
 477         timer jaeger_le_coultre;
 478
 479         cout << "timing computation of antipodes in Yukawa theory" << flush;
 480         clog << "-------computation of antipodes in Yukawa theory" << endl;
 481
 482         if (do_test) {
 483                 jaeger_le_coultre.start();
 484                 result += test_tree(tree1);  cout << '.' << flush;
 485                 result += test_tree(tree2);  cout << '.' << flush;
 486                 result += test_tree(tree3);  cout << '.' << flush;
 487                 result += test_tree(tree4);  cout << '.' << flush;
 488                 result += test_tree(tree5);  cout << '.' << flush;
 489                 result += test_tree(tree6);  cout << '.' << flush;
 490
 491                 if (!result) {
 492                         cout << " passed ";
 493                         clog << "(no output)" << endl;
 494                 } else {
 495                         cout << " failed ";
 496                 }
 497                 cout << int(1000*jaeger_le_coultre.read())*0.001 << "s (total)" << endl;
 498         } else {
 499                 cout << " disabled" << endl;
 500                 clog << "(no output)" << endl;
 501         }
 502
 503         return result;
 504 }