knotkit/algebra/module.h
Wojciech Politarczyk 203477c3a8 First shot at periodicity congruences
Verification of Przytycki's congruence works fine, however there are
some small modifications that could be made.

The, naive, approach to the verification of periodicity congruence for
the Khovanov polynomial didn't work. There are just too many cases to
check. Generation of all of them can exhaust the memory.
2016-11-24 15:28:31 +01:00

2152 lines
52 KiB
C++

#ifndef _KNOTKIT_ALGEBRA_MODULE_H
#define _KNOTKIT_ALGEBRA_MODULE_H
template<class R> class mod_map;
template<class R> class mod_span;
template<class R> class free_submodule;
template<class R> class quotient_module;
template<class R> class direct_sum;
template<class R> class tensor_product;
template<class R> class hom_module;
/* `module' is a bigraded module over a ring R. */
template<class R>
class module : public refcounted
{
private:
friend class reader;
friend class writer;
unsigned id;
static unsigned id_counter;
static map<unsigned, ptr<const module<R> > > reader_id_module;
static map<basedvector<unsigned, 1>,
ptr<const direct_sum<R> > > direct_sum_idx;
static map<basedvector<unsigned, 1>,
ptr<const tensor_product<R> > > tensor_product_idx;
static map<pair<unsigned, unsigned>,
ptr<const hom_module<R> > > hom_module_idx;
public:
module ()
{
id_counter ++;
id = id_counter;
}
module (const module &) = delete;
virtual ~module () { }
module &operator = (const module &) = delete;
public:
// the number of generators; n
virtual unsigned dim () const = 0;
// r; 1 <= r <= n
virtual unsigned free_rank () const = 0;
virtual grading generator_grading (unsigned i) const = 0;
virtual void show_generator (unsigned i) const = 0;
// r < i <= n
virtual R generator_ann (unsigned i) const = 0;
basedvector<grading, 1> grading_vector () const;
set<grading> gradings () const;
bool is_free () const { return dim () == free_rank (); }
bool is_zero (R c, unsigned i) const
{
if (i <= free_rank ())
return c == 0;
else
{
// ???
abort ();
// return generator_ann (i) | c;
}
}
R annihilator (R c, unsigned i) const
{
R iann = generator_ann (i);
return iann.div (iann.gcd (c));
}
bool isomorphic (ptr<const module<R> > m) const;
ptr<const quotient_module<R> > quotient (const mod_span<R> &span) const;
ptr<const quotient_module<R> > quotient (ptr<const free_submodule<R> > m) const;
ptr<const free_submodule<R> > submodule (const mod_span<R> &span) const;
multivariate_laurentpoly<Z> free_poincare_polynomial () const;
template<unsigned p>
multivariate_laurentpoly<Zp<p>> free_poincare_polynomial () const;
multivariate_laurentpoly<Z> free_delta_poincare_polynomial () const;
multivariate_laurentpoly<Z> free_ell_poincare_polynomial () const;
ptr<const direct_sum<R> >
add (basedvector<ptr<const module<R> >, 1> compound_summands) const;
ptr<const direct_sum<R> > add (ptr<const module<R> > m) const
{
basedvector<ptr<const module<R> >, 1> summands (2);
summands[1] = this;
summands[2] = m;
return add (summands);
}
// g -> (g, 0)
unsigned inject_1 (unsigned g, ptr<const module<R> > m) const
{
return g;
}
// g -> (0, g)
unsigned inject_2 (ptr<const module<R> > m, unsigned g) const
{
return dim () + g;
}
pair<unsigned, unsigned> project (ptr<const module<R> > m, unsigned g) const
{
if (g <= dim ())
return pair<unsigned, unsigned> (1, g);
else
return pair<unsigned, unsigned> (2, g - dim ());
}
ptr<const tensor_product<R> > tensor (ptr<const module<R> > m) const
{
basedvector<ptr<const module<R> >, 1> factors (2);
factors[1] = this;
factors[2] = m;
return tensor (factors);
}
ptr<const hom_module<R> > hom (ptr<const module<R> > to) const;
pair<unsigned, unsigned>
generator_factors (ptr<const module<R> > m, unsigned g) const
{
pair<unsigned, unsigned> p ((g - 1) % dim () + 1,
(g - 1) / dim () + 1);
assert (g == tensor_generators (p.first, m, p.second));
return p;
}
static ptr<const tensor_product<R> > tensor (basedvector<ptr<const module<R> >, 1> compound_factors);
unsigned tensor_generators (unsigned i, ptr<const module<R> > m, unsigned j) const
{
return (i - 1) + (j - 1) * dim () + 1;
}
virtual void append_direct_summands (basedvector<ptr<const module<R> >, 1> &summands) const
{
summands.append (this);
}
virtual void append_tensor_factors (basedvector<ptr<const module<R> >, 1> &factors) const
{
factors.append (this);
}
ptr<const free_submodule<R> > graded_piece (grading hq) const;
void write_self (writer &w) const { w.write_mod<R> (this); }
void show_self () const;
void display_self () const;
};
template<class R> unsigned module<R>::id_counter = 0;
template<class R> map<unsigned, ptr<const module<R> > > module<R>::reader_id_module;
template<class R> map<basedvector<unsigned, 1>,
ptr<const direct_sum<R> > > module<R>::direct_sum_idx;
template<class R> map<basedvector<unsigned, 1>,
ptr<const tensor_product<R> > > module<R>::tensor_product_idx;
template<class R> map<pair<unsigned, unsigned>,
ptr<const hom_module<R> > > module<R>::hom_module_idx;
template<class R>
class direct_sum : public module<R>
{
unsigned n;
basedvector<ptr<const module<R> >, 1> summands;
public:
direct_sum (basedvector<ptr<const module<R> >, 1> summands_)
: n(0),
summands(summands_)
{
#ifndef NDEBUG
for (unsigned i = 1; i <= summands.size (); i ++)
assert (summands[i]->is_free ());
#endif
for (unsigned i = 1; i <= summands.size (); i ++)
n += summands[i]->dim ();
}
~direct_sum () { }
unsigned dim () const { return n; }
unsigned free_rank () const { return n; }
grading generator_grading (unsigned i) const;
void show_generator (unsigned i) const;
R generator_ann (unsigned i) const { return R (0); }
void append_direct_summands (basedvector<ptr<const module<R> >, 1> &psummands) const
{
for (unsigned i = 1; i <= summands.size (); i ++)
psummands.append (summands[i]);
}
unsigned inject (unsigned i, unsigned g) const;
pair<unsigned, unsigned> project (unsigned g) const;
};
template<class R> grading
direct_sum<R>::generator_grading (unsigned i) const
{
pair<unsigned, unsigned> p = project (i);
return summands[p.first]->generator_grading (p.second);
}
template<class R> void
direct_sum<R>::show_generator (unsigned i) const
{
pair<unsigned, unsigned> p = project (i);
printf ("%d:", p.first);
return summands[p.first]->show_generator (p.second);
}
template<class R> unsigned
direct_sum<R>::inject (unsigned i, unsigned g) const
{
assert (i >= 1 && i <= summands.size ());
for (unsigned j = 1; j < i; j ++)
g += summands[j]->dim ();
return g;
}
template<class R> pair<unsigned, unsigned>
direct_sum<R>::project (unsigned g) const
{
assert (g <= n);
unsigned g0 = g;
for (unsigned j = 1; j <= summands.size (); j ++)
{
if (g <= summands[j]->dim ())
{
assert (inject (j, g) == g0);
return pair<unsigned, unsigned> (j, g);
}
else
g -= summands[j]->dim ();
}
abort (); // shouldn't get here
}
template<class R> ptr<const direct_sum<R> >
module<R>::add (basedvector<ptr<const module<R> >, 1> compound_summands) const
{
basedvector<ptr<const module<R> >, 1> summands;
for (unsigned i = 1; i <= compound_summands.size (); i ++)
compound_summands[i]->append_direct_summands (summands);
basedvector<unsigned, 1> summand_ids (summands.size ());
for (unsigned i = 1; i <= summands.size (); i ++)
summand_ids[i] = summands[i]->id;
pair<ptr<const direct_sum<R> > &, bool> p = direct_sum_idx.find (summand_ids);
if (!p.second)
p.first = new direct_sum<R> (summands);
return p.first;
}
template<class R>
class tensor_product : public module<R>
{
unsigned n;
basedvector<ptr<const module<R> >, 1> factors;
basedvector<unsigned, 1> generator_factors (unsigned g) const;
public:
tensor_product (basedvector<ptr<const module<R> >, 1> factors_)
: n(1),
factors(factors_)
{
#ifndef NDEBUG
for (unsigned i = 1; i <= factors.size (); i ++)
assert (factors[i]->is_free ());
#endif
for (unsigned i = 1; i <= factors.size (); i ++)
n *= factors[i]->dim ();
}
~tensor_product () { }
unsigned dim () const { return n; }
unsigned free_rank () const { return n; }
grading generator_grading (unsigned i) const;
void show_generator (unsigned i) const;
R generator_ann (unsigned i) const { return R (0); }
unsigned tensor_generators (basedvector<unsigned, 1> gs) const;
void append_tensor_factors (basedvector<ptr<const module<R> >, 1> &pfactors) const
{
for (unsigned i = 1; i <= factors.size (); i ++)
pfactors.append (factors[i]);
}
};
template<class R> ptr<const tensor_product<R> >
module<R>::tensor (basedvector<ptr<const module<R> >, 1> compound_factors)
{
basedvector<ptr<const module<R> >, 1> factors;
for (unsigned i = 1; i <= compound_factors.size (); i ++)
compound_factors[i]->append_tensor_factors (factors);
basedvector<unsigned, 1> factor_ids (factors.size ());
for (unsigned i = 1; i <= factors.size (); i ++)
factor_ids[i] = factors[i]->id;
pair<ptr<const tensor_product<R> > &, bool> p = tensor_product_idx.find (factor_ids);
if (!p.second)
p.first = new tensor_product<R> (factors);
return p.first;
}
template<class R> grading
tensor_product<R>::generator_grading (unsigned i) const
{
basedvector<unsigned, 1> gs = generator_factors (i);
assert (gs.size () == factors.size ());
grading gr;
for (unsigned i = 1; i <= factors.size (); i ++)
gr += factors[i]->generator_grading (gs[i]);
return gr;
}
template<class R> void
tensor_product<R>::show_generator (unsigned i) const
{
basedvector<unsigned, 1> gs = generator_factors (i);
assert (gs.size () == factors.size ());
printf ("o(");
for (unsigned i = 1; i <= factors.size (); i ++)
{
if (i > 1)
printf (",");
factors[i]->show_generator (gs[i]);
}
printf (")");
}
template<class R> unsigned
tensor_product<R>::tensor_generators (basedvector<unsigned, 1> gs) const
{
assert (gs.size () == factors.size ());
unsigned r = gs[gs.size ()] - 1;
for (unsigned i = gs.size () - 1; i >= 1; i --)
{
r *= factors[i]->dim ();
r += gs[i] - 1;
}
r ++;
return r;
}
template<class R> basedvector<unsigned, 1>
tensor_product<R>::generator_factors (unsigned g) const
{
basedvector<unsigned, 1> r (factors.size ());
unsigned g0 = g;
g --;
for (unsigned i = 1; i <= factors.size (); i ++)
{
r[i] = (g % factors[i]->dim ()) + 1;
g /= factors[i]->dim ();
}
assert (g == 0);
assert (tensor_generators (r) == g0);
return r;
}
template<class R>
class hom_module : public module<R>
{
public:
unsigned n;
ptr<const module<R> > from;
ptr<const module<R> > to;
public:
hom_module (ptr<const module<R> > from_,
ptr<const module<R> > to_)
: from(from_), to(to_)
{
assert (from->is_free ()
&& to->is_free ());
n = from->dim () * to->dim ();
}
~hom_module () { }
// e_ij -> ij
pair<unsigned, unsigned> generator_indices (unsigned g) const
{
unsigned d = from->dim ();
unsigned g0 = g;
g --;
pair<unsigned, unsigned> p ((g % d) + 1,
(g / d) + 1);
assert (generator (p.first, p.second) == g0);
return p;
}
// ij -> e_ij
unsigned generator (unsigned i, unsigned j) const
{
return (i - 1) + (j - 1) * from->dim () + 1;
}
unsigned dim () const { return n; }
unsigned free_rank () const { return n; }
grading generator_grading (unsigned i) const;
void show_generator (unsigned i) const;
R generator_ann (unsigned i) const { return R (0); }
linear_combination<R> map_as_element (const mod_map<R> &m) const;
};
template<class R> ptr<const hom_module<R> >
module<R>::hom (ptr<const module<R> > to) const
{
pair<ptr<const hom_module<R> > &, bool> p = hom_module_idx.find (pair<unsigned, unsigned>
(id, to->id));
if (!p.second)
p.first = new hom_module<R> (this, to);
return p.first;
}
template<class R> grading
hom_module<R>::generator_grading (unsigned i) const
{
pair<unsigned, unsigned> p = generator_indices (i);
return (to->generator_grading (p.second)
- from->generator_grading (p.first));
}
template<class R> void
hom_module<R>::show_generator (unsigned i) const
{
pair<unsigned, unsigned> p = generator_indices (i);
printf ("(");
from->show_generator (p.first);
printf (" -> ");
to->show_generator (p.second);
printf (")");
}
template<class R, class G>
class base_module : public module<R>
{
private:
G g;
public:
base_module () = delete;
base_module (const G &g_) : g(g_) { }
~base_module () { }
base_module &operator = (const base_module &) = delete;
unsigned dim () const { return g.dim (); }
unsigned free_rank () const { return g.free_rank (); }
grading generator_grading (unsigned i) const { return g.generator_grading (i); }
void show_generator (unsigned i) const { g.show_generator (i); }
R generator_ann (unsigned i) const { return g.generator_ann (i); }
};
template<class R>
class explicit_module : public module<R>
{
unsigned r;
basedvector<R, 1> ann;
basedvector<grading, 1> hq;
public:
explicit_module () = delete;
explicit_module (unsigned r_,
basedvector<R, 1> ann_,
basedvector<grading, 1> hq_)
: r(r_), ann(ann_), hq(hq_)
{
assert (hq.size () == r + ann.size ());
}
explicit explicit_module (unsigned r_, basedvector<grading, 1> hq_) : r(r_), hq(hq_) { }
~explicit_module () { }
explicit_module &operator = (const explicit_module &) = delete;
unsigned dim () const { return r + ann.size (); }
unsigned free_rank () const { return r; }
grading generator_grading (unsigned i) const { return hq[i]; }
void show_generator (unsigned i) const { printf ("%d", i); }
R generator_ann (unsigned i) const { return ann[i - r]; }
};
template<class R>
class free_submodule : public module<R>
{
friend class module<R>;
ptr<const module<R> > parent;
basedvector<linear_combination<R>, 1> gens;
basedvector<unsigned, 1> pivots;
public:
free_submodule (ptr<const module<R> > parent_,
basedvector<linear_combination<R>, 1> gens_,
basedvector<unsigned, 1> pivots_)
: parent(parent_),
gens(gens_),
pivots(pivots_)
{ }
~free_submodule () { }
free_submodule &operator = (const free_submodule &) = delete;
ptr<const module<R> > parent_module () const { return parent; }
unsigned dim () const { return gens.size (); }
unsigned free_rank () const { return gens.size (); }
grading generator_grading (unsigned i) const { return gens[i].hq (); }
void show_generator (unsigned i) const { show (gens[i]); }
R generator_ann (unsigned i) const { abort (); }
linear_combination<R> inject_generator (unsigned i) const { return gens[i]; }
linear_combination<R> inject (linear_combination<R> v) const;
mod_map<R> injection_map () const;
linear_combination<R> restrict (linear_combination<R> v0) const;
ptr<const free_submodule<R> > restrict_submodule (ptr<const free_submodule<R> > m) const;
ptr<const free_submodule<R> > intersection (ptr<const free_submodule<R> > m) const;
ptr<const free_submodule<R> > plus (ptr<const free_submodule<R> > m) const;
};
template<class R>
class quotient_module : public module<R>
{
friend class module<R>;
ptr<const module<R> > parent;
// note: these get filled in by module::quotient
basedvector<R, 1> ann;
// parent lc representing i
basedvector<linear_combination<R>, 1> rep;
// map from parent generator to lc in quotient
basedvector<map<unsigned, R>, 1> pi;
public:
quotient_module (const quotient_module &) = delete;
quotient_module (ptr<const module<R> > parent_)
: parent(parent_)
{ }
~quotient_module () { }
quotient_module &operator = (const quotient_module &) = delete;
ptr<const module<R> > parent_module () const { return parent; }
unsigned dim () const { return rep.size (); }
unsigned free_rank () const { return rep.size () - ann.size (); }
grading generator_grading (unsigned i) const { return rep[i].hq (); }
void show_generator (unsigned i) const
{
show (rep[i]);
printf ("/~");
}
R generator_ann (unsigned i) const
{
unsigned r = free_rank ();
assert (i > r);
return ann[i - r];
}
linear_combination<R> project_generator (unsigned i) const
{
assert (i >= 1 && i <= parent->dim ());
linear_combination<R> r (this);
for (typename map<unsigned, R>::const_iter j = pi[i]; j; j ++)
r.muladd (j.val (), j.key ());
return r;
}
linear_combination<R> project (linear_combination<R> v) const
{
assert (v.m == parent);
linear_combination<R> r (this);
for (linear_combination_const_iter<R> i = v; i; i ++)
r.muladd (i.val (), project_generator (i.key ()));
return r;
}
linear_combination<R> generator_rep (unsigned i) const { return rep[i]; }
};
template<class R>
class map_impl : public refcounted
{
public:
ptr<const module<R> > from;
ptr<const module<R> > to;
public:
map_impl (const map_impl &) = delete;
map_impl (ptr<const module<R> > fromto)
: from(fromto), to(fromto)
{ }
map_impl (ptr<const module<R> > from_, ptr<const module<R> > to_)
: from(from_), to(to_)
{ }
virtual ~map_impl () { }
map_impl &operator = (const map_impl &) = delete;
virtual const linear_combination<R> column (unsigned i) const = 0;
virtual const linear_combination<R> column_copy (unsigned i) const { return column (i); }
linear_combination<R> map (const linear_combination<R> &lc) const
{
linear_combination<R> r (this->to);
for (linear_combination_const_iter<R> i = lc; i; i ++)
r.muladd (i.val (), column (i.key ()));
return r;
};
};
template<class R>
class explicit_map_impl : public map_impl<R>
{
basedvector<linear_combination<R>, 1> columns;
public:
explicit_map_impl (ptr<const module<R> > fromto,
basedvector<linear_combination<R>, 1> columns_)
: map_impl<R>(fromto),
columns(columns_)
{ }
explicit_map_impl (ptr<const module<R> > from, ptr<const module<R> > to,
basedvector<linear_combination<R>, 1> columns_)
: map_impl<R>(from, to),
columns(columns_)
{ }
~explicit_map_impl () { }
const linear_combination<R> column (unsigned i) const { return columns[i]; }
const linear_combination<R> column_copy (unsigned i) const
{
return linear_combination<R> (COPY, columns[i]);
}
};
template<class R>
class zero_map_impl : public map_impl<R>
{
public:
zero_map_impl (ptr<const module<R> > fromto) : map_impl<R>(fromto) { }
zero_map_impl (ptr<const module<R> > from, ptr<const module<R> > to) : map_impl<R>(from, to) { }
const linear_combination<R> column (unsigned i) const { return linear_combination<R> (this->to); }
};
template<class R>
class id_map_impl : public map_impl<R>
{
public:
id_map_impl (ptr<const module<R> > fromto) : map_impl<R>(fromto) { }
id_map_impl (ptr<const module<R> > from, ptr<const module<R> > to) : map_impl<R>(from, to) { }
const linear_combination<R> column (unsigned i) const
{
linear_combination<R> r (this->to);
r.muladd (1, i);
return r;
}
};
template<class R>
class composition_impl : public map_impl<R>
{
// f(g(x))
ptr<const map_impl<R> > f, g;
public:
composition_impl (ptr<const map_impl<R> > f_, ptr<const map_impl<R> > g_)
: map_impl<R>(g_->from, f_->to),
f(f_),
g(g_)
{
assert (g->to == f->from);
}
const linear_combination<R> column (unsigned i) const
{
return f->map (g->column (i));
}
};
template<class R>
class direct_sum_impl : public map_impl<R>
{
// f\oplus g
ptr<const map_impl<R> > f, g;
public:
direct_sum_impl (ptr<const map_impl<R> > f_, ptr<const map_impl<R> > g_)
: map_impl<R>(f_->from->add (g_->from),
f_->to->add (g_->to)),
f(f_),
g(g_)
{
}
const linear_combination<R> column (unsigned i) const
{
pair<unsigned, unsigned> p = f->from->project (g->from, i);
linear_combination<R> r (this->to);
if (p.first == 1)
{
for (linear_combination_const_iter<R> j = f->column (p.second); j; j ++)
r.muladd (j.val (), f->to->inject_1 (j.key (), g->to));
}
else
{
assert (p.first == 2);
for (linear_combination_const_iter<R> j = g->column (p.second); j; j ++)
r.muladd (j.val (), f->to->inject_2 (g->to, j.key ()));
}
return r;
}
};
template<class R>
class tensor_impl : public map_impl<R>
{
// f\otimes g
ptr<const map_impl<R> > f, g;
public:
tensor_impl (ptr<const map_impl<R> > f_, ptr<const map_impl<R> > g_)
: map_impl<R>(f_->from->tensor (g_->from),
f_->to->tensor (g_->to)),
f(f_),
g(g_)
{
}
const linear_combination<R> column (unsigned i) const
{
pair<unsigned, unsigned> p = f->from->generator_factors (g->from, i);
// ??
return f->column (p.first).tensor (g->column (p.second));
}
};
template<class R>
class map_builder
{
public:
ptr<const module<R> > from, to;
basedvector<linear_combination<R>, 1> columns;
void init ();
public:
map_builder (ptr<const module<R> > fromto) : from(fromto), to(fromto) { init (); }
map_builder (ptr<const module<R> > fromto, int i) : from(fromto), to(fromto)
{
init ();
if (i == 1)
{
for (unsigned i = 1; i <= from->dim (); i ++)
columns[i].muladd (1, i);
}
else
assert (i == 0);
}
map_builder (ptr<const module<R> > from_, ptr<const module<R> > to_)
: from(from_), to(to_)
{
init ();
}
linear_combination<R> &operator [] (unsigned i) { return columns[i]; }
const linear_combination<R> &operator [] (unsigned i) const { return columns[i]; }
};
template<class R> void
map_builder<R>::init ()
{
columns.resize (from->dim ());
for (unsigned i = 1; i <= from->dim (); i ++)
columns[i] = linear_combination<R> (to);
}
template<class R>
class mod_map
{
// ???
enum impl_ctor { IMPL };
ptr<const map_impl<R> > impl;
mod_map (impl_ctor, ptr<const map_impl<R> > impl_) : impl(impl_) { }
public:
mod_map () { }
mod_map (const mod_map &m) : impl(m.impl) { }
mod_map (ptr<const module<R> > fromto) : impl(new zero_map_impl<R>(fromto)) { }
mod_map (ptr<const module<R> > fromto, int i)
{
if (i == 1)
impl = new id_map_impl<R> (fromto);
else
{
assert (i == 0);
impl = new zero_map_impl<R> (fromto);
}
}
mod_map (ptr<const module<R> > from, ptr<const module<R> > to)
: impl(new zero_map_impl<R> (from, to))
{ }
mod_map (ptr<const module<R> > fromto,
basedvector<linear_combination<R>, 1> columns)
: impl(new explicit_map_impl<R> (fromto, columns))
{ }
mod_map (ptr<const module<R> > from,
ptr<const module<R> > to,
basedvector<linear_combination<R>, 1> columns)
: impl(new explicit_map_impl<R> (from, to, columns))
{ }
mod_map (const map_builder<R> &b)
: impl(new explicit_map_impl<R> (b.from, b.to, b.columns))
{ }
mod_map (reader &r)
{
ptr<const module<R> > from = r.read_mod<R> ();
ptr<const module<R> > to = r.read_mod<R> ();
basedvector<linear_combination<R>, 1> columns (r);
impl = new explicit_map_impl<R> (from, to, columns);
}
~mod_map () { }
mod_map &operator = (const mod_map &m) { impl = m.impl; return *this; }
ptr<const module<R> > domain () const { return impl->from; }
ptr<const module<R> > codomain () const { return impl->to; }
bool operator == (const mod_map &m) const
{
assert (impl->from == m.impl->from);
assert (impl->to == m.impl->to);
for (unsigned i = 1; i <= impl->from->dim (); i ++)
{
if (impl->column (i) != m.impl->column (i))
return 0;
}
return 1;
}
bool operator != (const mod_map &m) const { return !operator == (m); }
bool operator == (int x) const
{
R c (x);
assert (c == 0);
for (unsigned i = 1; i <= impl->from->dim (); i ++)
{
if (impl->column (i) != 0)
return 0;
}
return 1;
}
bool operator != (int x) const { return !operator == (x); }
const linear_combination<R> column (unsigned i) const { return impl->column (i); }
const linear_combination<R> operator [] (unsigned i) const { return impl->column (i); }
const linear_combination<R> column_copy (unsigned i) const { return impl->column_copy (i); }
linear_combination<R> map (const linear_combination<R> &lc) const { return impl->map (lc); }
mod_map compose (const mod_map &m) const
{
return mod_map (IMPL,
new composition_impl<R> (impl, m.impl));
}
// ??? in the sense of direct sum
mod_map add (const mod_map &m) const
{
return mod_map (IMPL,
new direct_sum_impl<R> (impl, m.impl));
}
mod_map tensor (const mod_map &m) const
{
return mod_map (IMPL,
new tensor_impl<R> (impl, m.impl));
}
// ?? add and other map operations should not be explicit
mod_map operator + (const mod_map &m) const;
mod_map operator * (const R &c) const;
bool homogeneous () const;
void check_grading (grading delta) const;
// inj : ker -> from
ptr<const free_submodule<R> > kernel () const;
// inj : im -> to
ptr<const free_submodule<R> > image () const;
ptr<const free_submodule<R> > image (basedvector<linear_combination<R>, 1> vs) const;
ptr<const quotient_module<R> > cokernel () const;
ptr<const quotient_module<R> > homology () const;
mod_map restrict_from (ptr<const free_submodule<R> > new_from) const;
mod_map restrict_to (ptr<const free_submodule<R> > new_to) const;
mod_map restrict (ptr<const free_submodule<R> > new_from,
ptr<const free_submodule<R> > new_to) const;
mod_map restrict (ptr<const free_submodule<R> > new_fromto) const
{
return restrict (new_fromto, new_fromto);
}
mod_map induced_map_to (ptr<const quotient_module<R> > new_to);
mod_map induced_map (ptr<const quotient_module<R> > new_fromto);
mod_map graded_piece (grading hq) const;
// ???
basedvector<linear_combination<R>, 1> explicit_columns () const;
void write_self (writer &w) const
{
// write explicitly
write (w, *impl->from);
write (w, *impl->to);
write (w, explicit_columns ());
}
void show_self () const;
void display_self () const;
};
template<class R> linear_combination<R>
hom_module<R>::map_as_element (const mod_map<R> &m) const
{
assert (from == m.domain ()
&& to == m.codomain ());
linear_combination<R> r (this);
for (unsigned i = 1; i <= from->dim (); i ++)
{
for (linear_combination_const_iter<R> j = m.column (i); j; j ++)
r.muladd (j.val (), generator (i, j.key ()));
}
return r;
}
template<class R>
class mod_span
{
public:
basedvector<linear_combination<R>, 1> gens;
basedvector<unsigned, 1> pivots;
public:
mod_span (ptr<const module<R> > mod, basedvector<linear_combination<R>, 1> xs);
~mod_span () { }
};
template<class R>
class quotient_helper
{
public:
ptr<const module<R> > mod;
// rows of the presentation matrix
basedvector<linear_combination<R>, 1> rows;
basedvector<linear_combination<R>, 1> generators;
basedvector<linear_combination<R>, 1> generators_inv;
bool improve_pivot_row (unsigned i, unsigned j, unsigned i2);
bool improve_pivot_column (unsigned i, unsigned j, unsigned j2);
void improve_pivot (unsigned i, unsigned j);
public:
quotient_helper (ptr<const module<R> > mod_, basedvector<linear_combination<R>, 1> rows);
void normalize ();
};
template<class R> bool
module<R>::isomorphic (ptr<const module<R> > m) const
{
if (dim () != m->dim ()
|| free_rank () != m->free_rank ())
return 0;
set<grading> gradings;
for (unsigned i = 1; i <= dim (); i ++)
gradings += generator_grading (i);
for (set_const_iter<grading> i = gradings; i; i ++)
{
grading hq = i.val ();
unsigned hq_r = 0,
m_hq_r = 0;
basedvector<R, 1> hq_ann, m_hq_ann;
for (unsigned i = 1; i <= dim (); i ++)
{
if (generator_grading (i) == hq)
{
if (i <= free_rank ())
hq_r ++;
else
{
hq_ann.append (generator_ann (i));
}
}
if (m->generator_grading (i) == hq)
{
if (i <= free_rank ())
m_hq_r ++;
else
{
m_hq_ann.append (m->generator_ann (i));
}
}
}
if (hq_r != m_hq_r)
return 0;
if (hq_ann.size () != m_hq_ann.size ())
return 0;
for (unsigned i = 1; i <= hq_ann.size (); i ++)
{
R a = hq_ann[i],
m_a = m_hq_ann[i];
if (! (a | m_a && m_a | a))
return 0;
}
}
return 1;
}
template<class R>
quotient_helper<R>::quotient_helper (ptr<const module<R> > mod_,
basedvector<linear_combination<R>, 1> rows_)
: mod(mod_),
rows(rows_),
generators(mod->dim ()),
generators_inv(mod->dim ())
{
assert (mod->dim () == mod->free_rank ());
for (unsigned i = 1; i <= mod->dim (); i ++)
{
linear_combination<R> v (mod);
v.muladd (1, i);
generators[i] = v;
linear_combination<R> vinv (mod);
vinv.muladd (1, i);
generators_inv[i] = vinv;
}
}
template<class R> bool
quotient_helper<R>::improve_pivot_row (unsigned i, unsigned j, unsigned i2)
{
assert (i != i2);
const linear_combination<R> &r = rows[i],
&r2 = rows[i2];
R rc = r(j),
r2c = r2(j);
if (r2c == 0)
return 0;
#if 0
if (rc | r2c)
{
R q = r2c.div (rc);
rows[i2].mulsub (q, r);
assert (rows[i2](j) == 0);
return 0;
}
#endif
tuple<R, R, R> t = rc.extended_gcd (r2c);
assert (get<0> (t) == rc*get<1> (t) + get<2> (t)*r2c);
rows[i] = r*get<1> (t) + r2*get<2> (t);
rows[i2] = (rc.div (get<0> (t)))*r2 - (r2c.div (get<0> (t)))*r;
assert ((rc | r2c) == rc.divides (get<0> (t)));
assert (!rc.divides (get<0> (t)) || rows[i2](j) == 0);
return !rc.divides (get<0> (t));
}
template<class R> bool
quotient_helper<R>::improve_pivot_column (unsigned i, unsigned j, unsigned j2)
{
assert (j != j2);
#if 0
basedvector<linear_combination<R>, 1> orig_row_image (rows.size ());
for (unsigned k = 1; k <= rows.size (); k ++)
{
linear_combination<R> r (mod);
for (linear_combination_const_iter<R> l = rows[k]; l; l ++)
r.muladd (l.val (), generators[l.key ()]);
orig_row_image[k] = r;
}
#endif
const linear_combination<R> &r = rows[i];
R rc = r(j),
rc2 = r(j2);
assert (rc != 0);
if (rc2 == 0)
return 0;
assert (generators[j].hq () == generators[j2].hq ());
#if 0
if (rc | rc2)
{
R q = rc2.div (rc);
for (unsigned k = 1; k <= rows.size (); k ++)
{
linear_combination<R> &rk = rows[k];
rk.mulsub (rk(j) * q, j2);
}
assert (r(j2) == 0);
return 0;
}
#endif
tuple<R, R, R> t = rc.extended_gcd (rc2);
assert (get<0> (t) == rc*get<1> (t) + get<2> (t)*rc2);
for (unsigned k = 1; k <= rows.size (); k ++)
{
linear_combination<R> &rk = rows[k];
R rkc = rk(j),
rkc2 = rk(j2);
rk.set_coeff (rkc*get<1> (t) + rkc2*get<2> (t),
j);
rk.set_coeff (rkc2*(rc.div (get<0> (t))) - rkc*(rc2.div (get<0> (t))),
j2);
}
linear_combination<R> g = generators[j],
g2 = generators[j2];
assert (g.hq () == g2.hq ());
generators[j] = (rc.div (get<0> (t))) * g + (rc2.div (get<0> (t))) * g2;
generators[j2] = get<1> (t) * g2 - get<2> (t) * g;
#if 0
for (unsigned k = 1; k <= rows.size (); k ++)
{
linear_combination<R> r2 (mod);
for (linear_combination_const_iter<R> l = rows[k]; l; l ++)
r2.muladd (l.val (), generators[l.key ()]);
assert (r2 == orig_row_image[k]);
}
#endif
for (unsigned k = 1; k <= mod->dim (); k ++)
{
linear_combination<R> &ginv = generators_inv[k];
R d = ginv(j),
d2 = ginv(j2);
ginv.set_coeff (get<1> (t)*d + get<2> (t)*d2, j);
ginv.set_coeff (rc.div (get<0> (t)) * d2 - rc2.div (get<0> (t)) * d, j2);
}
#if 0
for (unsigned k = 1; k <= mod->dim (); k ++)
{
linear_combination<R> r (mod);
r.muladd (1, k);
linear_combination<R> r2 (mod);
for (linear_combination_const_iter<R> l = generators_inv[k]; l; l ++)
r2.muladd (l.val (), generators[l.key ()]);
assert (r == r2);
}
#endif
assert ((rc | rc2) == rc.divides (get<0> (t)));
assert (!rc.divides (get<0> (t)) || r(j2) == 0);
return !rc.divides (get<0> (t));
}
template<class R> void
quotient_helper<R>::improve_pivot (unsigned i, unsigned j)
{
for (;;)
{
bool changed = 0;
for (unsigned k = 1; k <= rows.size (); k ++)
{
if (k == i)
continue;
if (improve_pivot_row (i, j, k))
changed = 1;
}
for (unsigned k = 1; k <= mod->dim (); k ++)
{
if (k == j)
continue;
if (improve_pivot_column (i, j, k))
changed = 1;
}
#if 0
L:
for (linear_combination_const_iter<R> k = rows[i]; k; k ++)
{
if (k.key () != j)
{
if (improve_pivot_column (i, j, k.key ()))
changed = 1;
goto L;
}
}
#endif
if (!changed)
return;
}
}
template<class R> void
quotient_helper<R>::normalize ()
{
for (unsigned i = 1; i <= rows.size (); i ++)
{
if (rows[i] == 0)
continue;
pair<unsigned, R> p = rows[i].head ();
improve_pivot (i, p.first);
}
for (unsigned i = 1; i <= rows.size (); i ++)
{
if (rows[i] == 0)
continue;
assert (rows[i].card () == 1);
pair<unsigned, R> p = rows[i].head ();
grading ihq = generators[p.first].hq ();
for (unsigned j = i; j <= rows.size (); j ++)
{
if (rows[j] == 0)
continue;
assert (rows[j].card () == 1);
pair<unsigned, R> q = rows[j].head ();
grading jhq = generators[q.first].hq ();
if (ihq != jhq)
continue;
if (p.second | q.second)
continue;
rows[i] += rows[j];
improve_pivot (i, p.first);
#ifndef NDEBUG
assert (rows[i].card () == 1);
pair<unsigned, R> p2 = rows[i].head ();
assert (p2.first == p.first);
assert (rows[j].card () == 1);
pair<unsigned, R> q2 = rows[j].head ();
assert (q2.first == q.first);
assert (p2.second | q2.second);
#endif
}
}
}
template<class R> ptr<const quotient_module<R> >
module<R>::quotient (const mod_span<R> &span) const
{
unsigned n = dim ();
quotient_helper<R> h (this, span.gens);
h.normalize ();
basedvector<R, 1> ann (n);
for (unsigned i = 1; i <= h.rows.size (); i ++)
{
if (h.rows[i] == 0)
continue;
assert (h.rows[i].card () == 1);
pair<unsigned, R> p = h.rows[i].head ();
ann[p.first] = p.second;
}
ptr<quotient_module<R> > Q = new quotient_module<R> (this);
unsigned quot_r = 0;
basedvector<linear_combination<R>, 1> quot_rep;
basedvector<unsigned, 1> generator_quot_gen (n);
for (unsigned i = 1; i <= n; i ++)
generator_quot_gen[i] = 0;
for (unsigned i = 1; i <= n; i ++)
{
if (ann[i] == 0)
{
quot_r ++;
quot_rep.append (h.generators[i]);
generator_quot_gen[i] = quot_r;
}
}
unsigned quot_n = quot_r;
basedvector<R, 1> quot_ann;
for (unsigned i = 1; i <= n; i ++)
{
R a = ann[i];
if (a != 0 && !a.is_unit ())
{
quot_n ++;
quot_rep.append (h.generators[i]);
assert (generator_quot_gen[i] == 0);
generator_quot_gen[i] = quot_n;
quot_ann.append (a);
}
}
basedvector<map<unsigned, R>, 1> quot_pi (n);
for (unsigned i = 1; i <= n; i ++)
{
map<unsigned, R> v;
for (linear_combination_const_iter<R> j = h.generators_inv[i]; j; j ++)
{
unsigned qg = generator_quot_gen[j.key ()];
if (qg)
v.push (qg, j.val ());
}
quot_pi[i] = v;
}
assert (quot_n >= quot_r);
assert (quot_rep.size () == quot_n);
assert (quot_ann.size () == quot_n - quot_r);
Q->rep = quot_rep;
Q->pi = quot_pi;
Q->ann = quot_ann;
assert (Q->dim () == quot_n);
assert (Q->free_rank () == quot_r);
return Q;
}
template<class R> ptr<const quotient_module<R> >
module<R>::quotient (ptr<const free_submodule<R> > m) const
{
assert (m->parent_module () == this);
// ??? duplicate
mod_span<R> span (this, m->gens);
return quotient (span);
}
template<class R> ptr<const free_submodule<R> >
module<R>::submodule (const mod_span<R> &span) const
{
assert (free_rank () == dim ());
return new free_submodule<R> (this,
span.gens,
span.pivots);
}
template<class R> multivariate_laurentpoly<Z>
module<R>::free_poincare_polynomial () const
{
multivariate_laurentpoly<Z> r;
for (unsigned i = 1; i <= free_rank (); i ++)
{
grading hq = generator_grading (i);
multivariate_laurent_monomial m;
m.push_exponent (1, hq.h);
m.push_exponent (2, hq.q);
r.muladdeq (1, m);
}
return r;
}
template<class R> template <unsigned p> multivariate_laurentpoly<Zp<p>>
module<R>::free_poincare_polynomial() const {
multivariate_laurentpoly<Zp<p>> r;
for (unsigned i = 1; i <= free_rank (); i ++)
{
grading hq = generator_grading (i);
multivariate_laurent_monomial m;
m.push_exponent (1, hq.h);
m.push_exponent (2, hq.q);
r.muladdeq (1, m);
}
return r;
}
template<class R> multivariate_laurentpoly<Z>
module<R>::free_delta_poincare_polynomial () const
{
multivariate_laurentpoly<Z> r;
for (unsigned i = 1; i <= free_rank (); i ++)
{
grading hq = generator_grading (i);
multivariate_laurent_monomial m;
m.push_exponent (1, hq.q - 2*hq.h);
r.muladdeq (1, m);
}
return r;
}
template<class R> multivariate_laurentpoly<Z>
module<R>::free_ell_poincare_polynomial () const
{
multivariate_laurentpoly<Z> r;
for (unsigned i = 1; i <= free_rank (); i ++)
{
grading hq = generator_grading (i);
multivariate_laurent_monomial m;
m.push_exponent (1, hq.h - hq.q);
r.muladdeq (1, m);
}
return r;
}
template<class R> basedvector<grading, 1>
module<R>::grading_vector () const
{
basedvector<grading, 1> v (dim ());
for (unsigned i = 1; i <= dim (); i ++)
v[i] = generator_grading (i);
return v;
}
template<class R> set<grading>
module<R>::gradings () const
{
set<grading> gs;
for (unsigned i = 1; i <= dim (); i ++)
gs += generator_grading (i);
return gs;
}
template<class R> ptr<const free_submodule<R> >
module<R>::graded_piece (grading hq) const
{
basedvector<linear_combination<R>, 1> s;
for (unsigned i = 1; i <= dim (); i ++)
{
grading ihq = generator_grading (i);
if (ihq.h == hq.h
&& ihq.q == hq.q)
{
linear_combination<R> v (this);
v.muladd (1, i);
s.append (v);
}
}
mod_span<R> span (this, s);
return submodule (span);
}
template<class R> void
module<R>::show_self () const
{
printf ("module/");
R::show_ring ();
printf (" %p %d", this, dim ());
}
template<class R> void
module<R>::display_self () const
{
show_self (); newline ();
printf (" free_rank %d\n", free_rank ());
for (unsigned i = 1; i <= dim (); i ++)
{
printf (" %d ", i);
show (generator_grading (i));
printf (" ");
show_generator (i);
if (i > free_rank ())
{
printf (": ");
R iann = generator_ann (i);
if (iann.is_unit ())
printf ("0\n");
else
{
R::show_ring ();
if (iann != 0)
{
printf ("/");
show (iann);
R::show_ring ();
}
}
}
newline ();
}
}
template<class R> linear_combination<R>
free_submodule<R>::inject (linear_combination<R> v) const
{
linear_combination<R> r (parent);
for (linear_combination_const_iter<R> i = v; i; i ++)
r.muladd (i.val (), gens[i.key ()]);
return r;
}
template<class R> mod_map<R>
free_submodule<R>::injection_map () const
{
return mod_map<R> (this, parent, gens);
}
template<class R> mod_map<R>
mod_map<R>::induced_map_to (ptr<const quotient_module<R> > new_to)
{
assert (new_to->parent_module () == impl->to);
basedvector<linear_combination<R>, 1> v (impl->from->dim ());
for (unsigned i = 1; i <= impl->from->dim (); i ++)
v[i] = new_to->project (column (i));
return new explicit_map_impl<R> (impl->from, new_to, v);
}
template<class R> mod_map<R>
mod_map<R>::induced_map (ptr<const quotient_module<R> > new_fromto)
{
assert (impl->from == new_fromto->parent_module ());
assert (impl->to == new_fromto->parent_module ());
// ??? doesn't check induced map is well-defined
basedvector<linear_combination<R>, 1> v (new_fromto->dim ());
for (unsigned i = 1; i <= new_fromto->dim (); i ++)
v[i] = new_fromto->project (map (new_fromto->generator_rep (i)));
return new explicit_map_impl<R> (new_fromto, v);
}
template<class R> mod_map<R>
mod_map<R>::graded_piece (grading hq) const
{
basedvector<linear_combination<R>, 1> v (impl->from->dim ());
for (unsigned i = 1; i <= impl->from->dim (); i ++)
{
grading ihq = impl->from->generator_grading (i);
linear_combination<R> c = column (i);
linear_combination<R> d (impl->to);
for (linear_combination_const_iter<R> j = c; j; j ++)
{
grading jhq = impl->from->generator_grading (j.key ());
if (jhq.h - ihq.h == hq.h
&& jhq.q - ihq.q == hq.q)
d.muladd (j.val (), j.key ());
}
v[i] = d;
}
return mod_map (IMPL, new explicit_map_impl<R> (impl->from, impl->to, v));
}
template<class R> mod_map<R>
mod_map<R>::restrict_from (ptr<const free_submodule<R> > new_from) const
{
assert (new_from->parent_module () == impl->from);
basedvector<linear_combination<R>, 1> v (new_from->dim ());
for (unsigned i = 1; i <= new_from->dim (); i ++)
v[i] = map (new_from->inject_generator (i));
return mod_map (IMPL, new explicit_map_impl<R> (new_from, impl->to, v));
}
template<class R> mod_map<R>
mod_map<R>::restrict_to (ptr<const free_submodule<R> > new_to) const
{
assert (new_to->parent_module () == impl->to);
basedvector<linear_combination<R>, 1> v (impl->from->dim ());
for (unsigned i = 1; i <= impl->from->dim (); i ++)
v[i] = new_to->restrict (column (i));
return new explicit_map_impl<R> (impl->from, new_to, v);
}
template<class R> mod_map<R>
mod_map<R>::restrict (ptr<const free_submodule<R> > new_from,
ptr<const free_submodule<R> > new_to) const
{
assert (new_from->parent_module () == impl->from);
assert (new_to->parent_module () == impl->to);
basedvector<linear_combination<R>, 1> v (new_from->dim ());
for (unsigned i = 1; i <= new_from->dim (); i ++)
v[i] = new_to->restrict (map (new_from->inject_generator (i)));
return mod_map (IMPL, new explicit_map_impl<R> (new_from, new_to, v));
}
template<class R> linear_combination<R>
free_submodule<R>::restrict (linear_combination<R> v0) const
{
assert (v0.m == parent);
linear_combination<R> v (COPY, v0);
linear_combination<R> r (this);
for (unsigned i = 1; i <= gens.size (); i ++)
{
unsigned j = pivots[i];
R vc = v(j);
if (vc != 0)
{
const linear_combination<R> &g = gens[i];
R gc = g(j);
assert (gc | vc);
R q = vc.div (gc);
v.mulsub (q, g);
r.muladd (q, i);
}
}
assert (v == 0);
assert (inject (r) == v0);
return r;
}
template<class R> ptr<const free_submodule<R> >
free_submodule<R>::restrict_submodule (ptr<const free_submodule<R> > m) const
{
assert (m->parent == parent);
basedvector<linear_combination<R>, 1> span (m->dim ());
for (unsigned i = 1; i <= m->dim (); i ++)
span[i] = restrict (m->inject_generator (i));
mod_span<R> span2 (this, span);
return this->submodule (span2);
}
template<class R> ptr<const free_submodule<R> >
free_submodule<R>::intersection (ptr<const free_submodule<R> > m) const
{
assert (parent == m->parent);
unsigned md = m->dim (),
d = dim ();
basedvector<linear_combination<R>, 1> intr;
basedvector<linear_combination<R>, 1> hperp,
hproj;
basedvector<unsigned, 1> hpivots;
for (unsigned i = 1; i <= md; i ++)
{
linear_combination<R> perp (COPY, m->gens[i]),
proj (parent);
for (unsigned j = 1; j <= d; j ++)
{
unsigned k = pivots[j];
if (perp % k)
{
const linear_combination<R> &g = gens[j];
R c = g(k);
R d = perp(k);
assert (c | d);
R q = d.div (c);
perp.mulsub (q, g);
proj.mulsub (q, g);
assert (! (perp % k));
}
}
for (unsigned j = 1; j <= hpivots.size (); j ++)
{
unsigned k = hpivots[j];
if (perp % k)
{
const linear_combination<R> &h = hperp[j];
R c = h(k);
R d = perp(k);
assert (c | d);
R q = d.div (c);
perp.mulsub (q, h);
proj.mulsub (q, hproj[j]);
assert (! (perp % k));
}
}
if (perp == 0)
intr.append (proj);
else
{
hperp.append (perp);
hproj.append (proj);
hpivots.append (perp.head ().first);
}
}
mod_span<R> span (parent, intr);
return parent->submodule (span);
}
template<class R> ptr<const free_submodule<R> >
free_submodule<R>::plus (ptr<const free_submodule<R> > m) const
{
assert (parent == m->parent);
basedvector<linear_combination<R>, 1> s;
for (unsigned i = 1; i <= dim (); i ++)
s.append (gens[i]);
for (unsigned i = 1; i <= m->dim (); i ++)
s.append (m->gens[i]);
mod_span<R> span (parent, s);
return parent->submodule (span);
}
template<class R> bool
mod_map<R>::homogeneous () const
{
for (unsigned i = 1; i <= impl->from->dim (); i ++)
{
if (column (i) != 0)
{
grading dhq = column (i).hq () - impl->from->generator_grading (i);
for (unsigned j = i + 1; j <= impl->from->dim (); j ++)
{
if (column (j) != 0
&& dhq != column (j).hq () - impl->from->generator_grading (j))
return 0;
}
return 1;
}
}
return 1;
}
template<class R> void
mod_map<R>::check_grading (grading delta) const
{
for (unsigned i = 1; i <= impl->from->dim (); i ++)
{
if (column (i) != 0)
assert (column (i).hq () - impl->from->generator_grading (i) == delta);
}
}
template<class R> mod_map<R>
mod_map<R>::operator * (const R &c) const
{
basedvector<linear_combination<R>, 1> v (impl->from->dim ());
for (unsigned i = 1; i <= impl->from->dim (); i ++)
v[i] = c*column (i);
return mod_map (IMPL, new explicit_map_impl<R> (impl->from, impl->to, v));
}
template<class R> mod_map<R>
mod_map<R>::operator + (const mod_map &m) const
{
assert (impl->from == m.impl->from && impl->to == m.impl->to);
basedvector<linear_combination<R>, 1> v (impl->from->dim ());
for (unsigned i = 1; i <= m.impl->from->dim (); i ++)
v[i] = column (i) + m.column (i);
return mod_map (IMPL, new explicit_map_impl<R> (impl->from, impl->to, v));
}
template<class R> ptr<const free_submodule<R> >
mod_map<R>::kernel () const
{
ptr<const module<R> > from = impl->from,
to = impl->to;
basedvector<linear_combination<R>, 1> from_xs (from->dim ());
for (unsigned i = 1; i <= from->dim (); i ++)
{
linear_combination<R> x (from);
x.muladd (1, i);
from_xs[i] = x;
}
basedvector<linear_combination<R>, 1> to_xs (COPY2, explicit_columns ());
for (unsigned i = 1; i <= to->dim (); i ++)
{
linear_combination<R> from_v (from),
to_v (to);
for (unsigned j = 1; j <= to_xs.size (); j ++)
{
R to_vc = to_v(i);
if (to_vc.is_unit ())
break;
linear_combination<R> &to_x = to_xs[j],
&from_x = from_xs[j];
R to_xc = to_x(i);
if (to_xc == 0)
continue;
if (to_vc == 0 && to_xc != 0)
{
to_v += to_x;
from_v += from_x;
}
else if (! (to_vc | to_xc))
{
tuple<R, R, R> t = to_vc.extended_gcd (to_xc);
assert (get<0> (t) == to_vc*get<1> (t) + get<2> (t)*to_xc);
to_v = get<1> (t)*to_v + get<2> (t)*to_x;
from_v = get<1> (t)*from_v + get<2> (t)*from_x;
assert (to_v(i) != 0);
}
}
R to_vc = to_v(i);
if (to_vc != 0)
{
for (unsigned j = 1; j <= to_xs.size (); j ++)
{
linear_combination<R> &to_x = to_xs[j],
&from_x = from_xs[j];
R to_xc = to_x(i);
if (to_xc != 0)
{
assert (to_vc | to_xc);
R q = to_xc.div (to_vc);
to_x.mulsub (q, to_v);
from_x.mulsub (q, from_v);
}
assert (to_x(i) == 0);
}
}
}
mod_span<R> span (from, from_xs);
return from->submodule (span);
}
template<class R> ptr<const free_submodule<R> >
mod_map<R>::image () const
{
mod_span<R> span (impl->to, explicit_columns ());
return impl->to->submodule (span);
}
template<class R> ptr<const free_submodule<R> >
mod_map<R>::image (basedvector<linear_combination<R>, 1> vs) const
{
mod_span<R> span (impl->from, vs);
ptr<const free_submodule<R> > s = impl->from->submodule (span);
mod_map<R> r = restrict_from (s);
return r.image ();
}
template<class R> ptr<const quotient_module<R> >
mod_map<R>::cokernel () const
{
mod_span<R> span (impl->to, explicit_columns ());
return impl->to->quotient (span);
}
template<class R> ptr<const quotient_module<R> >
mod_map<R>::homology () const
{
ptr<const free_submodule<R> > ker = kernel (),
im = image ();
// display ("ker:", *ker);
// display ("im:", *im);
ptr<const free_submodule<R> > im2 = ker->restrict_submodule (im);
// display ("im2:", *im2);
return ker->quotient (im2);
}
template<class R> basedvector<linear_combination<R>, 1>
mod_map<R>::explicit_columns () const
{
basedvector<linear_combination<R>, 1> v (impl->from->dim ());
for (unsigned i = 1; i <= impl->from->dim (); i ++)
v[i] = column (i);
return v;
}
template<class R>
mod_span<R>::mod_span (ptr<const module<R> > mod,
basedvector<linear_combination<R>, 1> xs0)
{
assert (mod->free_rank () == mod->dim ());
basedvector<linear_combination<R>, 1> xs (COPY2, xs0);
for (unsigned i = 1; i <= mod->dim (); i ++)
{
linear_combination<R> v (mod);
for (unsigned j = 1; j <= xs.size (); j ++)
{
R vc = v(i);
if (vc.is_unit ())
break;
linear_combination<R> &x = xs[j];
R xc = x(i);
if (xc == 0)
continue;
if (vc == 0)
{
v += x;
}
else if (! (vc | xc))
{
tuple<R, R, R> t = vc.extended_gcd (xc);
assert (get<0> (t) == vc*get<1> (t) + get<2> (t)*xc);
v = get<1> (t)*v + get<2> (t)*x;
assert (v(i) != 0);
}
}
R vc = v(i);
if (vc != 0)
{
for (unsigned j = 1; j <= xs.size (); j ++)
{
linear_combination<R> &x = xs[j];
R xc = x(i);
if (xc != 0)
{
assert (vc | xc);
R q = xc.div (vc);
x.mulsub (q, v);
}
assert (x(i) == 0);
}
pivots.append (i);
gens.append (v);
}
}
}
template<class R> void
mod_map<R>::show_self () const
{
printf ("mod_map ");
show (*impl->from);
printf (" -> ");
show (*impl->to);
}
template<class R> void
mod_map<R>::display_self () const
{
show_self (); newline ();
for (unsigned i = 1; i <= impl->from->dim (); i ++)
{
printf (" %d ", i);
impl->from->show_generator (i);
printf (" ");
show (impl->from->generator_grading (i));
printf (": ");
show (column (i));
newline ();
}
}
// ??? io
template<class R> void
writer::write_mod (ptr<const module<R> > m)
{
pair<unsigned &, bool> p = aw->id_io_id.find (m->id);
if (p.second)
{
write_int ((int)p.first);
}
else
{
++ aw->io_id_counter;
unsigned io_id = aw->io_id_counter;
p.first = io_id;
write_int (- (int)io_id);
unsigned n = m->dim (),
r = m->free_rank ();
write_unsigned (n);
write_unsigned (r);
for (unsigned i = 1; i <= n; i ++)
write (*this, m->generator_grading (i));
for (unsigned i = r + 1; i <= n; i ++)
write (*this, m->generator_ann (i));
}
}
template<class R> ptr<const module<R> >
reader::read_mod ()
{
int io_id = read_int ();
if (io_id < 0)
{
unsigned n = read_unsigned ();
unsigned r = read_unsigned ();
basedvector<grading, 1> gr (n);
for (unsigned i = 1; i <= n; i ++)
gr[i] = grading (*this);
basedvector<R, 1> ann (n - r);
for (unsigned i = r + 1; i <= n; i ++)
ann[i - r] = R (*this);
ptr<const module<R> > m = new explicit_module<R> (r, ann, gr);
ar->io_id_id.push ((unsigned)(-io_id), m->id);
module<R>::reader_id_module.push (m->id, m);
return m;
}
else
{
unsigned id = ar->io_id_id(io_id);
return module<R>::reader_id_module(id);
}
}
#endif // _KNOTKIT_ALGEBRA_MODULE_H