diff --git a/src/math/polysat/saturation.cpp b/src/math/polysat/saturation.cpp
index 01a1279ee..63b7f868c 100644
--- a/src/math/polysat/saturation.cpp
+++ b/src/math/polysat/saturation.cpp
@@ -870,19 +870,24 @@ namespace polysat {
         if (s.try_eval(p, val))
             return val == 0 ? N : val.trailing_zeros();
         
-        if (!p.is_var() && p.is_monomial()) {
-            // it's just a product => sum them up
-            dd::pdd_monomial monomial = *p.begin();
-            unsigned parity_sum = monomial.coeff.trailing_zeros();
-            for (pvar c : monomial.vars)
-                parity_sum += min_parity(m.mk_var(c));
-            return std::min(N, parity_sum);
+        unsigned min = 0;
+        if (!p.is_var()) {
+            // parity of a product => sum of parities
+            // parity of sum => minimum of monomial's minimal parities
+            min = UINT32_MAX;
+            for (const auto& monomial : p) {
+                unsigned parity_sum = monomial.coeff.trailing_zeros();
+                for (pvar c : monomial.vars)
+                    parity_sum += min_parity(m.mk_var(c));
+                min = std::min(min, parity_sum);
+            }
         }
+        SASSERT(min <= N);
         
-        for (unsigned j = N; j > 0; --j)
+        for (unsigned j = N; j > min; --j)
             if (is_forced_true(s.parity_at_least(p, j)))
                 return j;
-        return 0;
+        return min;
     }
 
     unsigned saturation::max_parity(pdd const& p) {
@@ -892,19 +897,21 @@ namespace polysat {
         if (s.try_eval(p, val))
             return val == 0 ? N : val.trailing_zeros();
 
+        unsigned max = N;
         if (!p.is_var() && p.is_monomial()) {
             // it's just a product => sum them up
+            // the case of a sum is harder as the lower bound (because of carry bits)
+            // ==> restricted for now to monomials
             dd::pdd_monomial monomial = *p.begin();
-            unsigned parity_sum = monomial.coeff.trailing_zeros();
+            max = monomial.coeff.trailing_zeros();
             for (pvar c : monomial.vars)
-                parity_sum += max_parity(m.mk_var(c));
-            return std::min(N, parity_sum);
+                max += max_parity(m.mk_var(c));
         }
 
-        for (unsigned j = 0; j < N; ++j) 
+        for (unsigned j = 0; j < max; ++j) 
             if (is_forced_true(s.parity_at_most(p, j)))
                 return j;
-        return N;
+        return max;
     }
 
     bool saturation::try_parity(pvar x, conflict& core, inequality const& axb_l_y) {
@@ -1220,7 +1227,7 @@ namespace polysat {
                 m_lemma.reset();
                 m_lemma.insert(~c);
                 m_lemma.insert_eval(~s.eq(y));
-                for (auto& sc_lhs : side_condition_lhs) // TODO: Do we really need the path as a side-condition in case of parity elimination?
+                for (auto& sc_lhs : side_condition_lhs) // the "path to get the parities"
                     m_lemma.insert(sc_lhs);
                 for (auto& sc_rhs : side_condition_rhs)
                     m_lemma.insert(sc_rhs);
diff --git a/src/math/polysat/variable_elimination.cpp b/src/math/polysat/variable_elimination.cpp
index c969437f8..0a563f65e 100644
--- a/src/math/polysat/variable_elimination.cpp
+++ b/src/math/polysat/variable_elimination.cpp
@@ -537,6 +537,32 @@ namespace polysat {
         return id;
     }
     
+    pdd parity_tracker::get_pseudo_inverse(const pdd &p, unsigned parity) {
+        SASSERT(parity < p.power_of_2()); // parity == p.power_of_two() does not make sense. It would mean a * a' = 0
+        LOG("Getting pseudo-inverse of " << p << " for parity " << parity);
+        if (p.is_val())
+            return p.manager().mk_val(get_pseudo_inverse_val(p.val(), parity, p.power_of_2()));
+        
+        unsigned v = get_id(p);
+        if (m_pseudo_inverse.size() <= v)
+            m_pseudo_inverse.setx(v, unsigned_vector(p.power_of_2(), -1), {});
+        
+        if (m_pseudo_inverse[v].size() <= parity)
+            m_pseudo_inverse[v].reserve(p.power_of_2(), -1);
+        
+        pvar p_inv = m_pseudo_inverse[v][parity];
+        
+        if (p_inv == -1) {
+            p_inv = s.add_var(p.power_of_2());
+            m_pseudo_inverse[v][parity] = p_inv;
+            // NB: Strictly speaking this condition does not say that "p_inv" is the pseudo-inverse.
+            // However, the variable elimination lemma stays valid even if "p_inv" is not really the pseudo-inverse anymore (e.g., because of parity was reduced)
+            s.add_clause(~s.parity_at_most(p, parity), s.eq(p * s.var(p_inv), rational::power_of_two(parity)), true); // TODO: Depends on the definition of redundancy
+            verbose_stream() << "Pseudo-Inverse v" << p_inv << " of " << p << " introduced\n";
+        }
+        return s.var(p_inv);
+    }
+    
     pdd parity_tracker::get_inverse(const pdd &p) {
         LOG("Getting inverse of " << p);
         if (p.is_val()) {
@@ -550,7 +576,8 @@ namespace polysat {
             return s.var(m_inverse[v]);
         
         pvar inv = s.add_var(p.power_of_2());
-        pdd inv_pdd = p.manager().mk_var(inv);
+        verbose_stream() << "Inverse v" << inv << " of " << p << " introduced\n";
+        pdd inv_pdd = s.var(inv);
         m_inverse.setx(v, inv, -1);
         s.add_clause(s.eq(inv_pdd * p, p.manager().one()), false);
         return inv_pdd;
@@ -576,6 +603,7 @@ namespace polysat {
         }
         else {
             odd_v = s.add_var(p.power_of_2());
+            verbose_stream() << "Odd part v" << odd_v << " of " << p << " introduced\n";
             m_odd.setx(v, optional<std::pair<pvar, bool_vector>>({ odd_v, bool_vector(p.power_of_2(), false) }), optional<std::pair<pvar, bool_vector>>::undef());
         }
         
@@ -593,12 +621,14 @@ namespace polysat {
                 path.push_back(~c);
                 if (needs_propagate)
                     m_builder.insert(~c);
+                verbose_stream() << "Side-condition: " << ~c << "\n";
             }
             else {
                 upper = middle;
                 path.push_back(c);
                 if (needs_propagate)
                     m_builder.insert(c);
+                verbose_stream() << "Side-condition: " << c << "\n";
             }
             LOG("Its in [" << lower << "; " << upper << ")");
         }
@@ -638,29 +668,47 @@ namespace polysat {
         if (is_multiple == l_true) // we multiply with a factor to make them equal
             return { b1 - b * mul_fac, true, {} };
         
-        #if 1
+#if 1
         return { p, false, {} };
-        #else
+#else
+        
+        if (!a.is_monomial() || !a1.is_monomial())
+            return { p , false, {} };
+        
         if (!a1.is_var() && !a1.is_val()) {
+            // TODO: Compromise: Maybe only monomials...? Does this make sense?
             //return { p, false, {} };
-            LOG("Warning: Inverting " << a1 << " although it is not a single variable - might not be a good idea"); // TODO: Compromise: Maybe only monomials...?
+            LOG("Warning: Inverting " << a1 << " although it is not a single variable - might not be a good idea");
         }
         if (!a.is_var() && !a.is_val()) {
             //return { p, false, {} };
             LOG("Warning: Inverting " << a << " although it is not a single variable - might not be a good idea");
         }
         
-        if (!a.is_monomial() || !a1.is_monomial())
-            return { p , false, {} };
-        
         // We don't know whether it will work. Use the parity of the assignment
         
+#if 1
+        unsigned a_parity;
+        if ((a_parity = saturation.min_parity(a)) != saturation.max_parity(a) || saturation.min_parity(a1) < a_parity)
+            return { p, false, {} }; // We need the parity of a and this has to be for sure less than the parity of a1
+        
+        svector<signed_constraint> precondition;
+        pdd a_pi = get_pseudo_inverse(a, a_parity);
+        
+        if (a_parity > 0) {
+            pdd shift = s.lshr(a1, a1.manager().mk_val(a_parity));
+            precondition.push_back(s.eq(rational::power_of_two(a_parity) * shift, a1)); // TODO: Or s.parity_at_least(a1, a_parity) but we want to reuse the variable introduced by the shift
+            return { a_pi * (-b) * shift + b1, true, {std::move(precondition)} };
+        }
+        // Special case: If it is already odd we can directly use the pseudo inverse (as it is the inverse in this case!)
+        return { a_pi * (-b) * a + b1, true, {std::move(precondition)} };
+#else
         unsigned a_parity;
         unsigned a1_parity;
-        
+            
         if ((a_parity = saturation.min_parity(a)) != saturation.max_parity(a) || (a1_parity = saturation.min_parity(a1)) != saturation.max_parity(a1))
             return { p, false, {} }; // We need the parity, but we failed to get it precisely
-            
+        
         if (a_parity > a1_parity) {
             SASSERT(false); // get_multiple should have excluded this case already
             return { p, false, {} };
@@ -673,9 +721,10 @@ namespace polysat {
         pdd inv_odd_a = get_inverse(odd_a);
         
         LOG("Forced elimination: " << odd_a1 * inv_odd_a * rational::power_of_two(a1_parity - a_parity) * b + b1);
-        verbose_stream() << "Forced elimination: " << odd_a1 * inv_odd_a * rational::power_of_two(a1_parity - a_parity) * b + b1 << "\n";
-        verbose_stream() << "From: " << "eliminated v" << x << " with a = " << a << "; b = " << b << "; p = " << p << "\n";
-        return { odd_a1 * inv_odd_a * rational::power_of_two(a1_parity - a_parity) * b + b1, true, {std::move(precondition)} };
+        verbose_stream() << "Forced elimination: " << odd_a1 * inv_odd_a * rational::power_of_two(a1_parity - a_parity) * (-b) + b1 << "\n";
+        verbose_stream() << "From: " << "eliminated v" << x << " with a = " << a << "; -b = " << -b << "; p = " << p << "\n";
+        return { odd_a1 * inv_odd_a * rational::power_of_two(a1_parity - a_parity) * (-b) + b1, true, {std::move(precondition)} };
+#endif
 #endif
     }
 }
diff --git a/src/math/polysat/variable_elimination.h b/src/math/polysat/variable_elimination.h
index 3b1bde455..ee403b1b3 100644
--- a/src/math/polysat/variable_elimination.h
+++ b/src/math/polysat/variable_elimination.h
@@ -62,6 +62,7 @@ namespace polysat {
         
         vector<optional<std::pair<pvar, bool_vector>>> m_odd;
         unsigned_vector m_inverse;
+        vector<unsigned_vector> m_pseudo_inverse;
         
         struct optional_pdd_hash {
             unsigned operator()(optional<pdd> const& args) const {
@@ -72,12 +73,19 @@ namespace polysat {
         
         pdd_to_id m_pdd_to_id; // if we want to use arbitrary pdds instead of pvars
         
+        rational get_pseudo_inverse_val(const rational &p, unsigned parity, unsigned power_of_two) {
+            rational iv;
+            VERIFY((p / rational::power_of_two(parity)).mult_inverse(power_of_two - parity, iv));
+            return iv;
+        }
+        
         unsigned get_id(const pdd& p);
         
     public:
         
         parity_tracker(solver& s) : s(s), m_builder(s) {}
         
+        pdd get_pseudo_inverse(const pdd& p, unsigned parity);
         pdd get_inverse(const pdd& p);
         pdd get_odd(const pdd& p, unsigned parity, svector<signed_constraint>& pat);