X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=check%2Ftime_antipode.cpp;h=43d48f64d85c03b6ccdd0fff0f65c1da0c45f2d5;hp=2d3a13ab706d594228b1a2a792f18d5ac630c806;hb=8247d44aac20d266d1b135468dacf521ea655f20;hpb=dad107ff48f68d45e72469a8716df375ae145cf3

diff --git a/check/time_antipode.cpp b/check/time_antipode.cpp
index 2d3a13ab..43d48f64 100644
--- a/check/time_antipode.cpp
+++ b/check/time_antipode.cpp
@@ -9,12 +9,14 @@
  *  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
@@ -33,11 +35,18 @@
 
 #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;
 
@@ -46,13 +55,9 @@ 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);
+	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);
@@ -79,9 +84,11 @@ 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; 
+	virtual vertex* copy() const = 0;
+	virtual ijpair get_increment() const { return indices; }
+	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;
 };
@@ -92,28 +99,44 @@ protected:
  */
 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;
+	Sigma(ijpair ij = ijpair(0,0)) : vertex(ij) { }
+	vertex* copy() const { return new Sigma(*this); }
+	ijpair get_increment() 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() const { return new Sigma_flipped(*this); }
+	ijpair get_increment() 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:
+};
 
 
 /*
@@ -122,17 +145,31 @@ ex Sigma::evaluate(const symbol &x) const
 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;
+	vertex* copy() const { return new Gamma(*this); }
+  	ijpair get_increment() const { return ijpair(indices.first+indices.second+1, 0); }
+	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;
 }
 
 
@@ -142,17 +179,31 @@ ex Gamma::evaluate(const symbol &x) const
 class Vacuum : public vertex {
 public:
 	Vacuum(ijpair ij = ijpair(0,0)) : vertex(ij) { }
-	vertex* copy(void) const { return new Vacuum(*this); }
+	vertex* copy() 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;
 }
 
 
@@ -163,44 +214,99 @@ class node {
 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;
-	unsigned total_edges(void) const;
+	const ex evaluate(const symbol &x, unsigned grad) const;
+	unsigned total_edges() 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 node::total_edges() 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
@@ -209,27 +315,81 @@ unsigned node::total_edges(void) const
  * 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 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 node mytree(unsigned cuts=0)
+static const node tree4(unsigned cuts=0)
 {
 	return (Gamma()
 	        + child(Sigma()
@@ -237,59 +397,98 @@ static node mytree(unsigned cuts=0)
 	                + 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)));
+}
 
-static unsigned test(void)
+/*               Vacuum
+ *               /    \
+ *           Sigma    Sigma0
+ *             |        |
+ *           Sigma    Sigma
+ *                      |
+ *                    Vacuum
+ */
+static const node tree6(unsigned cuts=0)
 {
-	const symbol x("x");
-	
-	const unsigned edges = mytree().total_edges();
-   	const unsigned vertices = edges+1;
+	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_tree(const node tree_generator(unsigned))
+{
+	const int edges = tree_generator(0).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 time_antipode()
 {
 	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 ";
@@ -297,7 +496,7 @@ unsigned time_antipode(void)
 		} 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;