diff --git a/src/sat/smt/xor_gaussian.cpp b/src/sat/smt/xor_gaussian.cpp
index d2da78d4b..5c8d5f245 100644
--- a/src/sat/smt/xor_gaussian.cpp
+++ b/src/sat/smt/xor_gaussian.cpp
@@ -358,16 +358,33 @@ bool EGaussian::full_init(bool& created) {
     return !inconsistent();
 }
 
-static void print_matrix(ostream& out, PackedMatrix& mat) {
-    for (unsigned rowIdx = 0; rowIdx < mat.num_rows(); rowIdx++) {
-        const PackedRow& row = mat[rowIdx];
-        for(int i = 0; i < row.get_size() * 64; i++) {
+std::ostream& PackedMatrix::display_dense(std::ostream& out) const {
+    for (unsigned rowIdx = 0; rowIdx < num_rows(); rowIdx++) {
+        const PackedRow& row = (*this)[rowIdx];
+        for(int i = 0; i < row.get_size() * 64; i++) 
             out << (int)row[i];
-        }
         out << " -- rhs: " << row.rhs() << " -- row: " << rowIdx << "\n";
     }
+    return out;
 }
 
+std::ostream& PackedMatrix::display_sparse(std::ostream& out) const {
+    for (auto const& row : *this) {
+        bool first = true;
+        for (int i = 0; i < row.get_size() * 64; ++i) {
+            if (row[i]) {
+                if (first && row.rhs())
+                    out << -i-1 << " ";
+                else
+                    out << i+1 << " ";
+                first = false;
+            }
+        }
+    }
+    return out;
+}
+
+
 // Applies Gaussian-Jordan elimination (search level). This function does not add conflicts/propagate/... Just reduces the matrix
 void EGaussian::eliminate() {
     SASSERT(m_solver.s().at_search_lvl());
@@ -376,10 +393,9 @@ void EGaussian::eliminate() {
     unsigned row_i = 0;
     unsigned col = 0;
 
-    print_matrix(std::cout, m_mat);
-    std::cout << std::endl;
-    TRACE("xor", print_matrix(tout, m_mat));
-
+    m_mat.display_dense(std::cout) << std::endl;
+    TRACE("xor", m_mat.display_dense(tout) << "\n");
+    
     // Gauss-Jordan Elimination
     while (row_i != m_num_rows && col != m_num_cols) {
         unsigned row_with_1_in_col = row_i;
@@ -409,9 +425,8 @@ void EGaussian::eliminate() {
             row_i++;
         }
         col++;
-        TRACE("xor", print_matrix(tout, m_mat));
-        print_matrix(std::cout, m_mat);
-        std::cout << std::endl;
+        TRACE("xor", m_mat.display_dense(tout) << "\n");
+        m_mat.display_dense(std::cout) << "\n";        
     }
     std::cout << "-------------" << std::endl;
 }
@@ -811,10 +826,10 @@ void EGaussian::eliminate_column(unsigned p, gauss_data& gqd) {
     unsigned row_i = 0;
 
     elim_called++;
-    
-    print_matrix(std::cout, m_mat);
-    std::cout << std::endl;
-    TRACE("xor", print_matrix(tout, m_mat));
+
+
+    m_mat.display_dense(std::cout) << "\n";
+    TRACE("xor", m_mat.display_dense(tout) << "\n");
     
     while (row_i < row_size) {
         //Row has a '1' in eliminating column, and it's not the row responsible
@@ -945,10 +960,9 @@ void EGaussian::eliminate_column(unsigned p, gauss_data& gqd) {
             }
         }
         row_i++;
-        
-        print_matrix(std::cout, m_mat);
-        std::cout << std::endl;
-        TRACE("xor", print_matrix(tout, m_mat));
+
+        m_mat.display_dense(std::cout) << "\n";
+        TRACE("xor", m_mat.display_dense(tout) << "\n");        
     }
 
 }
@@ -1082,3 +1096,5 @@ void EGaussian::move_back_xor_clauses() {
     for (const auto& x: m_xorclauses) 
         m_solver.m_xorclauses.push_back(std::move(x));    
 }
+
+
diff --git a/src/sat/smt/xor_gaussian.h b/src/sat/smt/xor_gaussian.h
index 1635feec8..358b89fa7 100644
--- a/src/sat/smt/xor_gaussian.h
+++ b/src/sat/smt/xor_gaussian.h
@@ -530,11 +530,11 @@ namespace xr {
             }
         };
     
-        iterator begin() {
+        iterator begin() const {
             return iterator(mp, numCols);
         }
     
-        iterator end() {
+        iterator end() const {
             return iterator(mp + numRows * (numCols + 1), numCols);
         }
         
@@ -545,6 +545,10 @@ namespace xr {
         unsigned num_cols() const {
             return numCols;
         }
+
+        std::ostream& display_dense(std::ostream& out) const;
+        
+        std::ostream& display_sparse(std::ostream& out) const;
         
     private:
     
@@ -609,6 +613,8 @@ namespace xr {
             }
             out << std::endl;
         }
+
+        std::ostream& display(std::ostream& out) const { return m_mat.display_sparse(out); }
     
     private:
         xr::solver& m_solver;   // original sat solver
diff --git a/src/sat/smt/xor_solver.cpp b/src/sat/smt/xor_solver.cpp
index 9e3f8a9e2..06b601073 100644
--- a/src/sat/smt/xor_solver.cpp
+++ b/src/sat/smt/xor_solver.cpp
@@ -310,9 +310,8 @@ namespace xr {
     }
 
     std::ostream& solver::display(std::ostream& out) const {
-        out << "xor clauses: " << m_xorclauses.size() << "\n";
-        for (auto const& x : m_xorclauses)
-            out << x << "\n";
+        for (auto* mat : m_gmatrices) 
+            mat->display(out);        
         return out;
     }