From: Richard Kreckel Date: Mon, 5 Mar 2001 20:02:00 +0000 (+0000) Subject: * New funny timing added: compute an antipode in Yukawa theory. X-Git-Tag: release_0-8-0~36 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=dad107ff48f68d45e72469a8716df375ae145cf3 * New funny timing added: compute an antipode in Yukawa theory. --- diff --git a/check/Makefile.am b/check/Makefile.am index 7932f1c6..a5e781b9 100644 --- a/check/Makefile.am +++ b/check/Makefile.am @@ -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 index 00000000..2d3a13ab --- /dev/null +++ b/check/time_antipode.cpp @@ -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 +#include + +// whether to run this beast or not: +static const bool do_test = true; + +typedef pair ijpair; +typedef pair 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 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::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::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 counter; + for (unsigned i=0; i<(1U<::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; +} diff --git a/check/times.cpp b/check/times.cpp index 1e7475c7..246a5123 100644 --- a/check/times.cpp +++ b/check/times.cpp @@ -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) { diff --git a/check/times.h b/check/times.h index 28ce4a75..9127ca52 100644 --- a/check/times.h +++ b/check/times.h @@ -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 diff --git a/check/times.ref b/check/times.ref index d5196fb2..7764bda8 100644 --- a/check/times.ref +++ b/check/times.ref @@ -46,3 +46,5 @@ (no output) -------Lewis-Wester test Q' (charpoly(P')) (no output) +-------computation of an antipode in Yukawa theory +(no output)