Q uses gmp. Deleted README2. Z::muladdeq uses mpz_add_mul.

Makefile

@@ -17,7 +17,7 @@ CXXFLAGS = $(OPTFLAGS) -Wall -Wno-unused $(INCLUDES)
 
 LIB_OBJS = lib/refcount.o \
   lib/lib.o lib/smallbitset.o lib/bitset.o lib/setcommon.o lib/io.o lib/directed_multigraph.o
-ALGEBRA_OBJS = algebra/algebra.o algebra/Q.o algebra/grading.o algebra/polynomial.o
+ALGEBRA_OBJS = algebra/algebra.o algebra/grading.o algebra/polynomial.o
 KNOTKIT_OBJS = planar_diagram.o dt_code.o knot_diagram.o cube.o spanning_tree_complex.o \
   smoothing.o cobordism.o knot_tables.o sseq.o \
   knot_parser/knot_parser.o knot_parser/knot_scanner.o \
@@ -34,8 +34,7 @@ ALGEBRA_HEADERS = algebra/algebra.h algebra/grading.h algebra/module.h \
   algebra/Z2.h algebra/linear_combination.h \
   algebra/Z.h algebra/Zp.h algebra/Q.h \
   algebra/polynomial.h algebra/multivariate_polynomial.h \
-  algebra/multivariate_laurentpoly.h \
-  algebra/laurentpoly.h algebra/fraction_field.h
+  algebra/multivariate_laurentpoly.h algebra/fraction_field.h
 KNOTKIT_HEADERS = knotkit.h planar_diagram.h dt_code.h knot_diagram.h \
   smoothing.h cobordism.h cube.h spanning_tree_complex.h cube_impl.h sseq.h
 

README2

@@ -1,3 +0,0 @@
-Read this!
-
-and make another change

algebra/Q.cpp

@@ -1,19 +0,0 @@
-
-#include <algebra/algebra.h>
-
-void
-Q::reduce ()
-{
-  unsigned c = unsigned_gcd ((unsigned)std::abs (n), d);
-  n /= (int)c;
-  d /= c;
-}
-
-void
-Q::show_self () const
-{
-  if (d == 1)
-    printf ("%d", n);
-  else
-    printf ("%d/%u", n, d);
-}

algebra/Q.h

