From 9bed7282ac90e0c6142b2602e4cadd77f28ce1f7 Mon Sep 17 00:00:00 2001
From: Cotton Seed <cotton@alum.mit.edu>
Date: Tue, 28 Feb 2012 22:08:52 -0500
Subject: [PATCH] Added tensor_product and hom_module to algebra/module.h.

---
 algebra/linear_combination.h |  19 +-
 algebra/module.h             | 366 ++++++++++++++++++++++++++++++++++-
 main.cpp                     |  50 +++++
 3 files changed, 425 insertions(+), 10 deletions(-)

diff --git a/algebra/linear_combination.h b/algebra/linear_combination.h
index 92e62dc..dda27ce 100644
--- a/algebra/linear_combination.h
+++ b/algebra/linear_combination.h
@@ -145,6 +145,17 @@ class linear_combination
   
   unsigned card () const { return v.card (); }
   
+  linear_combination<R> tensor (const linear_combination<R> &lc) const
+  {
+    linear_combination<R> r (m->tensor (lc.m));
+    for (linear_combination_const_iter<R> i = *this; i; i ++)
+      for (linear_combination_const_iter<R> j = lc; j; j ++)
+	{
+	  r.muladd (i.val () * j.val (), m->tensor_generators (i.key (), lc.m, j.key ()));
+	}
+    return r;
+  }
+  
 #ifndef NDEBUG
   void check () const
   {
@@ -292,7 +303,10 @@ linear_combination<R>::show_self () const
       else
 	printf (" + ");
       show (i.val ());
-      printf ("*%d", i.key ());
+      
+      // printf ("*%d", i.key ());
+      printf ("*");
+      m->show_generator (i.key ());
     }
 }
 
@@ -445,7 +459,8 @@ linear_combination<Z2>::show_self () const
 	first = 0;
       else
 	printf ("+");
-      printf ("%d", i.val ());
+      // printf ("%d", i.val ());
+      m->show_generator (i.val ());
     }
 }
 
diff --git a/algebra/module.h b/algebra/module.h
index ce333ed..80aaf59 100644
--- a/algebra/module.h
+++ b/algebra/module.h
@@ -1,15 +1,38 @@
 
 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 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:
+  unsigned id;
+  
+  static unsigned id_counter;
+  
+#if 0
+  static map<basedvector<unsigned, 1>,
+    ptr<const direct_sum<R> > direct_sum_idx;
+#endif
+  
+  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 () { }
+  module ()
+  {
+    id = id_counter;
+    id_counter ++;
+  }
   module (const module &); // doesn't exist
   virtual ~module () { }
   
@@ -28,6 +51,8 @@ class module : public refcounted
   // r < i <= n
   virtual R generator_ann (unsigned i) const = 0;
   
+  bool is_free () const { return dim () == free_rank (); }
+  
   bool is_zero (R c, unsigned i) const
   {
     if (i <= free_rank ())
@@ -53,10 +78,278 @@ class module : public refcounted
   multivariate_laurentpoly<Z> free_poincare_polynomial () const;
   multivariate_laurentpoly<Z> free_delta_poincare_polynomial () const;
   
+  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_tensor_factors (basedvector<ptr<const module<R> >, 1> &factors) const
+  {
+    factors.append (this);
+  }
+  
   void show_self () const;
   void display_self () const;
 };
 
+template<class R> unsigned module<R>::id_counter = 1;
+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;
+
+#if 0
+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 <= terms.size (); i ++)
+      assert (terms[i]->is_free ());
+#endif
+    
+    for (unsigned i = 1; i <= terms.size (); i ++)
+      n += terms[i]->dim ();
+  }
+  
+  unsigned dim () const { return n; }
+  unsigned free_rank () const { return n; }
+  
+  // ???
+  grading generator_grading (unsigned i) const { abort (); }
+  void show_generator (unsigned i) const { abort (); }
+};
+#endif
+
+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>
 {
@@ -287,10 +580,10 @@ template<class R>
 class composition_impl : public map_impl<R>
 {
   // f(g(x))
-  ptr<map_impl<R> > f, g;
+  ptr<const map_impl<R> > f, g;
   
  public:
-  composition_impl (ptr<map_impl<R> > f_, ptr<map_impl<R> > g_)
+  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_)
@@ -304,6 +597,30 @@ class composition_impl : public map_impl<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_->from)),
+      f(f_),
+      g(g_)
+  {
+  }
+  
+  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
 {
@@ -349,9 +666,12 @@ map_builder<R>::init ()
 template<class R>
 class mod_map
 {
-  ptr<map_impl<R> > impl;
+  // ???
+  enum impl_ctor { IMPL };
+
+  ptr<const map_impl<R> > impl;
   
-  mod_map (ptr<map_impl<R> > impl_) : impl(impl_) { }
+  mod_map (impl_ctor, ptr<const map_impl<R> > impl_) : impl(impl_) { }
   
  public:
   mod_map () { }
@@ -389,13 +709,16 @@ class 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->columns (i) != m.impl->columns (i))
+	if (impl->column (i) != m.impl->column (i))
 	  return 0;
       }
     return 1;
