From 8da9850d45b2fcff50b12743c610c0c249ad9948 Mon Sep 17 00:00:00 2001
From: Jakob Rath <git@jakobrath.eu>
Date: Wed, 21 Dec 2022 12:13:05 +0100
Subject: [PATCH] Add rational::pseudo_inverse

---
 src/util/mpq.h        |  8 ++++++++
 src/util/rational.cpp | 18 ++++++++++++++++++
 src/util/rational.h   |  7 +++++++
 3 files changed, 33 insertions(+)

diff --git a/src/util/mpq.h b/src/util/mpq.h
index f1afcbe50..7d5a6795a 100644
--- a/src/util/mpq.h
+++ b/src/util/mpq.h
@@ -490,6 +490,8 @@ public:
 
     void machine_div_rem(mpz const & a, mpz const & b, mpz & c, mpz & d) { mpz_manager<SYNCH>::machine_div_rem(a, b, c, d); }
 
+    void machine_div2k(mpz const & a, unsigned k, mpz & c) { mpz_manager<SYNCH>::machine_div2k(a, k, c); }
+
     void div(mpz const & a, mpz const & b, mpz & c) { mpz_manager<SYNCH>::div(a, b, c); }
     
     void rat_div(mpz const & a, mpz const & b, mpq & c) {
@@ -516,6 +518,12 @@ public:
         machine_div(a.m_num, b.m_num, c);
     }
 
+    void machine_idiv2k(mpq const & a, unsigned k, mpq & c) {
+        SASSERT(is_int(a));
+        machine_div2k(a.m_num, k, c.m_num);
+        reset_denominator(c);
+    }
+
     void idiv(mpq const & a, mpq const & b, mpq & c) {
         SASSERT(is_int(a) && is_int(b));
         div(a.m_num, b.m_num, c.m_num);
diff --git a/src/util/rational.cpp b/src/util/rational.cpp
index af3c89ced..531669dd4 100644
--- a/src/util/rational.cpp
+++ b/src/util/rational.cpp
@@ -153,3 +153,21 @@ bool rational::mult_inverse(unsigned num_bits, rational & result) const {
     return true;
 }
 
+/**
+ * Compute multiplicative pseudo-inverse modulo 2^num_bits:
+ *
+ *      mod(n * n.pseudo_inverse(bits), 2^bits) == 2^k,
+ *      where k is maximal such that 2^k divides n.
+ *
+ * Precondition: number is non-zero.
+ */
+rational rational::pseudo_inverse(unsigned num_bits) const {
+    rational result;
+    rational const& n = *this;
+    SASSERT(!n.is_zero());  // TODO: or we define pseudo-inverse of 0 as 0.
+    unsigned const k = n.trailing_zeros();
+    rational const odd = machine_div2k(n, k);
+    VERIFY(odd.mult_inverse(num_bits, result));
+    SASSERT_EQ(mod(n * result, rational::power_of_two(num_bits)), rational::power_of_two(k));
+    return result;
+}
diff --git a/src/util/rational.h b/src/util/rational.h
index 2ec822bb4..e45fbd4e3 100644
--- a/src/util/rational.h
+++ b/src/util/rational.h
@@ -229,6 +229,12 @@ public:
         return r;
     }
 
+    friend inline rational machine_div2k(rational const & r1, unsigned k) {
+        rational r;
+        rational::m().machine_idiv2k(r1.m_val, k, r.m_val);
+        return r;
+    }
+
     friend inline rational mod(rational const & r1, rational const & r2) {
         rational r;
         rational::m().mod(r1.m_val, r2.m_val, r.m_val);
@@ -355,6 +361,7 @@ public:
     }
 
     bool mult_inverse(unsigned num_bits, rational & result) const;
+    rational pseudo_inverse(unsigned num_bits) const;
 
     static rational const & zero() {
         return m_zero;