@@ -3,61 +3,182 @@ class Q
 {
  public:
   typedef ::linear_combination<Q> linear_combination;
-  // typedef linear_combination_iter<Q> linear_combination_iter;
   typedef ::linear_combination_const_iter<Q> linear_combination_const_iter;
-  
+
  private:
-  int n;
-  unsigned d;
+  enum steal { STEAL };
   
-  void reduce ();
+  class Q_impl : public refcounted
+  {
+  public:
+    mpq_t x;
+    
+  public:
+    Q_impl () { mpq_init (x); }
+    Q_impl (int init)
+    {
+      mpq_init (x);
+      mpq_set_si (x, init, 1);
+    }
+    
+    Q_impl (copy, mpq_srcptr init)
+    {
+      mpq_init (x);
+      mpq_set (x, init);
+    }
+    
+    Q_impl (steal, mpq_srcptr init) { x[0] = *init; }
+    Q_impl (reader &r)
+    {
+      mpq_init (x);
+      mpz_inp_raw (mpq_numref (x), r.fp);
+      mpz_inp_raw (mpq_denref (x), r.fp);
+    }
+    
+    ~Q_impl () { mpq_clear (x); }
+    
+    void write_self (writer &w) const
+    {
+      mpz_out_raw (w.fp, mpq_numref (x));
+      mpz_out_raw (w.fp, mpq_denref (x));
+    }
+  };
   
+  ptr<Q_impl> impl;
+  
+  Q (steal, mpq_srcptr init) : impl(new Q_impl (STEAL, init)) { }
+
  public:
-  Q () : n(0), d(1) { }
-  Q (int x) : n(x), d(1) { }
-  Q (int n_, unsigned d_) : n(n_), d(d_) { assert (d != 0); reduce (); }
-  Q (const Q &q) : n(q.n), d(q.d) { }
+  Q () : impl(new Q_impl) { }
+  Q (int init) : impl(new Q_impl (init)) { }
+  Q (const Q &q) : impl(q.impl) { }
+  Q (copy, const Q &q) : impl(new Q_impl (COPY, q.impl->x)) { }
+  Q (reader &r) : impl(new Q_impl (r)) { }
   ~Q () { }
   
-  Q &operator = (const Q &q) { n = q.n; d = q.d; return *this; }
-  Q &operator = (int x) { n = x; d = 1; return *this; }
+  Q &operator = (const Q &q) { impl = q.impl; return *this; }
+  Q &operator = (int x) { impl = new Q_impl (x); return *this; }
   
-  bool operator == (const Q &q) const { return n == q.n && d == q.d; }
+  bool operator == (const Q &q) const { return mpq_cmp (impl->x,q.impl->x) == 0; }
   
-  bool operator == (int x) const { return n == x && d == 1; }
-  bool operator != (int x) const { return !operator == (x); }
+  bool operator == (int r) const { return mpq_cmp_si (impl->x, r, 1) == 0; }
+  bool operator != (int r) const { return !operator == (r); }
   
-  bool is_unit () const { return n != 0; }
+  bool operator < (const Q &q) const { return mpq_cmp (impl->x, q.impl->x) < 0; }
   
-  Q operator + (const Q &x) const { return Q (n*x.d + x.n*d, d*x.d); }
-  Q operator - (const Q &x) const { return Q (n*x.d - x.n*d, d*x.d); }
-  Q operator * (const Q &x) const { return Q (n*x.n, d*x.d); }
-  Q operator / (const Q &x) const { return operator * (x.recip ()); }
-  
-  Q recip () const
+  bool is_unit () const
   {
-    assert (is_unit ());
-    if (n < 0)
-      return Q (-d, -n);
-    else
-      return Q (d, n);
+    return *this != 0;
   }
   
-  Q &operator += (const Q &x) { n = n*x.d + x.n*d; d *= x.d; reduce (); return *this; }
-  Q &operator -= (const Q &x) { n = n*x.d - x.n*d; d *= x.d; reduce (); return *this; }
-  Q &operator *= (const Q &x) { n *= x.n; d *= x.d; reduce (); return *this; }
-  Q &operator /= (const Q &x) { return operator *= (x.recip ()); }
+  Q operator - () const
+  {
+    mpq_t x;
+    mpq_init (x);
+    mpq_neg (x, impl->x);
+    return Q (STEAL, x);
+  }
+
+  Q recip () const
+  {
+    mpq_t x;
+    mpq_init (x);
+    mpq_inv (x, impl->x);
+    return Q (STEAL, x);
+  }
   
-  bool divides (const Q &q) const { return (n != 0) || (q.n == 0); }
-  bool operator | (const Q &q) const { return divides (q); }
+  Q operator + (const Q &q) const
+  {
+    mpq_t x;
+    mpq_init (x);
+    mpq_add (x, impl->x, q.impl->x);
+    return Q (STEAL, x);
+  }
+  
+  Q operator - (const Q &q) const
+  {
+    mpq_t x;
+    mpq_init (x);
+    mpq_sub (x, impl->x, q.impl->x);
+    return Q (STEAL, x);
+  }
+  
+  Q operator * (const Q &q) const
+  {
+    mpq_t x;
+    mpq_init (x);
+    mpq_mul (x, impl->x, q.impl->x);
+    return Q (STEAL, x);
+  }
+  
+  Q operator / (const Q &q) const
+  {
+    assert (q != 0);
+    
+    mpq_t x;
+    mpq_init (x);
+    mpq_div (x, impl->x, q.impl->x);
+    return Q (STEAL, x);
+  }
+  
+  Q &muladdeq (const Q &q1, const Q &q2)
+  {
+    // ??? do inline saves refcount overhead
+    return operator += (q1 * q2);
+  }
+  
+  Q &operator += (const Q &q)
+  {
+    mpq_add (impl->x, impl->x, q.impl->x);
+    return *this;
+  }
+  
+  Q &operator -= (const Q &q)
+  {
+    mpq_sub (impl->x, impl->x, q.impl->x);
+    return *this;
+  }
+  
+  Q &operator *= (const Q &q)
+  {
+    mpq_mul (impl->x, impl->x, q.impl->x);
+    return *this;
+  }
+  
+  Q &operator /= (const Q &q)
+  {
+    assert (q != 0);
+    
+    mpq_div (impl->x, impl->x, q.impl->x);
+    return *this;
+  }
+  
+  // d | n, d.divides (n)
+  bool divides (const Q &num) const
+  {
+    return *this != 0 || num == 0;
+  }
+  
+  bool operator | (const Q &num) const { return divides (num); }
+  
+  Q div (const Q &d) const { return operator / (d); }
+  
+  triple<Q, Q, Q> extended_gcd (const Q &q) const
+  {
+    if (*this != 0)
+      return triple<Q, Q, Q> (*this, 1, 0);
+    else
+      return triple<Q, Q, Q> (q, 0, 1);
+  }
   
   Q gcd (const Q &q) const
   {
-    assert (n != 0 && q.n != 0);
-    return Q (1);
+    assert (*this != 0 || q != 0);
+    return 1;
   }
   
   static void show_ring () { printf ("Q"); }
-  void show_self () const;
+  void show_self () const { mpq_out_str (stdout, 10, impl->x); }
   void display_self () const { show_self (); newline (); }
+  void write_self (writer &w) const { write (w, *impl); }
 };