@@ -420,7 +743,17 @@ class mod_map
   linear_combination<R> operator [] (unsigned i) const { return impl->column (i); }
   
   linear_combination<R> map (const linear_combination<R> &lc) const { return impl->map (lc); }
-  mod_map compose (const mod_map &m) const { return new composition_impl<R> (impl, m.impl); }
+  mod_map compose (const mod_map &m) const
+  {
+    return mod_map (IMPL,
+		    new composition_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;
@@ -457,6 +790,21 @@ class mod_map
   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
 {
@@ -1319,7 +1667,9 @@ mod_map<R>::display_self () const
   show_self (); newline ();
   for (unsigned i = 1; i <= impl->from->dim (); i ++)
     {
-      printf ("  %d: ", i);
+      printf ("  %d ", i);
+      impl->from->show_generator (i);
+      printf (": ");
       show (column (i));
       newline ();
     }
diff --git a/main.cpp b/main.cpp
index 539c61c..2df7032 100644
--- a/main.cpp
+++ b/main.cpp
@@ -165,6 +165,55 @@ check (const dt_code &dt)
 int
 main ()
 {
+  ptr<const explicit_module<Q> > A
+    = new explicit_module<Q> (2, basedvector<Q, 1> (), basedvector<grading, 1> (2));
+  ptr<const explicit_module<Q> > B
+    = new explicit_module<Q> (3, basedvector<Q, 1> (), basedvector<grading, 1> (3));
+  ptr<const explicit_module<Q> > C
+    = new explicit_module<Q> (3, basedvector<Q, 1> (), basedvector<grading, 1> (3));
+  ptr<const explicit_module<Q> > D
+    = new explicit_module<Q> (2, basedvector<Q, 1> (), basedvector<grading, 1> (2));
+  ptr<const explicit_module<Q> > E
+    = new explicit_module<Q> (2, basedvector<Q, 1> (), basedvector<grading, 1> (3));
+  ptr<const explicit_module<Q> > F
+    = new explicit_module<Q> (2, basedvector<Q, 1> (), basedvector<grading, 1> (2));
+  
+  map_builder<Q> fb (A, B);
+  fb[1].muladd (2, 1);
+  fb[1].muladd (3, 2);
+  fb[2].muladd (-5, 2);
+  fb[2].muladd (4, 3);
+  mod_map<Q> f (fb);
+  display ("f:\n", f);
+  
+  map_builder<Q> gb (C, D);
+  gb[1].muladd (1, 1);
+  gb[2].muladd (3, 1);
+  gb[2].muladd (-2, 2);
+  gb[3].muladd (-6, 2);
+  mod_map<Q> g (gb);
+  display ("g:\n", g);
+
+  map_builder<Q> hb (E, F);
+  hb[1].muladd (3, 2);
+  hb[2].muladd (-3, 1);
+  mod_map<Q> h (hb);
+  display ("h:\n", h);
+  
+  mod_map<Q> fg = f.tensor (g);
+  display ("fg:\n", fg);
+  
+  ptr<const module<Q> > AB_C = (A->tensor (B))->tensor (C),
+    A_BC = A->tensor (B->tensor (C));
+  assert (AB_C == A_BC);
+  
+  assert ((f.tensor (g)).tensor (h) == f.tensor (g.tensor (h)));
+  
+  ptr<const hom_module<Q> > homAB = A->hom (B);
+  linear_combination<Q> x = homAB->map_as_element (f);
+  display ("x:\n", x);
+
+#if 0 
   for (unsigned i = 1; i <= 14; i ++)
     {
       for (unsigned j = 1; j <= mt_links (i, 0); j ++)
@@ -173,6 +222,7 @@ main ()
       for (unsigned j = 1; j <= mt_links (i, 1); j ++)
 	check (mt_link (i, 1, j));
     }
+#endif
   
 #if 0
   knot_diagram kd (rolfsen_knot (8, 19));