--- a/examples/cohomology/rips-pairwise-cohomology.cpp Thu May 14 14:04:43 2009 -0700
+++ b/examples/cohomology/rips-pairwise-cohomology.cpp Thu Jul 09 00:59:32 2009 -0700
@@ -6,6 +6,7 @@
#include <utilities/containers.h> // for BackInsertFunctor
#include <utilities/timer.h>
+#include <utilities/log.h>
#include <string>
@@ -15,28 +16,36 @@
typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
typedef PairDistances::DistanceType DistanceType;
typedef PairDistances::IndexType Vertex;
-
-typedef CohomologyPersistence<DistanceType> Persistence;
-typedef Persistence::SimplexIndex Index;
-typedef Persistence::Death Death;
-
+
typedef Rips<PairDistances> Generator;
typedef Generator::Simplex Smplx;
+typedef std::vector<Smplx> SimplexVector;
+typedef SimplexVector::const_iterator SV_const_iterator;
-typedef std::map<Smplx, Index,
- Smplx::VertexComparison> Complex;
-typedef std::vector<Smplx> SimplexVector;
+typedef boost::tuple<Dimension, DistanceType> BirthInfo;
+typedef CohomologyPersistence<BirthInfo, SV_const_iterator> Persistence;
+typedef Persistence::SimplexIndex Index;
+typedef Persistence::Death Death;
+typedef std::map<Smplx, Index, Smplx::VertexComparison> Complex;
-void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance);
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime);
int main(int argc, char* argv[])
{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
Dimension skeleton;
DistanceType max_distance;
+ ZpField::Element prime;
std::string infilename;
- program_options(argc, argv, infilename, skeleton, max_distance);
+ program_options(argc, argv, infilename, skeleton, max_distance, prime);
+ Timer total_timer; total_timer.start();
PointContainer points;
read_points(infilename, points);
@@ -46,30 +55,62 @@
SimplexVector v;
Complex c;
+ Timer rips_timer; rips_timer.start();
rips.generate(skeleton, max_distance, make_push_back_functor(v));
std::sort(v.begin(), v.end(), Generator::Comparison(distances));
+ rips_timer.stop();
std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
Timer persistence_timer; persistence_timer.start();
- Persistence p;
+ ZpField zp(prime);
+ Persistence p(zp);
for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
{
std::vector<Index> boundary;
- for (Smplx::BoundaryIterator bcur = cur->boundary_begin();
- bcur != cur->boundary_end(); ++bcur)
+ for (Smplx::BoundaryIterator bcur = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
boundary.push_back(c[*bcur]);
Index idx; Death d;
- boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), size(*cur));
- c[*cur] = idx;
- if (d && (size(*cur) - *d) > 0)
- std::cout << (cur->dimension() - 1) << " " << *d << " " << size(*cur) << std::endl;
+ bool store = cur->dimension() < skeleton;
+ boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, cur);
+
+ // c[*cur] = idx;
+ if (store)
+ c.insert(std::make_pair(*cur, idx));
+
+ if (d && (size(*cur) - d->get<1>()) > 0)
+ {
+ AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
+ std::cout << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
+ }
}
+ // output infinte persistence cocycles
+ for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
+ std::cout << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
persistence_timer.stop();
+
+
+ // p.show_cocycles();
+ // Output alive cocycles
+ for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
+ {
+ std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
+ for (Persistence::ZColumn::const_iterator zcur = cur->zcolumn.begin(); zcur != cur->zcolumn.end(); ++zcur)
+ {
+ const Smplx& s = **(zcur->si);
+ std::cout << zcur->coefficient << " ";
+ for (Smplx::VertexContainer::const_iterator vcur = s.vertices().begin(); vcur != s.vertices().end(); ++vcur)
+ std::cout << *vcur << " ";
+ std::cout << std::endl;
+ }
+ }
+ total_timer.stop();
+ rips_timer.check("Rips timer");
persistence_timer.check("Persistence timer");
+ total_timer.check("Total timer");
}
-void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance)
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime)
{
namespace po = boost::program_options;
@@ -81,7 +122,13 @@
visible.add_options()
("help,h", "produce help message")
("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2), "Dimension of the Rips complex we want to compute")
+ ("prime,p", po::value<ZpField::Element>(&prime)->default_value(2), "Prime p for the field F_p")
("max-distance,m", po::value<DistanceType>(&max_distance)->default_value(Infinity), "Maximum value for the Rips complex construction");
+#if LOGGING
+ std::vector<std::string> log_channels;
+ visible.add_options()
+ ("log,l", po::value< std::vector<std::string> >(&log_channels), "log channels to turn on (info, debug, etc)");
+#endif
po::positional_options_description pos;
pos.add("input-file", 1);
@@ -94,6 +141,11 @@
options(all).positional(pos).run(), vm);
po::notify(vm);
+#if LOGGING
+ for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
+ stderrLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
+#endif
+
if (vm.count("help") || !vm.count("input-file"))
{
std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl;
--- a/examples/cohomology/triangle-cohomology.cpp Thu May 14 14:04:43 2009 -0700
+++ b/examples/cohomology/triangle-cohomology.cpp Thu Jul 09 00:59:32 2009 -0700
@@ -54,7 +54,7 @@
stderrLog.subscribeTo(RLOG_CHANNEL("info"));
stderrLog.subscribeTo(RLOG_CHANNEL("error"));
- //stderrLog.subscribeTo(RLOG_CHANNEL("topology"));
+ stderrLog.subscribeTo(RLOG_CHANNEL("topology"));
#endif
SimplexVector v;
@@ -66,10 +66,13 @@
// Compute persistence
Complex c;
- Persistence p;
+ ZpField zp(11);
+ Persistence p(zp);
unsigned i = 0;
for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
{
+ std::cout << "-------" << std::endl;
+
std::vector<Index> boundary;
for (Smplx::BoundaryIterator bcur = cur->boundary_begin();
bcur != cur->boundary_end(); ++bcur)
@@ -84,7 +87,7 @@
// (i.e. when a 1 class kills 0, it's really that in cohomology forward 0 kills 1,
// in cohomology backward 1 kills 0, and in homology 1 kills 0)
- //p.show_cycles();
+ p.show_cocycles();
}
}
--- a/include/topology/cohomology-persistence.h Thu May 14 14:04:43 2009 -0700
+++ b/include/topology/cohomology-persistence.h Thu Jul 09 00:59:32 2009 -0700
@@ -16,6 +16,7 @@
#include <list>
#include <utility>
+#include <topology/field-arithmetic.h>
#include "utilities/types.h"
#include <boost/optional.hpp>
@@ -23,86 +24,105 @@
namespace bi = boost::intrusive;
-template<class BirthInfo_, class SimplexData_ = Empty<> >
+template<class BirthInfo_, class SimplexData_ = Empty<>, class Field_ = ZpField>
class CohomologyPersistence
{
public:
typedef BirthInfo_ BirthInfo;
typedef SimplexData_ SimplexData;
+ typedef Field_ Field;
+
+ typedef typename Field::Element FieldElement;
- struct SNode;
- typedef bi::list<SNode, bi::constant_time_size<false> > ZRow;
+ CohomologyPersistence(const Field& field = Field()):
+ field_(field) {}
+
- // Simplex representation
- struct SHead: public SimplexData
- {
- SHead(const SHead& other):
- SimplexData(other), order(other.order) {} // don't copy row since we can't
- SHead(const SimplexData& sd, unsigned o):
- SimplexData(sd), order(o) {}
+ // An entry in a cocycle column
+ struct SNode; // members: si, coefficient, ci
+ typedef s::vector<SNode> ZColumn;
+ typedef bi::list<SNode, bi::constant_time_size<false> > ZRow;
+ class CompareSNode;
- // intrusive list corresponding to row of s in Z^*, not ordered in any particular order
- ZRow row;
- unsigned order;
- };
-
+ struct SHead; // members: row, order
typedef s::list<SHead> Simplices;
typedef typename Simplices::iterator SimplexIndex;
- struct Cocycle;
+ struct Cocycle; // members: zcolumn, birth, order
typedef s::list<Cocycle> Cocycles;
typedef typename Cocycles::iterator CocycleIndex;
-
- // An entry in a cocycle column; it's also an element in an intrusive list, hence the list_base_hook<>
- typedef bi::list_base_hook<bi::link_mode<bi::auto_unlink> > auto_unlink_hook;
- struct SNode: public auto_unlink_hook
- {
- SNode() {}
- SNode(SimplexIndex sidx): si(sidx) {}
-
- // eventually store a field element
-
- SimplexIndex si;
- CocycleIndex cocycle; // TODO: is there no way to get rid of this overhead?
-
- void unlink() { auto_unlink_hook::unlink(); }
- };
- class CompareSNode;
+ typedef std::pair<CocycleIndex, FieldElement> CocycleCoefficientPair;
- typedef s::vector<SNode> ZColumn;
- struct Cocycle
- {
- Cocycle(const BirthInfo& b, unsigned o):
- birth(b), order(o) {}
-
- ZColumn cocycle;
- BirthInfo birth;
- unsigned order;
-
- bool operator<(const Cocycle& other) const { return order > other.order; }
- bool operator==(const Cocycle& other) const { return order == other.order; }
- };
-
- typedef boost::optional<BirthInfo> Death;
- typedef std::pair<SimplexIndex,
- Death> IndexDeathPair;
+ typedef boost::optional<BirthInfo> Death;
+ typedef std::pair<SimplexIndex, Death> IndexDeathPair;
// return either a SimplexIndex or a Death
// BI = BoundaryIterator; it should dereference to a SimplexIndex
template<class BI>
- IndexDeathPair add(BI begin, BI end, BirthInfo b, const SimplexData& sd = SimplexData());
+ IndexDeathPair add(BI begin, BI end, BirthInfo b, bool store = true, const SimplexData& sd = SimplexData());
- void show_cycles() const;
-
+ void show_cocycles() const;
+ CocycleIndex begin() { return cocycles_.begin(); }
+ CocycleIndex end() { return cocycles_.end(); }
private:
- void add_cocycle(Cocycle& z1, Cocycle& z2);
+ void add_cocycle(CocycleCoefficientPair& z1, CocycleCoefficientPair& z2);
private:
Simplices simplices_;
Cocycles cocycles_;
+ Field field_;
};
+
+// Simplex representation
+template<class BirthInfo_, class SimplexData_, class Field_>
+struct CohomologyPersistence<BirthInfo_, SimplexData_, Field_>::SHead: public SimplexData
+{
+ SHead(const SHead& other):
+ SimplexData(other), order(other.order) {} // don't copy row since we can't
+ SHead(const SimplexData& sd, unsigned o):
+ SimplexData(sd), order(o) {}
+
+ // intrusive list corresponding to row of s in Z^*, not ordered in any particular order
+ ZRow row;
+ unsigned order;
+};
+
+// An entry in a cocycle column; it's also an element in an intrusive list, hence the list_base_hook<>
+typedef bi::list_base_hook<bi::link_mode<bi::auto_unlink> > auto_unlink_hook;
+template<class BirthInfo_, class SimplexData_, class Field_>
+struct CohomologyPersistence<BirthInfo_, SimplexData_, Field_>::SNode: public auto_unlink_hook
+{
+ SNode(const SNode& other):
+ si(other.si), coefficient(other.coefficient),
+ ci(other.ci) {}
+
+ SNode(SimplexIndex sidx, FieldElement coef, CocycleIndex cidx):
+ si(sidx), coefficient(coef), ci(cidx) {}
+
+ SimplexIndex si;
+ FieldElement coefficient;
+
+ CocycleIndex ci; // TODO: is there no way to get rid of this overhead?
+
+ void unlink() { auto_unlink_hook::unlink(); }
+};
+
+template<class BirthInfo_, class SimplexData_, class Field_>
+struct CohomologyPersistence<BirthInfo_, SimplexData_, Field_>::Cocycle
+{
+ Cocycle(const BirthInfo& b, unsigned o):
+ birth(b), order(o) {}
+
+ ZColumn zcolumn;
+ BirthInfo birth;
+ unsigned order;
+
+ bool operator<(const Cocycle& other) const { return order > other.order; }
+ bool operator==(const Cocycle& other) const { return order == other.order; }
+};
+
#include "cohomology-persistence.hpp"
--- a/include/topology/cohomology-persistence.hpp Thu May 14 14:04:43 2009 -0700
+++ b/include/topology/cohomology-persistence.hpp Thu Jul 09 00:59:32 2009 -0700
@@ -9,52 +9,62 @@
static rlog::RLogChannel* rlCohomology = DEF_CHANNEL("topology/cohomology", rlog::Log_Debug);
#endif
-template<class BirthInfo, class SimplexData>
-class CohomologyPersistence<BirthInfo, SimplexData>::CompareSNode
+template<class BirthInfo, class SimplexData, class Field>
+class CohomologyPersistence<BirthInfo, SimplexData, Field>::CompareSNode
{
public:
bool operator()(const SNode& s1, const SNode& s2) const { return s1.si->order < s2.si->order; }
};
-template<class BirthInfo, class SimplexData>
+template<class BirthInfo, class SimplexData, class Field>
template<class BI>
-typename CohomologyPersistence<BirthInfo, SimplexData>::IndexDeathPair
-CohomologyPersistence<BirthInfo, SimplexData>::
-add(BI begin, BI end, BirthInfo birth, const SimplexData& sd)
+typename CohomologyPersistence<BirthInfo, SimplexData, Field>::IndexDeathPair
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+add(BI begin, BI end, BirthInfo birth, bool store, const SimplexData& sd)
{
// Create simplex representation
simplices_.push_back(SHead(sd, simplices_.empty() ? 0 : (simplices_.back().order + 1)));
SimplexIndex si = boost::prior(simplices_.end());
// Find out if there are cocycles that evaluate to non-zero on the new simplex
- typedef std::list<CocycleIndex> Candidates;
+ typedef std::list<CocycleCoefficientPair> Candidates;
Candidates candidates, candidates_bulk;
rLog(rlCohomology, "Boundary");
+ bool sign = true; // TODO: this is very crude, fix later
for (BI cur = begin; cur != end; ++cur)
{
- rLog(rlCohomology, " %d", (*cur)->order);
+ rLog(rlCohomology, " %d %d", (*cur)->order, sign);
for (typename ZRow::const_iterator zcur = (*cur)->row.begin(); zcur != (*cur)->row.end(); ++zcur)
- candidates_bulk.push_back(zcur->cocycle);
+ candidates_bulk.push_back(std::make_pair(zcur->ci,
+ (sign ? zcur->coefficient : field_.neg(zcur->coefficient))));
+ sign = !sign;
}
- candidates_bulk.sort(make_indirect_comparison(std::less<Cocycle>()));
+ candidates_bulk.sort(make_first_comparison(make_indirect_comparison(std::less<Cocycle>())));
+#if LOGGING
rLog(rlCohomology, " Candidates bulk");
for (typename Candidates::iterator cur = candidates_bulk.begin();
cur != candidates_bulk.end(); ++cur)
- rLog(rlCohomology, " %d", (*cur)->order);
+ rLog(rlCohomology, " %d %d", cur->first->order, cur->second);
+#endif
- // Remove duplicates --- this is really Z_2, we need a more sophisticated
+ // Remove duplicates
{
typename Candidates::const_iterator cur = candidates_bulk.begin();
while (cur != candidates_bulk.end())
{
typename Candidates::const_iterator next = cur;
- unsigned count = 0;
- while (next != candidates_bulk.end() && *next == *cur) { ++next; ++count; }
+ FieldElement sum = field_.zero();
+ while (next != candidates_bulk.end() && next->first == cur->first)
+ {
+ sum = field_.add(sum, next->second);
+ ++next;
+ }
- if (count % 2)
- candidates.push_back(*cur);
+ rLog(rlCohomology, " Bulk candidate %d sum %d", cur->first->order, sum);
+ if (!field_.is_zero(sum))
+ candidates.push_back(CocycleCoefficientPair(cur->first, sum));
cur = next;
}
@@ -63,16 +73,22 @@
// Birth
if (candidates.empty())
{
- rLog(rlCohomology, "Birth");
+ if (!store)
+ {
+ simplices_.pop_back();
+ return std::make_pair(simplices_.begin(), Death()); // TODO: shouldn't return front
+ }
unsigned order = cocycles_.empty() ? 0 : cocycles_.front().order + 1;
cocycles_.push_front(Cocycle(birth, order));
+
+ rLog(rlCohomology, "Birth: %d", cocycles_.front().order);
// set up the cocycle
- ZColumn& cocycle = cocycles_.front().cocycle;
- cocycle.push_back(si);
- cocycle.front().cocycle = cocycles_.begin();
- si->row.push_back(cocycles_.front().cocycle.front());
+ ZColumn& cocycle = cocycles_.front().zcolumn;
+ cocycle.push_back(SNode(si, field_.id(), cocycles_.begin()));
+ si->row.push_back(cocycles_.front().zcolumn.front());
+ rLog(rlCohomology, " Cocyle: %d", si->order);
return std::make_pair(si, Death());
}
@@ -80,67 +96,97 @@
// Death
rLog(rlCohomology, "Death");
-#if 0
+#if LOGGING
// Debug only, output candidates
rLog(rlCohomology, " Candidates");
- for (typename Candidates::iterator cur = candidates.begin();
- cur != candidates.end(); ++cur)
- rLog(rlCohomology, " %d", (*cur)->order);
+ for (typename Candidates::iterator cur = candidates.begin(); cur != candidates.end(); ++cur)
+ rLog(rlCohomology, " %d %d", cur->first->order, cur->second);
#endif
- Cocycle& z = *candidates.front();
- Death d = z.birth;
+ CocycleCoefficientPair& z = candidates.front();
+ Death d = z.first->birth;
// add z to everything else in candidates
for (typename Candidates::iterator cur = boost::next(candidates.begin());
cur != candidates.end(); ++cur)
- add_cocycle(**cur, z);
+ add_cocycle(*cur, z);
- for (typename ZColumn::iterator cur = z.cocycle.begin(); cur != z.cocycle.end(); ++cur)
+ for (typename ZColumn::iterator cur = z.first->zcolumn.begin(); cur != z.first->zcolumn.end(); ++cur)
cur->unlink();
- cocycles_.erase(candidates.front());
+ cocycles_.erase(candidates.front().first);
return std::make_pair(si, d);
}
-template<class BirthInfo, class SimplexData>
+template<class BirthInfo, class SimplexData, class Field>
void
-CohomologyPersistence<BirthInfo, SimplexData>::
-show_cycles() const
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+show_cocycles() const
{
- std::cout << "Cocycles" << std::endl;
+ std::cout << "Cocycles: " << cocycles_.size() << std::endl;
for (typename Cocycles::const_iterator cur = cocycles_.begin(); cur != cocycles_.end(); ++cur)
{
- std::cout << cur->order << ": ";
- for (typename ZColumn::const_iterator zcur = cur->cocycle.begin(); zcur != cur->cocycle.end(); ++zcur)
- std::cout << zcur->si->order << ", ";
+ // std::cout << cur->order << " (" << cur->birth << "): ";
+ for (typename ZColumn::const_iterator zcur = cur->zcolumn.begin(); zcur != cur->zcolumn.end(); ++zcur)
+ std::cout << zcur->coefficient << " * " << zcur->si->order << ", ";
std::cout << std::endl;
}
}
-template<class BirthInfo, class SimplexData>
+template<class BirthInfo, class SimplexData, class Field>
void
-CohomologyPersistence<BirthInfo, SimplexData>::
-add_cocycle(Cocycle& to, Cocycle& from)
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+add_cocycle(CocycleCoefficientPair& to, CocycleCoefficientPair& from)
{
- rLog(rlCohomology, "Adding cocycle %d to %d", from.order, to.order);
+ rLog(rlCohomology, "Adding cocycle %d to %d", from.first->order, to.first->order);
+
+ ZColumn nw;
+ FieldElement multiplier = field_.neg(field_.div(to.second, from.second));
+ CocycleIndex ci = to.first;
- ZColumn nw;
+ typename ZColumn::iterator tcur = to.first->zcolumn.begin();
+ typename ZColumn::iterator fcur = from.first->zcolumn.begin();
+ CompareSNode cmp;
+ while (tcur != to.first->zcolumn.end() && fcur != from.first->zcolumn.end())
+ {
+ rLog(rlCohomology, " %d %d", tcur->si->order, fcur->si->order);
+ if (cmp(*tcur, *fcur))
+ {
+ nw.push_back(*tcur);
+ ++tcur;
+ }
+ else if (cmp(*fcur, *tcur))
+ {
+ nw.push_back(SNode(fcur->si, field_.mul(multiplier, fcur->coefficient), ci));
+ ++fcur;
+ }
+ else // equality
+ {
+ FieldElement res = field_.mul(multiplier, tcur->coefficient);
+ res = field_.add(fcur->coefficient, res);
+ if (!field_.is_zero(res))
+ nw.push_back(SNode(fcur->si, res, ci));
+ ++tcur; ++fcur;
+ }
+ }
+ for (; tcur != to.first->zcolumn.end(); ++tcur)
+ {
+ rLog(rlCohomology, " %d", tcur->si->order);
+ nw.push_back(SNode(*tcur));
+ }
+ for (; fcur != from.first->zcolumn.end(); ++fcur)
+ {
+ rLog(rlCohomology, " %d", fcur->si->order);
+ nw.push_back(SNode(fcur->si, field_.mul(multiplier, fcur->coefficient), ci));
+ }
- std::set_symmetric_difference(to.cocycle.begin(), to.cocycle.end(), // this is catered to Z_2
- from.cocycle.begin(), from.cocycle.end(),
- std::back_insert_iterator<ZColumn>(nw),
- CompareSNode());
- for (typename ZColumn::iterator cur = to.cocycle.begin(); cur != to.cocycle.end(); ++cur)
+ for (typename ZColumn::iterator cur = to.first->zcolumn.begin(); cur != to.first->zcolumn.end(); ++cur)
cur->unlink();
- to.cocycle.swap(nw);
+ to.first->zcolumn.swap(nw);
- for (typename ZColumn::iterator cur = to.cocycle.begin(); cur != to.cocycle.end(); ++cur)
- {
+ for (typename ZColumn::iterator cur = to.first->zcolumn.begin(); cur != to.first->zcolumn.end(); ++cur)
cur->si->row.push_back(*cur);
- cur->cocycle = nw[0].cocycle;
- }
}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/field-arithmetic.h Thu Jul 09 00:59:32 2009 -0700
@@ -0,0 +1,64 @@
+#ifndef __FIELD_ARITHMETIC_H__
+#define __FIELD_ARITHMETIC_H__
+
+#include <vector>
+#include <gmpxx.h>
+
+class ZpField
+{
+ public:
+ typedef int Element;
+
+ ZpField(Element p = 2);
+
+ Element id() const { return 1; }
+ Element zero() const { return 0; }
+
+ Element neg(Element a) const { return p_ - a; }
+ Element add(Element a, Element b) const { return (a+b) % p_; }
+
+ Element inv(Element a) const { return inverses_[a]; }
+ Element mul(Element a, Element b) const { return (a*b) % p_; }
+ Element div(Element a, Element b) const { return mul(a, inv(b)); }
+
+ bool is_zero(Element a) const { return (a % p_) == 0; }
+
+ private:
+ Element p_;
+ std::vector<Element> inverses_;
+};
+
+ZpField::
+ZpField(Element p):
+ p_(p), inverses_(p_)
+{
+ for (Element i = 1; i < p_; ++i)
+ for (Element j = 1; j < p_; ++j)
+ if (mul(i,j) == 1)
+ {
+ inverses_[i] = j;
+ break;
+ }
+}
+
+class QField
+{
+ public:
+ typedef mpq_class Element;
+
+ QField() {}
+
+ Element id() const { return 1; }
+ Element zero() const { return 0; }
+
+ Element neg(Element a) const { return -a; }
+ Element add(Element a, Element b) const { return (a+b); }
+
+ Element inv(Element a) const { return id()/a; }
+ Element mul(Element a, Element b) const { return (a*b); }
+ Element div(Element a, Element b) const { return a/b; }
+
+ bool is_zero(Element a) const { return a == 0; }
+};
+
+#endif // __FIELD_ARITHMETIC_H__
--- a/include/utilities/indirect.h Thu May 14 14:04:43 2009 -0700
+++ b/include/utilities/indirect.h Thu Jul 09 00:59:32 2009 -0700
@@ -28,6 +28,27 @@
{ return IndirectComparison<Comparison>(cmp); }
+template<class Comparison_>
+class FirstComparison
+{
+ public:
+ typedef Comparison_ Comparison;
+
+ FirstComparison(Comparison cmp):
+ cmp_(cmp) {}
+
+ template<class Pair>
+ bool operator()(const Pair& a, const Pair& b) const
+ { return cmp_(a.first, b.first); }
+
+ private:
+ Comparison cmp_;
+};
+template<class Comparison>
+FirstComparison<Comparison> make_first_comparison(const Comparison& cmp)
+{ return FirstComparison<Comparison>(cmp); }
+
+
template<class Comparison>
struct ThreeOutcomeCompare: public Comparison
{