* New funny timing added: compute an antipode in Yukawa theory.
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Mon, 5 Mar 2001 20:02:00 +0000 (20:02 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Mon, 5 Mar 2001 20:02:00 +0000 (20:02 +0000)
check/Makefile.am
check/time_antipode.cpp [new file with mode: 0644]
check/times.cpp
check/times.h
check/times.ref

index 7932f1c63bce3927ae8f2e06ef9a01912fbf8657..a5e781b95ac3fdba4626d9117457d0fd67ef8e92 100644 (file)
@@ -14,7 +14,7 @@ times_SOURCES = time_dennyfliegner.cpp time_gammaseries.cpp \
   time_lw_C.cpp time_lw_D.cpp time_lw_E.cpp time_lw_F.cpp time_lw_G.cpp \
   time_lw_H.cpp time_lw_IJKL.cpp time_lw_M1.cpp time_lw_M2.cpp time_lw_N.cpp \
   time_lw_O.cpp time_lw_P.cpp time_lw_Pprime.cpp time_lw_Q.cpp \
-  time_lw_Qprime.cpp timer.cpp times.cpp times.h
+  time_lw_Qprime.cpp time_antipode.cpp timer.cpp times.cpp times.h
 times_LDADD = ../ginac/libginac.la
 INCLUDES = -I$(srcdir)/../ginac -I../ginac
 CLEANFILES = exams.out checks.out times.out
diff --git a/check/time_antipode.cpp b/check/time_antipode.cpp
new file mode 100644 (file)
index 0000000..2d3a13a
--- /dev/null
@@ -0,0 +1,307 @@
+/** @file time_antipode.cpp
+ *
+ *  This is a beautiful example that calculates the counterterm for the
+ *  overall divergence of some special sorts of Feynman diagrams in a massless
+ *  Yukawa theory.  For this end it computes the antipode of the corresponding
+ *  decorated rooted tree using dimensional regularization in the parameter
+ *  x==-(D-4)/2, which leads to a Laurent series in x.  The renormalization
+ *  scheme used is the minimal subtraction scheme (MS).  From an efficiency
+ *  point of view it boils down to power series expansion.  It also has quite
+ *  an intriguing check for consistency, which is why we include it here.
+ *
+ *  This program is based on work by Isabella Bierenbaum and Dirk Kreimer.
+ *  For details, please see the diploma theses of Isabella Bierenbaum.
+ */
+
+/*
+ *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "times.h"
+#include <utility>
+#include <list>
+
+// whether to run this beast or not:
+static const bool do_test = true;
+
+typedef pair<unsigned, unsigned> ijpair;
+typedef pair<class node, bool> child;
+
+const constant TrOne("Tr[One]", numeric(4));
+
+/* Extract only the divergent part of a series and discard the rest. */
+static ex div_part(const ex &exarg, const symbol &x, unsigned grad)
+{
+       unsigned order = grad;
+       ex exser;
+       // maybe we have to generate more terms on the series (obnoxious):
+       do {
+               exser = exarg.series(x==0, order);
+               ++order;
+       } while (exser.degree(x) < 0);
+       ex exser_trunc;
+       for (int i=exser.ldegree(x); i<0; ++i)
+               exser_trunc += exser.coeff(x,i)*pow(x,i);
+       // NB: exser_trunc is by construction collected in x.
+       return exser_trunc;
+}
+
+/* F_ab(a, i, b, j, "x") is a common pattern in all vertex evaluators. */
+static ex F_ab(int a, int i, int b, int j, const symbol &x)
+{
+       if ((i==0 && a<=0) || (j==0 && b<=0))
+               return 0;
+       else
+               return (tgamma(2-a-(i+1)*x)*
+                       tgamma(2-b-(1+j)*x)*
+                       tgamma(a+b-2+(1+i+j)*x)/
+                       tgamma(a+i*x)/
+                       tgamma(b+j*x)/tgamma(4-a-b-(2+i+j)*x));
+}
+
+/* Abstract base class (ABC) for all types of vertices. */
+class vertex {
+public:
+       vertex(ijpair ij = ijpair(0,0)) : indices(ij) { }
+       void increment_indices(const ijpair &ind) { indices.first += ind.first; indices.second += ind.second; }
+       virtual ~vertex() { }
+       virtual vertex* copy(void) const = 0;
+       virtual ijpair get_increment(void) const { return indices; }
+       virtual ex evaluate(const symbol &x) const = 0; 
+protected:
+       ijpair indices;
+};
+
+
+/*
+ * Class of vertices of type Sigma.
+ */
+class Sigma : public vertex {
+public:
+       Sigma(ijpair ij = ijpair(0,0), bool f = true) : vertex(ij), flag(f) { }
+       vertex* copy(void) const { return new Sigma(*this); }
+       ijpair get_increment(void) const;
+       ex evaluate(const symbol &x) const;
+private:
+       bool flag;
+};
+
+ijpair Sigma::get_increment(void) const
+{
+       if (flag == true)
+           return  ijpair(indices.first+1, indices.second);
+       else
+           return  ijpair(indices.first, indices.second+1);
+}
+
+ex Sigma::evaluate(const symbol &x) const
+{
+       int i = indices.first;
+       int j = indices.second;
+       return (F_ab(0,i,1,j,x)+F_ab(1,i,1,j,x)-F_ab(1,i,0,j,x))/2;
+}
+
+
+/*
+ *Class of vertices of type Gamma.
+ */
+class Gamma : public vertex {
+public:
+       Gamma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
+       vertex* copy(void) const { return new Gamma(*this); }
+       ijpair get_increment(void) const { return ijpair(indices.first+indices.second+1, 0); }
+       ex evaluate(const symbol &x) const;
+private:
+};
+
+ex Gamma::evaluate(const symbol &x) const
+{
+       int i = indices.first;
+       int j = indices.second;
+       return F_ab(1,i,1,j,x);
+}
+
+
+/*
+ * Class of vertices of type Vacuum.
+ */
+class Vacuum : public vertex {
+public:
+       Vacuum(ijpair ij = ijpair(0,0)) : vertex(ij) { }
+       vertex* copy(void) const { return new Vacuum(*this); }
+       ijpair get_increment() const { return ijpair(0, indices.first+indices.second+1); }
+       ex evaluate(const symbol &x) const;
+private:
+};
+
+ex Vacuum::evaluate(const symbol &x) const
+{
+       int i = indices.first;
+       int j = indices.second;
+       return (-TrOne*(F_ab(0,i,1,j,x)-F_ab(1,i,1,j,x)+F_ab(1,i,0,j,x)))/2;
+}
+
+
+/*
+ * Class of nodes (or trees or subtrees), including list of children.
+ */
+class node {
+public:
+       node(const vertex &v) { vert = v.copy(); }
+       node(const node &n) { vert = (n.vert)->copy(); children = n.children; }
+       ~node() { delete vert; }
+       void add_child(const node &, bool = false);
+       ex evaluate(const symbol &x, unsigned grad) const;
+       unsigned total_edges(void) const;
+private:
+       vertex *vert;
+       list<child> children;
+};
+
+void node::add_child(const node &childnode, bool cut)
+{
+       children.push_back(child(childnode, cut));
+       if(!cut)
+               vert->increment_indices(childnode.vert->get_increment());
+}
+
+ex node::evaluate(const symbol &x, unsigned grad) const
+{
+       ex product = 1;
+       for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
+               if (!i->second)
+                       product *= i->first.evaluate(x,grad);
+               else
+                       product *= -div_part(i->first.evaluate(x,grad),x,grad);
+       }
+       return (product * vert->evaluate(x));
+}
+
+unsigned node::total_edges(void) const
+{
+       unsigned accu = 0;
+       for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
+               accu += i->first.total_edges();
+               ++accu;
+       }
+       return accu;
+}
+
+
+/*
+ * These operators let us write down trees in an intuitive way, by adding
+ * arbitrarily complex children to a given vertex.  The eye candy that can be
+ * produced with it makes detection of errors much simpler than with code
+ * written using calls to node's method add_child() because it allows for
+ * editor-assisted indentation.
+ */
+node operator+(const node &n, const child &c)
+{
+       node nn(n);
+       nn.add_child(c.first, c.second);
+       return nn;
+}
+void operator+=(node &n, const child &c)
+{
+       n.add_child(c.first, c.second);
+}
+
+/*
+ * Build this sample rooted tree characterized by a certain combination of
+ * cut or uncut edges as specified by the unsigned parameter:
+ *              Gamma
+ *              /   \
+ *         Sigma     Vacuum
+ *        /   \       /   \
+ *    Sigma Sigma  Sigma0 Sigma
+ */
+static node mytree(unsigned cuts=0)
+{
+       return (Gamma()
+               + child(Sigma()
+                       + child(Sigma(), bool(cuts & 1))
+                       + child(Sigma(), bool(cuts & 2)),
+                       bool(cuts & 4))
+               + child(Vacuum()
+                       + child(Sigma(ijpair(0,0),false), bool(cuts & 8))
+                       + child(Sigma(), bool(cuts & 16)),
+                       bool(cuts & 32)));
+}
+
+
+static unsigned test(void)
+{
+       const symbol x("x");
+       
+       const unsigned edges = mytree().total_edges();
+       const unsigned vertices = edges+1;
+       
+       // fill a vector of all possible 2^edges combinations of cuts...
+       vector<node> counter;
+       for (unsigned i=0; i<(1U<<edges); ++i)
+               counter.push_back(mytree(i));
+       
+       // ...the sum, when evaluated, is the antipode...
+       ex accu = 0;
+       for (vector<node>::iterator i=counter.begin(); i!=counter.end(); ++i)
+               accu += i->evaluate(x,vertices);
+       
+       // ...which is only interesting term-wise in the series expansion...
+       ex result = accu.series(x==0,vertices).expand().normal();
+       
+       // ...and has the nice property that in each term all the Eulers cancel:
+       if (result.has(Euler)) {
+               clog << "The antipode was miscalculated\nAntipode==" << result
+                    << "\nshould not have any occurrence of Euler" << endl;
+               return 1;
+       }
+       return 0;
+}
+
+unsigned time_antipode(void)
+{
+       unsigned result = 0;
+       unsigned count = 0;
+       timer jaeger_le_coultre;
+       double time = .0;
+       
+       cout << "timing computation of an antipode in Yukawa theory" << flush;
+       clog << "-------computation of an antipode in Yukawa theory" << endl;
+       
+       if (do_test) {
+               jaeger_le_coultre.start();
+               // correct for very small times:
+               do {
+                       result = test();
+                       ++count;
+               } while ((time=jaeger_le_coultre.read())<0.1 && !result);
+               cout << '.' << flush;
+               
+               if (!result) {
+                       cout << " passed ";
+                       clog << "(no output)" << endl;
+               } else {
+                       cout << " failed ";
+               }
+               cout << int(1000*(time/count))*0.001 << 's' << endl;
+       } else {
+               cout << " disabled" << endl;
+               clog << "(no output)" << endl;
+       }
+       
+       return result;
+}
index 1e7475c768156d4dd0a4f6854e50d699bceca700..246a51235649053b9279b6c1e983ba36b304b145 100644 (file)
@@ -176,6 +176,13 @@ int main()
                ++result;
        }
        
+       try {
+               result += time_antipode();
+       } catch (const exception &e) {
+               cout << "Error: caught exception " << e.what() << endl;
+               ++result;
+       }
+       
        if (result) {
                cout << "Error: something went wrong. ";
                if (result == 1) {
index 28ce4a758fb8db5bc0a288b08b99eb403691dfe2..9127ca52104edf06eebfaf2b4277979f5bccaad2 100644 (file)
@@ -72,5 +72,6 @@ unsigned time_lw_P();
 unsigned time_lw_Pprime();
 unsigned time_lw_Q();
 unsigned time_lw_Qprime();
+unsigned time_antipode();
 
 #endif // ndef CHECKS_H
index d5196fb224263f63a6b8e68a91d2a6cad9bccb7a..7764bda84d2a4ebe895b344fc748cc3f521729f3 100644 (file)
@@ -46,3 +46,5 @@
 (no output)
 -------Lewis-Wester test Q' (charpoly(P'))
 (no output)
+-------computation of an antipode in Yukawa theory
+(no output)