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 Isabella Bierenbaum and Dirk Kreimer.
  13  *  For details, please see the diploma theses of Isabella Bierenbaum.
  14  */
  15
  16 /*
  17  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
  18  *
  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.
  23  *
  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.
  28  *
  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
  32  */
  33
  34 #include "times.h"
  35 #include <utility>
  36 #include <list>
  37
  38 // whether to run this beast or not:
  39 static const bool do_test = true;
  40
  41 typedef pair<unsigned, unsigned> ijpair;
  42 typedef pair<class node, bool> child;
  43
  44 const constant TrOne("Tr[One]", numeric(4));
  45
  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)
  48 {
  49         unsigned order = grad;
  50         ex exser;
  51         // maybe we have to generate more terms on the series (obnoxious):
  52         do {
  53                 exser = exarg.series(x==0, order);
  54                 ++order;
  55         } while (exser.degree(x) < 0);
  56         ex exser_trunc;
  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.
  60         return exser_trunc;
  61 }
  62
  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)
  65 {
  66         if ((i==0 && a<=0) || (j==0 && b<=0))
  67                 return 0;
  68         else
  69                 return (tgamma(2-a-(i+1)*x)*
  70                         tgamma(2-b-(1+j)*x)*
  71                         tgamma(a+b-2+(1+i+j)*x)/
  72                         tgamma(a+i*x)/
  73                         tgamma(b+j*x)/tgamma(4-a-b-(2+i+j)*x));
  74 }
  75
  76 /* Abstract base class (ABC) for all types of vertices. */
  77 class vertex {
  78 public:
  79         vertex(ijpair ij = ijpair(0,0)) : indices(ij) { }
  80         void increment_indices(const ijpair &ind) { indices.first += ind.first; indices.second += ind.second; }
  81         virtual ~vertex() { }
  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;
  85 protected:
  86         ijpair indices;
  87 };
  88
  89
  90 /*
  91  * Class of vertices of type Sigma.
  92  */
  93 class Sigma : public vertex {
  94 public:
  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;
  99 private:
 100         bool flag;
 101 };
 102
 103 ijpair Sigma::get_increment(void) const
 104 {
 105         if (flag == true)
 106             return  ijpair(indices.first+1, indices.second);
 107         else
 108             return  ijpair(indices.first, indices.second+1);
 109 }
 110
 111 ex Sigma::evaluate(const symbol &x) const
 112 {
 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;
 116 }
 117
 118
 119 /*
 120  *Class of vertices of type Gamma.
 121  */
 122 class Gamma : public vertex {
 123 public:
 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;
 128 private:
 129 };
 130
 131 ex Gamma::evaluate(const symbol &x) const
 132 {
 133         int i = indices.first;
 134         int j = indices.second;
 135         return F_ab(1,i,1,j,x);
 136 }
 137
 138
 139 /*
 140  * Class of vertices of type Vacuum.
 141  */
 142 class Vacuum : public vertex {
 143 public:
 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;
 148 private:
 149 };
 150
 151 ex Vacuum::evaluate(const symbol &x) const
 152 {
 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;
 156 }
 157
 158
 159 /*
 160  * Class of nodes (or trees or subtrees), including list of children.
 161  */
 162 class node {
 163 public:
 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;
 171 private:
 172         vertex *vert;
 173         list<child> children;
 174 };
 175
 176 const node & node::operator=(const node &n)
 177 {
 178         delete vert;
 179         vert = (n.vert)->copy();
 180         children = n.children;
 181 }
 182
 183 void node::add_child(const node &childnode, bool cut)
 184 {
 185         children.push_back(child(childnode, cut));
 186         if(!cut)
 187                 vert->increment_indices(childnode.vert->get_increment());
 188 }
 189
 190 ex node::evaluate(const symbol &x, unsigned grad) const
 191 {
 192         ex product = 1;
 193         for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
 194                 if (!i->second)
 195                         product *= i->first.evaluate(x,grad);
 196                 else
 197                         product *= -div_part(i->first.evaluate(x,grad),x,grad);
 198         }
 199         return (product * vert->evaluate(x));
 200 }
 201
 202 unsigned node::total_edges(void) const
 203 {
 204         unsigned accu = 0;
 205         for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
 206                 accu += i->first.total_edges();
 207                 ++accu;
 208         }
 209         return accu;
 210 }
 211
 212
 213 /*
 214  * These operators let us write down trees in an intuitive way, by adding
 215  * arbitrarily complex children to a given vertex.  The eye candy that can be
 216  * produced with it makes detection of errors much simpler than with code
 217  * written using calls to node's method add_child() because it allows for
 218  * editor-assisted indentation.
 219  */
 220 const node operator+(const node &n, const child &c)
 221 {
 222         node nn(n);
 223         nn.add_child(c.first, c.second);
 224         return nn;
 225 }
 226 void operator+=(node &n, const child &c)
 227 {
 228         n.add_child(c.first, c.second);
 229 }
 230
 231 /*
 232  * Build this sample rooted tree characterized by a certain combination of
 233  * cut or uncut edges as specified by the unsigned parameter:
 234  *              Gamma
 235  *              /   \
 236  *         Sigma     Vacuum
 237  *        /   \       /   \
 238  *    Sigma Sigma  Sigma0 Sigma
 239  */
 240 static node mytree(unsigned cuts=0)
 241 {
 242         return (Gamma()
 243                 + child(Sigma()
 244                         + child(Sigma(), bool(cuts & 1))
 245                         + child(Sigma(), bool(cuts & 2)),
 246                         bool(cuts & 4))
 247                 + child(Vacuum()
 248                         + child(Sigma(ijpair(0,0),false), bool(cuts & 8))
 249                         + child(Sigma(), bool(cuts & 16)),
 250                         bool(cuts & 32)));
 251 }
 252
 253
 254 static unsigned test(void)
 255 {
 256         const symbol x("x");
 257
 258         const unsigned edges = mytree().total_edges();
 259         const unsigned vertices = edges+1;
 260
 261         // fill a vector of all possible 2^edges combinations of cuts...
 262         vector<node> counter;
 263         for (unsigned i=0; i<(1U<<edges); ++i)
 264                 counter.push_back(mytree(i));
 265
 266         // ...the sum, when evaluated, is the antipode...
 267         ex accu = 0;
 268         for (vector<node>::iterator i=counter.begin(); i!=counter.end(); ++i)
 269                 accu += i->evaluate(x,vertices);
 270
 271         // ...which is only interesting term-wise in the series expansion...
 272         ex result = accu.series(x==0,vertices).expand().normal();
 273
 274         // ...and has the nice property that in each term all the Eulers cancel:
 275         if (result.has(Euler)) {
 276                 clog << "The antipode was miscalculated\nAntipode==" << result
 277                      << "\nshould not have any occurrence of Euler" << endl;
 278                 return 1;
 279         }
 280         return 0;
 281 }
 282
 283 unsigned time_antipode(void)
 284 {
 285         unsigned result = 0;
 286         unsigned count = 0;
 287         timer jaeger_le_coultre;
 288         double time = .0;
 289
 290         cout << "timing computation of an antipode in Yukawa theory" << flush;
 291         clog << "-------computation of an antipode in Yukawa theory" << endl;
 292
 293         if (do_test) {
 294                 jaeger_le_coultre.start();
 295                 // correct for very small times:
 296                 do {
 297                         result = test();
 298                         ++count;
 299                 } while ((time=jaeger_le_coultre.read())<0.1 && !result);
 300                 cout << '.' << flush;
 301
 302                 if (!result) {
 303                         cout << " passed ";
 304                         clog << "(no output)" << endl;
 305                 } else {
 306                         cout << " failed ";
 307                 }
 308                 cout << int(1000*(time/count))*0.001 << 's' << endl;
 309         } else {
 310                 cout << " disabled" << endl;
 311                 clog << "(no output)" << endl;
 312         }
 313
 314         return result;
 315 }