* 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.
+ * This program is based on work by
+ * Isabella Bierenbaum <bierenbaum@thep.physik.uni-mainz.de> and
+ * Dirk Kreimer <dkreimer@bu.edu>.
+ * For details, please see <http://www.arXiv.org/abs/hep-th/0111192>.
*/
/*
- * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2003 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
#include "times.h"
#include <utility>
-#include <list>
+#include <vector>
+#include <set>
+#include <map>
+#include <typeinfo>
+#include <stdexcept>
// whether to run this beast or not:
static const bool do_test = true;
+// regularization parameter:
+static const symbol x("x");
+
typedef pair<unsigned, unsigned> ijpair;
typedef pair<class node, bool> child;
/* 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);
+ const ex exser = exarg.series(x==0, grad);
+ if (exser.degree(x)<0)
+ throw runtime_error("divergent part truncation disaster");
ex exser_trunc;
for (int i=exser.ldegree(x); i<0; ++i)
exser_trunc += exser.coeff(x,i)*pow(x,i);
virtual ~vertex() { }
virtual vertex* copy(void) const = 0;
virtual ijpair get_increment(void) const { return indices; }
- virtual ex evaluate(const symbol &x) const = 0;
+ virtual const ex evaluate(const symbol &x, const unsigned grad) const = 0;
+ bool operator==(const vertex &v) const { return (indices==v.indices); }
+ bool operator<(const vertex &v) const { return (indices<v.indices); }
protected:
ijpair indices;
};
*/
class Sigma : public vertex {
public:
- Sigma(ijpair ij = ijpair(0,0), bool f = true) : vertex(ij), flag(f) { }
+ Sigma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
vertex* copy(void) const { return new Sigma(*this); }
- ijpair get_increment(void) const;
- ex evaluate(const symbol &x) const;
+ ijpair get_increment(void) const { return ijpair(indices.first+indices.second+1, 0); }
+ const ex evaluate(const symbol &x, const unsigned grad) const;
private:
- bool flag;
};
-ijpair Sigma::get_increment(void) const
+const ex Sigma::evaluate(const symbol &x, const unsigned grad) const
{
- if (flag == true)
- return ijpair(indices.first+1, indices.second);
- else
- return ijpair(indices.first, indices.second+1);
+ // c.f. comments in node::evaluate()
+ static map<Sigma,ex> catalog;
+ static unsigned prev_grad = 0;
+ if (grad>prev_grad) {
+ catalog.clear();
+ prev_grad = grad;
+ }
+ map<Sigma,ex>::iterator pos = catalog.find(*this);
+ if (pos==catalog.end()) {
+ int i = indices.first;
+ int j = indices.second;
+ 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();
+ pos = catalog.insert(map<Sigma,ex>::value_type(*this,result)).first;
+ if (grad<prev_grad)
+ prev_grad = grad;
+ }
+ return pos->second;
}
-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 Sigma_flipped, sitting in the upper fermionline of Vacuum; no consequences for Gamma. */
+class Sigma_flipped : public Sigma {
+public:
+ Sigma_flipped(ijpair ij = ijpair(0,0)) : Sigma(ij) { }
+ vertex* copy(void) const { return new Sigma_flipped(*this); }
+ ijpair get_increment(void) const { return ijpair(0, indices.first+indices.second+1); }
+ const ex evaluate(const symbol &x, const unsigned grad) const { return Sigma::evaluate(x, grad); }
+private:
+};
/*
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;
+ const ex evaluate(const symbol &x, const unsigned grad) const;
private:
};
-ex Gamma::evaluate(const symbol &x) const
+const ex Gamma::evaluate(const symbol &x, const unsigned grad) const
{
- int i = indices.first;
- int j = indices.second;
- return F_ab(1,i,1,j,x);
+ // c.f. comments in node::evaluate()
+ static map<Gamma,ex> catalog;
+ static unsigned prev_grad = 0;
+ if (prev_grad<grad) {
+ catalog.clear();
+ prev_grad = grad;
+ }
+ map<Gamma,ex>::iterator pos = catalog.find(*this);
+ if (pos==catalog.end()) {
+ int i = indices.first;
+ int j = indices.second;
+ const ex result = (F_ab(1,i,1,j,x)).series(x==0,grad).expand();
+ pos = catalog.insert(map<Gamma,ex>::value_type(*this,result)).first;
+ if (grad<prev_grad)
+ prev_grad = grad;
+ }
+ return pos->second;
}
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;
+ const ex evaluate(const symbol &x, const unsigned grad) const;
private:
};
-ex Vacuum::evaluate(const symbol &x) const
+const ex Vacuum::evaluate(const symbol &x, const unsigned grad) 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;
+ // c.f. comments in node::evaluate()
+ static map<Vacuum,ex> catalog;
+ static unsigned prev_grad = 0;
+ if (prev_grad<grad) {
+ catalog.clear();
+ prev_grad = grad;
+ }
+ map<Vacuum,ex>::iterator pos = catalog.find(*this);
+ if (pos==catalog.end()) {
+ int i = indices.first;
+ int j = indices.second;
+ 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();
+ pos = catalog.insert(map<Vacuum,ex>::value_type(*this,result)).first;
+ if (grad<prev_grad)
+ prev_grad = grad;
+ }
+ return pos->second;
}
public:
node(const vertex &v) { vert = v.copy(); }
node(const node &n) { vert = (n.vert)->copy(); children = n.children; }
+ const node & operator=(const node &);
~node() { delete vert; }
void add_child(const node &, bool = false);
- ex evaluate(const symbol &x, unsigned grad) const;
+ const ex evaluate(const symbol &x, unsigned grad) const;
unsigned total_edges(void) const;
+ bool operator==(const node &) const;
+ bool operator<(const node &) const;
private:
vertex *vert;
- list<child> children;
+ multiset<child> children;
};
+const node & node::operator=(const node &n)
+{
+ if (this!=&n) {
+ delete vert;
+ vert = (n.vert)->copy();
+ children = n.children;
+ }
+ return *this;
+}
+
void node::add_child(const node &childnode, bool cut)
{
- children.push_back(child(childnode, cut));
+ children.insert(child(childnode, cut));
if(!cut)
vert->increment_indices(childnode.vert->get_increment());
}
-ex node::evaluate(const symbol &x, unsigned grad) const
+const ex node::evaluate(const symbol &x, unsigned grad) const
{
- ex product = 1;
- for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
+ static map<node,ex> catalog; // lookup table for already evaluated subnodes
+ static unsigned prev_grad = 0; // grad argument that the catalog has been build for
+ if (grad>prev_grad) {
+ // We haven't computed far enough last time. Our catalog cannot cope with
+ // the demands for this value of grad so let's clear it.
+ catalog.clear();
+ prev_grad = grad;
+ }
+ ex product = 1; // accumulator for all the children
+ for (multiset<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
+ map<node,ex>::iterator pos = catalog.find(i->first);
+ if (pos==catalog.end()) {
+ pos = catalog.insert(map<node,ex>::value_type(i->first,i->first.evaluate(x,grad).series(x==0,grad).expand())).first;
+ if (grad<prev_grad) {
+ // We have just spoiled the catalog by inserting a series computed
+ // to lower grad than the others in it. So let's make sure next time
+ // we don't use one of the newly inserted ones by bumping prev_grad
+ // down to the current value of grad.
+ prev_grad = grad;
+ }
+ }
if (!i->second)
- product *= i->first.evaluate(x,grad);
+ product *= pos->second;
else
- product *= -div_part(i->first.evaluate(x,grad),x,grad);
+ product *= -div_part(pos->second,x,grad);
}
- return (product * vert->evaluate(x));
+ return (product * vert->evaluate(x,grad));
}
unsigned node::total_edges(void) const
{
unsigned accu = 0;
- for (list<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
+ for (multiset<child>::const_iterator i=children.begin(); i!=children.end(); ++i) {
accu += i->first.total_edges();
++accu;
}
return accu;
}
+bool node::operator==(const node &n) const
+{
+ // Are the types of the top-level node vertices equal?
+ if (typeid(*vert)!=typeid(*n.vert))
+ return false;
+ // Are the indices of the top-level nodes equal?
+ if (!(*vert==*n.vert))
+ return false;
+ // Are the sets of children equal, one by one?
+ return (children==n.children);
+}
+
+bool node::operator<(const node &n) const
+{
+ // Are the types of the top-level node vertices different?
+ if (typeid(*vert)!=typeid(*n.vert))
+ return typeid(*vert).before(typeid(*n.vert));
+ // Are the indices of the top-level nodes different?
+ if (!(*vert==*n.vert))
+ return (vert<n.vert);
+ // Are the sets of children different, one by one?
+ return (children<n.children);
+}
/*
* These operators let us write down trees in an intuitive way, by adding
* written using calls to node's method add_child() because it allows for
* editor-assisted indentation.
*/
-node operator+(const node &n, const child &c)
+const 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
+
+/* Gamma
+ * |
+ * Gamma
*/
-static node mytree(unsigned cuts=0)
+static const node tree1(unsigned cuts=0)
+{
+ return (Gamma()
+ + child(Gamma(),
+ bool(cuts & 1)));
+}
+
+/* Gamma
+ * / | \
+ * Vacuum Gamma Vacuum
+ * / | \
+ * Sigma Sigma Sigma0
+ */
+static const node tree2(unsigned cuts=0)
+{
+ return (Gamma()
+ + child(Vacuum()
+ + child(Sigma(), bool(cuts & 1))
+ + child(Sigma(), bool(cuts & 2))
+ + child(Sigma_flipped(), bool(cuts & 4)),
+ bool(cuts & 8))
+ + child(Gamma(), bool(cuts & 16))
+ + child(Gamma(), bool(cuts & 32)));
+}
+
+/* Gamma
+ * |
+ * Gamma
+ * |
+ * Gamma
+ * / \
+ * Vacuum Gamma
+ * / \ \
+ * Sigma Sigma Sigma
+ */
+static const node tree3(unsigned cuts=0)
+{
+ return (Gamma()
+ + child(Gamma()
+ + child(Gamma()
+ + child(Gamma()
+ + child(Sigma(), bool(cuts & 1)),
+ bool(cuts & 2))
+ + child(Vacuum()
+ + child(Sigma(), bool(cuts & 4))
+ + child(Sigma(), bool(cuts & 8)),
+ bool(cuts & 16)),
+ bool(cuts & 32)),
+ bool(cuts & 64)));
+}
+
+/* Gamma
+ * / \
+ * Sigma Vacuum
+ * / \ / \
+ * Sigma Sigma Sigma0 Sigma
+ */
+static const node tree4(unsigned cuts=0)
{
return (Gamma()
+ child(Sigma()
+ child(Sigma(), bool(cuts & 2)),
bool(cuts & 4))
+ child(Vacuum()
- + child(Sigma(ijpair(0,0),false), bool(cuts & 8))
+ + child(Sigma_flipped(), bool(cuts & 8))
+ + child(Sigma(), bool(cuts & 16)),
+ bool(cuts & 32)));
+}
+
+/* Sigma
+ * / | \
+ * Sigma Vacuum Vacuum
+ * / \ \
+ * Sigma Sigma0 Sigma
+ */
+static const node tree5(unsigned cuts=0)
+{
+ return (Sigma()
+ + child(Sigma(), bool(cuts & 1))
+ + child(Vacuum()
+ + child(Sigma(), bool(cuts & 2))
+ + child(Sigma_flipped(), bool(cuts & 4)),
+ bool(cuts & 8))
+ + child(Vacuum()
+ child(Sigma(), bool(cuts & 16)),
bool(cuts & 32)));
}
+/* Vacuum
+ * / \
+ * Sigma Sigma0
+ * | |
+ * Sigma Sigma
+ * |
+ * Vacuum
+ */
+static const node tree6(unsigned cuts=0)
+{
+ return (Vacuum()
+ + child(Sigma()
+ + child(Sigma(), bool(cuts & 1)),
+ bool(cuts & 2))
+ + child(Sigma_flipped()
+ + child(Sigma()
+ + child(Vacuum(), bool(cuts & 4)),
+ bool(cuts & 8)),
+ bool(cuts & 16)));
+}
-static unsigned test(void)
+static unsigned test_tree(const node (*tree_generator)(unsigned=0))
{
- const symbol x("x");
-
- const unsigned edges = mytree().total_edges();
- const unsigned vertices = edges+1;
+ const int edges = tree_generator().total_edges();
+ const int 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));
+ counter.push_back(tree_generator(i));
- // ...the sum, when evaluated, is the antipode...
- ex accu = 0;
+ // ...the sum, when evaluated and reexpanded, is the antipode...
+ ex result = 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();
+ result = (result+i->evaluate(x,vertices)).series(x==0,vertices).expand();
// ...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;
+ } else if (result.ldegree(x)!=-vertices || result.degree(x)!=0) {
+ clog << "The antipode was miscalculated\nAntipode==" << result
+ << "\nshould run from " << x << "^(" << -vertices << ") to "
+ << x << "^(0)" << "but it runs from " << x << "^("
+ << result.ldegree(x) << ")" << "to " << x << "^("
+ << result.degree(x) << ")" << 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;
+ cout << "timing computation of antipodes in Yukawa theory" << flush;
+ clog << "-------computation of antipodes 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;
+ result += test_tree(tree1); cout << '.' << flush;
+ result += test_tree(tree2); cout << '.' << flush;
+ result += test_tree(tree3); cout << '.' << flush;
+ result += test_tree(tree4); cout << '.' << flush;
+ result += test_tree(tree5); cout << '.' << flush;
+ result += test_tree(tree6); cout << '.' << flush;
if (!result) {
cout << " passed ";
} else {
cout << " failed ";
}
- cout << int(1000*(time/count))*0.001 << 's' << endl;
+ cout << int(1000*jaeger_le_coultre.read())*0.001 << "s (total)" << endl;
} else {
cout << " disabled" << endl;
clog << "(no output)" << endl;
git://www.ginac.de / ginac.git / blobdiff