algebra/Z.h

@@ -5,9 +5,9 @@ class Z
   typedef ::linear_combination<Z> linear_combination;
   typedef ::linear_combination_const_iter<Z> linear_combination_const_iter;
   
+ private:
   enum steal { STEAL };
   
- private:
   class Z_impl : public refcounted
   {
   public:
@@ -114,8 +114,8 @@ class Z
   
   Z &muladdeq (const Z &z1, const Z &z2)
   {
-    // ??? use muladd primitive
-    return operator += (z1 * z2);
+    mpz_addmul (impl->x, z1.impl->x, z2.impl->x);
+    return *this;
   }
   
   Z &operator += (const Z &z)

main.cpp

@@ -111,6 +111,7 @@ test_field ()
 int
 main ()
 {
+#if 0
   knot_diagram kd (rolfsen_knot (8, 19));
   cube<Z2> c (kd);
   sseq ss = compute_szabo_sseq (c);
@@ -140,15 +141,18 @@ main ()
     
     assert (q == p*p);
   }
+#endif
   
-#if 0
+#if 1
   test_ring<Z2> (2);
   test_ring<Z> (0);
+  test_ring<Q> (0);
   test_ring<Zp<2> > (2);
   test_ring<Zp<3> > (3);
   test_ring<Zp<5> > (5);
   test_ring<Zp<7> > (7);
   
+  test_field<Q> ();
   test_field<Zp<7> > ();
   test_field<Zp<5> > ();
   test_field<Zp<3> > ();

todo.txt

@@ -1,5 +1,6 @@
 in knotkit/
- - twisted theory over U
+ - my (or Kh) tangle algebra over Z
+ - enumerate 2-connect sum knots to test mutation invariance over Z
  - add homology orientations for odd Khovanov homology and spectral
    sequence over Z
  - unify smoothing, resolution_diagram wiring (and knot_diagram?)
@@ -16,9 +17,6 @@ in lib/
  - add make_pair, etc.
 
 in algebra/
- - Q should use gmp
- - laurentpoly interface needs to be cleaned up
- - add Zp
  - linear_combnation can be more efficient (eg no searching when
    coefficients die)
  - monomial ideals and maybe groebner bases