#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