From: Richard Kreckel Date: Tue, 14 Aug 2001 18:18:10 +0000 (+0000) Subject: - Updated to the optimized version. This allows us to X-Git-Tag: release_0-9-3~12 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=8710f8123e5f137fe32a34b5c4de9f5aca6e8d73;hp=4e3f330d382e88e2ae180adcaea71ba5067cbf50 - Updated to the optimized version. This allows us to a) add more tests for consistency b) actually believe the result of the timings because they were irreproducible before (due to different branches taken by different representations, which happens because eariler tests account for low timings by iterating an unknown number of times thus producing different hashvalues for symbols in later calculations.) --- diff --git a/check/time_antipode.cpp b/check/time_antipode.cpp index f44c0b98..2d64beef 100644 --- a/check/time_antipode.cpp +++ b/check/time_antipode.cpp @@ -13,8 +13,6 @@ * Isabella Bierenbaum and * Dirk Kreimer . * For details, please ask for the diploma theses of Isabella Bierenbaum. - * Also, this file is based on an early version of their program, they now - * have production-code that is somewhat more efficient than the stuff below. */ /* @@ -38,11 +36,17 @@ #include "times.h" #include #include +#include +#include +#include #include // whether to run this beast or not: static const bool do_test = true; +// regularization parameter: +static const symbol x("x"); + typedef pair ijpair; typedef pair child; @@ -82,7 +86,9 @@ public: 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 catalog; + static unsigned prev_grad = 0; + if (grad>prev_grad) { + catalog.clear(); + prev_grad = grad; + } + map::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::value_type(*this,result)).first; + if (gradsecond; } -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: +}; /* @@ -125,15 +147,29 @@ 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; + 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 catalog; + static unsigned prev_grad = 0; + if (prev_grad::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::value_type(*this,result)).first; + if (gradsecond; } @@ -145,15 +181,29 @@ 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; + 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 catalog; + static unsigned prev_grad = 0; + if (prev_grad::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::value_type(*this,result)).first; + if (gradsecond; } @@ -167,50 +217,96 @@ public: 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; - vector children; + multiset children; }; const node & node::operator=(const node &n) { - delete vert; - vert = (n.vert)->copy(); - children = n.children; + 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 (vector::const_iterator i=children.begin(); i!=children.end(); ++i) { + static map 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::const_iterator i=children.begin(); i!=children.end(); ++i) { + map::iterator pos = catalog.find(i->first); + if (pos==catalog.end()) { + pos = catalog.insert(map::value_type(i->first,i->first.evaluate(x,grad).series(x==0,grad).expand())).first; + if (gradsecond) - 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 (vector::const_iterator i=children.begin(); i!=children.end(); ++i) { + for (multiset::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 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(); + 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; } @@ -285,21 +474,19 @@ static unsigned test(void) 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 "; @@ -307,7 +494,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; diff --git a/check/times.ref b/check/times.ref index 7764bda8..08d2b9c4 100644 --- a/check/times.ref +++ b/check/times.ref @@ -46,5 +46,5 @@ (no output) -------Lewis-Wester test Q' (charpoly(P')) (no output) --------computation of an antipode in Yukawa theory +-------computation of antipodes in Yukawa theory (no output)