diff --git a/src/crab/interval.cpp b/src/crab/interval.cpp
index 8db510e2a..f37af40e7 100644
--- a/src/crab/interval.cpp
+++ b/src/crab/interval.cpp
@@ -4,6 +4,29 @@
 
 namespace crab {
 
+static interval_t make_dividend_when_both_nonzero(const interval_t& dividend, const interval_t& divisor) {
+    if (dividend.ub() >= 0) {
+        return dividend;
+    }
+    if (divisor.ub() < 0) {
+        return dividend + divisor + interval_t{1};
+    }
+    return dividend + interval_t{1} - divisor;
+}
+
+interval_t interval_t::operator*(const interval_t& x) const {
+    if (is_bottom() || x.is_bottom()) {
+        return bottom();
+    }
+    const auto [clb, cub] = std::minmax({
+        _lb * x._lb,
+        _lb * x._ub,
+        _ub * x._lb,
+        _ub * x._ub,
+    });
+    return interval_t{clb, cub};
+}
+
 // Signed division. eBPF has no instruction for this.
 interval_t interval_t::operator/(const interval_t& x) const {
     if (is_bottom() || x.is_bottom()) {
@@ -14,38 +37,37 @@ interval_t interval_t::operator/(const interval_t& x) const {
         //   the linear interval solver can perform many divisions where
         //   the divisor is a singleton interval. We optimize for this case.
         number_t c = *n;
-        if (c == number_t{1}) {
+        if (c == 1) {
             return *this;
-        } else if (c > number_t{0}) {
-            return interval_t(_lb / bound_t{c}, _ub / bound_t{c});
-        } else if (c < number_t{0}) {
-            return interval_t(_ub / bound_t{c}, _lb / bound_t{c});
+        } else if (c > 0) {
+            return interval_t{_lb / c, _ub / c};
+        } else if (c < 0) {
+            return interval_t{_ub / c, _lb / c};
         } else {
             // The eBPF ISA defines division by 0 as resulting in 0.
-            return interval_t(number_t(0));
+            return interval_t{0};
         }
     }
     if (x.contains(0)) {
         // The divisor contains 0.
-        interval_t l(x._lb, bound_t(-1));
-        interval_t u(bound_t(1), x._ub);
-        return (operator/(l) | operator/(u) | interval_t(number_t(0)));
+        interval_t l{x._lb, -1};
+        interval_t u{1, x._ub};
+        return operator/(l) | operator/(u) | interval_t{0};
     } else if (contains(0)) {
         // The dividend contains 0.
-        interval_t l(_lb, bound_t(-1));
-        interval_t u(bound_t(1), _ub);
-        return ((l / x) | (u / x) | interval_t(number_t(0)));
+        interval_t l{_lb, -1};
+        interval_t u{1, _ub};
+        return (l / x) | (u / x) | interval_t{0};
     } else {
         // Neither the dividend nor the divisor contains 0
-        interval_t a =
-            (_ub < number_t{0})
-                ? (*this + ((x._ub < number_t{0}) ? (x + interval_t(number_t(1))) : (interval_t(number_t(1)) - x)))
-                : *this;
-        bound_t ll = a._lb / x._lb;
-        bound_t lu = a._lb / x._ub;
-        bound_t ul = a._ub / x._lb;
-        bound_t uu = a._ub / x._ub;
-        return interval_t(std::min({ll, lu, ul, uu}), std::max({ll, lu, ul, uu}));
+        interval_t a = make_dividend_when_both_nonzero(*this, x);
+        const auto [clb, cub] = std::minmax({
+            a._lb / x._lb,
+            a._lb / x._ub,
+            a._ub / x._lb,
+            a._ub / x._ub,
+        });
+        return interval_t{clb, cub};
     }
 }
 
@@ -63,34 +85,33 @@ interval_t interval_t::SDiv(const interval_t& x) const {
             if (c == 1) {
                 return *this;
             } else if (c != 0) {
-                return interval_t(_lb / bound_t{c}, _ub / bound_t{c});
+                return interval_t{_lb / c, _ub / c};
             } else {
                 // The eBPF ISA defines division by 0 as resulting in 0.
-                return interval_t(number_t(0));
+                return interval_t{0};
             }
         }
     }
     if (x.contains(0)) {
         // The divisor contains 0.
-        interval_t l(x._lb, bound_t(-1));
-        interval_t u(bound_t(1), x._ub);
-        return (SDiv(l) | SDiv(u) | interval_t(number_t(0)));
+        interval_t l{x._lb, -1};
+        interval_t u{1, x._ub};
+        return SDiv(l) | SDiv(u) | interval_t{0};
     } else if (contains(0)) {
         // The dividend contains 0.
-        interval_t l(_lb, bound_t(-1));
-        interval_t u(bound_t(1), _ub);
-        return (l.SDiv(x) | u.SDiv(x) | interval_t(number_t(0)));
+        interval_t l{_lb, -1};
+        interval_t u{1, _ub};
+        return l.SDiv(x) | u.SDiv(x) | interval_t{0};
     } else {
         // Neither the dividend nor the divisor contains 0
-        interval_t a =
-            (_ub < number_t{0})
-                ? (*this + ((x._ub < number_t{0}) ? (x + interval_t(number_t(1))) : (interval_t(number_t(1)) - x)))
-                : *this;
-        bound_t ll = a._lb / x._lb;
-        bound_t lu = a._lb / x._ub;
-        bound_t ul = a._ub / x._lb;
-        bound_t uu = a._ub / x._ub;
-        return interval_t(std::min({ll, lu, ul, uu}), std::max({ll, lu, ul, uu}));
+        interval_t a = make_dividend_when_both_nonzero(*this, x);
+        const auto [clb, cub] = std::minmax({
+            a._lb / x._lb,
+            a._lb / x._ub,
+            a._ub / x._lb,
+            a._ub / x._ub,
+        });
+        return interval_t{clb, cub};
     }
 }
 
@@ -107,36 +128,35 @@ interval_t interval_t::UDiv(const interval_t& x) const {
             number_t c{n->cast_to<uint64_t>()};
             if (c == 1) {
                 return *this;
-            } else if (c > number_t{0}) {
-                return interval_t(_lb.UDiv(bound_t{c}), _ub.UDiv(bound_t{c}));
+            } else if (c > 0) {
+                return interval_t{_lb.UDiv(c), _ub.UDiv(c)};
             } else {
                 // The eBPF ISA defines division by 0 as resulting in 0.
-                return interval_t(number_t(0));
+                return interval_t{0};
             }
         }
     }
     if (x.contains(0)) {
         // The divisor contains 0.
-        interval_t l(x._lb, bound_t(-1));
-        interval_t u(bound_t(1), x._ub);
-        return UDiv(l) | UDiv(u) | interval_t(number_t(0));
+        interval_t l{x._lb, -1};
+        interval_t u{1, x._ub};
+        return UDiv(l) | UDiv(u) | interval_t{0};
     }
     if (contains(0)) {
         // The dividend contains 0.
-        interval_t l(_lb, bound_t(-1));
-        interval_t u(bound_t(1), _ub);
-        return l.UDiv(x) | u.UDiv(x) | interval_t(number_t(0));
+        interval_t l{_lb, -1};
+        interval_t u{1, _ub};
+        return l.UDiv(x) | u.UDiv(x) | interval_t{0};
     }
     // Neither the dividend nor the divisor contains 0
-    interval_t a =
-        (_ub < number_t{0})
-            ? (*this + ((x._ub < number_t{0}) ? (x + interval_t(number_t(1))) : (interval_t(number_t(1)) - x)))
-            : *this;
-    bound_t ll = a._lb.UDiv(x._lb);
-    bound_t lu = a._lb.UDiv(x._ub);
-    bound_t ul = a._ub.UDiv(x._lb);
-    bound_t uu = a._ub.UDiv(x._ub);
-    return interval_t(std::min({ll, lu, ul, uu}), std::max({ll, lu, ul, uu}));
+    interval_t a = make_dividend_when_both_nonzero(*this, x);
+    const auto [clb, cub] = std::minmax({
+        a._lb.UDiv(x._lb),
+        a._lb.UDiv(x._ub),
+        a._ub.UDiv(x._lb),
+        a._ub.UDiv(x._ub),
+    });
+    return interval_t{clb, cub};
 }
 
 // Signed remainder (modulo).
@@ -149,37 +169,37 @@ interval_t interval_t::SRem(const interval_t& x) const {
     if (const auto dividend = singleton()) {
         if (const auto divisor = x.singleton()) {
             if (*divisor == 0) {
-                return interval_t(*dividend);
+                return interval_t{*dividend};
             }
-            return interval_t(*dividend % *divisor);
+            return interval_t{*dividend % *divisor};
         }
     }
     if (x.contains(0)) {
         // The divisor contains 0.
-        interval_t l(x._lb, bound_t(-1));
-        interval_t u(bound_t(1), x._ub);
+        interval_t l{x._lb, -1};
+        interval_t u{1, x._ub};
         return SRem(l) | SRem(u) | *this;
     }
     if (x.ub().is_finite() && x.lb().is_finite()) {
-        number_t min_divisor = min(abs(*x.lb().number()), abs(*x.ub().number()));
-        number_t max_divisor = max(abs(*x.lb().number()), abs(*x.ub().number()));
+        auto [xlb, xub] = x.pair_number();
+        const auto [min_divisor, max_divisor] = std::minmax({xlb.abs(), xub.abs()});
 
         if (ub() < min_divisor && -lb() < min_divisor) {
             // The modulo operation won't change the destination register.
             return *this;
         }
 
-        if (lb() < number_t{0}) {
-            if (ub() > number_t{0}) {
-                return interval_t(-(max_divisor - 1), max_divisor - 1);
+        if (lb() < 0) {
+            if (ub() > 0) {
+                return interval_t{-(max_divisor - 1), max_divisor - 1};
             } else {
-                return interval_t(-(max_divisor - 1), number_t{0});
+                return interval_t{-(max_divisor - 1), 0};
             }
         }
-        return interval_t(number_t{0}, max_divisor - 1);
+        return interval_t{0, max_divisor - 1};
     }
     // Divisor has infinite range, so result can be anything between the dividend and zero.
-    return *this | interval_t(number_t(0));
+    return *this | interval_t{0};
 }
 
 // Unsigned remainder (modulo).
@@ -187,26 +207,28 @@ interval_t interval_t::URem(const interval_t& x) const {
     if (is_bottom() || x.is_bottom()) {
         return bottom();
     }
-    if (const auto n = x.singleton()) {
-        // Divisor is a singleton:
-        //   the linear interval solver can perform many divisions where
-        //   the divisor is a singleton interval. We optimize for this case.
-        number_t c = *n;
-        if (c > number_t{0}) {
-            return interval_t(_lb.URem(bound_t{c}), _ub.URem(bound_t{c}));
+    if (const auto dividend = singleton()) {
+        if (const auto divisor = x.singleton()) {
+            if (dividend->fits_cast_to<uint64_t>() && divisor->fits_cast_to<uint64_t>()) {
+                // The BPF ISA defines modulo by 0 as resulting in the original value.
+                if (*divisor == 0) {
+                    return interval_t{*dividend};
+                }
+                uint64_t dividend_val = dividend->cast_to<uint64_t>();
+                uint64_t divisor_val = divisor->cast_to<uint64_t>();
+                return interval_t{dividend_val % divisor_val};
+            }
         }
-        // The eBPF ISA defines modulo 0 as being unchanged.
-        return *this;
     }
     if (x.contains(0)) {
         // The divisor contains 0.
-        interval_t l(x._lb, bound_t(-1));
-        interval_t u(bound_t(1), x._ub);
+        interval_t l{x._lb, -1};
+        interval_t u{1, x._ub};
         return URem(l) | URem(u) | *this;
     } else if (contains(0)) {
         // The dividend contains 0.
-        interval_t l(_lb, bound_t(-1));
-        interval_t u(bound_t(1), _ub);
+        interval_t l{_lb, -1};
+        interval_t u{1, _ub};
         return l.URem(x) | u.URem(x) | *this;
     } else {
         // Neither the dividend nor the divisor contains 0
@@ -214,13 +236,13 @@ interval_t interval_t::URem(const interval_t& x) const {
             // Divisor is infinite. A "negative" dividend could result in anything except
             // a value between the upper bound and 0, so set to top.  A "positive" dividend
             // could result in anything between 0 and the dividend - 1.
-            return (_ub < number_t{0}) ? top() : ((*this - interval_t(number_t{1})) | interval_t(number_t(0)));
-        } else if (_ub.is_finite() && (_ub.number()->cast_to<uint64_t>() < x._lb.number()->cast_to<uint64_t>())) {
+            return _ub < 0 ? top() : (*this - interval_t{1}) | interval_t{0};
+        } else if (_ub.is_finite() && _ub.number()->cast_to<uint64_t>() < x._lb.number()->cast_to<uint64_t>()) {
             // Dividend lower than divisor, so the dividend is the remainder.
             return *this;
         } else {
             number_t max_divisor{x._ub.number()->cast_to<uint64_t>()};
-            return interval_t(number_t{0}, max_divisor - 1);
+            return interval_t{0, max_divisor - 1};
         }
     }
 }
@@ -230,12 +252,12 @@ interval_t interval_t::And(const interval_t& x) const {
     if (is_bottom() || x.is_bottom()) {
         return bottom();
     }
-    assert(is_top() || (lb() >= number_t{0}));
-    assert(x.is_top() || (x.lb() >= number_t{0}));
+    assert(is_top() || (lb() >= 0));
+    assert(x.is_top() || (x.lb() >= 0));
 
     if (const auto left = singleton()) {
         if (const auto right = x.singleton()) {
-            return interval_t(*left & *right);
+            return interval_t{*left & *right};
         }
     }
     if (x == interval_t{std::numeric_limits<uint32_t>::max()}) {
@@ -252,11 +274,11 @@ interval_t interval_t::And(const interval_t& x) const {
     if (x.contains(std::numeric_limits<uint64_t>::max())) {
         return truncate_to<uint64_t>();
     } else if (!is_top() && !x.is_top()) {
-        return interval_t(number_t{0}, std::min(ub(), x.ub()));
+        return interval_t{0, std::min(ub(), x.ub())};
     } else if (!x.is_top()) {
-        return interval_t(number_t{0}, x.ub());
+        return interval_t{0, x.ub()};
     } else if (!is_top()) {
-        return interval_t(number_t{0}, ub());
+        return interval_t{0, ub()};
     } else {
         return top();
     }
@@ -268,17 +290,16 @@ interval_t interval_t::Or(const interval_t& x) const {
     }
     if (const auto left_op = singleton()) {
         if (const auto right_op = x.singleton()) {
-            return interval_t(*left_op | *right_op);
+            return interval_t{*left_op | *right_op};
         }
     }
-    if (lb() >= number_t{0} && x.lb() >= number_t{0}) {
+    if (lb() >= 0 && x.lb() >= 0) {
         if (const auto left_ub = ub().number()) {
             if (const auto right_ub = x.ub().number()) {
-                const number_t m = *left_ub > *right_ub ? *left_ub : *right_ub;
-                return interval_t(number_t{0}, m.fill_ones());
+                return interval_t{0, std::max(*left_ub, *right_ub).fill_ones()};
             }
         }
-        return interval_t(number_t{0}, bound_t::plus_infinity());
+        return interval_t{0, bound_t::plus_infinity()};
     }
     return top();
 }
@@ -289,7 +310,7 @@ interval_t interval_t::Xor(const interval_t& x) const {
     }
     if (const auto left_op = singleton()) {
         if (const auto right_op = x.singleton()) {
-            return interval_t(*left_op ^ *right_op);
+            return interval_t{*left_op ^ *right_op};
         }
     }
     return Or(x);
@@ -301,13 +322,11 @@ interval_t interval_t::Shl(const interval_t& x) const {
     }
     if (const auto shift = x.singleton()) {
         const number_t k = *shift;
-        if (k < number_t{0}) {
-            // CRAB_ERROR("lshr shift operand cannot be negative");
+        if (k < 0) {
             return top();
         }
-        // Some crazy linux drivers generate shl instructions with
-        // huge shifts.  We limit the number of times the loop is run
-        // to avoid wasting too much time on it.
+        // Some crazy linux drivers generate shl instructions with huge shifts.
+        // We limit the number of times the loop is run to avoid wasting too much time on it.
         if (k <= 128) {
             number_t factor = 1;
             for (int i = 0; k > i; i++) {
@@ -325,13 +344,11 @@ interval_t interval_t::AShr(const interval_t& x) const {
     }
     if (const auto shift = x.singleton()) {
         const number_t k = *shift;
-        if (k < number_t{0}) {
-            // CRAB_ERROR("ashr shift operand cannot be negative");
+        if (k < 0) {
             return top();
         }
-        // Some crazy linux drivers generate ashr instructions with
-        // huge shifts.  We limit the number of times the loop is run
-        // to avoid wasting too much time on it.
+        // Some crazy linux drivers generate ashr instructions with huge shifts.
+        // We limit the number of times the loop is run to avoid wasting too much time on it.
         if (k <= 128) {
             number_t factor = 1;
             for (int i = 0; k > i; i++) {
@@ -348,20 +365,9 @@ interval_t interval_t::LShr(const interval_t& x) const {
         return bottom();
     }
     if (const auto shift = x.singleton()) {
-        const number_t k = *shift;
-        if (k < number_t{0}) {
-            // CRAB_ERROR("lshr shift operand cannot be negative");
-            return top();
-        }
-        // Some crazy linux drivers generate lshr instructions with
-        // huge shifts.  We limit the number of times the loop is run
-        // to avoid wasting too much time on it.
-        if (k <= 128) {
-            if (lb() >= number_t{0} && ub().is_finite() && shift) {
-                const number_t lb = *this->lb().number();
-                const number_t ub = *this->ub().number();
-                return interval_t(lb >> k, ub >> k);
-            }
+        if (*shift > 0 && lb() >= 0 && ub().is_finite()) {
+            const auto [lb, ub] = this->pair_number();
+            return interval_t{lb >> *shift, ub >> *shift};
         }
     }
     return top();
diff --git a/src/crab/interval.hpp b/src/crab/interval.hpp
index 657b8f694..c7de6dda6 100644
--- a/src/crab/interval.hpp
+++ b/src/crab/interval.hpp
@@ -36,12 +36,6 @@ class interval_t final {
   private:
     interval_t() : _lb(number_t{0}), _ub(-1) {}
 
-    static number_t abs(const number_t& x) { return x < 0 ? -x : x; }
-
-    static number_t max(const number_t& x, const number_t& y) { return x.operator<=(y) ? y : x; }
-
-    static number_t min(const number_t& x, const number_t& y) { return x.operator<(y) ? x : y; }
-
   public:
     interval_t(const bound_t& lb, const bound_t& ub)
         : _lb(lb > ub ? bound_t{number_t{0}} : lb), _ub(lb > ub ? bound_t{-1} : ub) {}
@@ -116,12 +110,12 @@ class interval_t final {
 
     [[nodiscard]]
     bool is_bottom() const {
-        return (_lb > _ub);
+        return _lb > _ub;
     }
 
     [[nodiscard]]
     bool is_top() const {
-        return (_lb.is_infinite() && _ub.is_infinite());
+        return _lb.is_infinite() && _ub.is_infinite();
     }
 
     bool operator==(const interval_t& x) const {
@@ -225,17 +219,7 @@ class interval_t final {
 
     interval_t& operator-=(const interval_t& x) { return operator=(operator-(x)); }
 
-    interval_t operator*(const interval_t& x) const {
-        if (is_bottom() || x.is_bottom()) {
-            return bottom();
-        } else {
-            bound_t ll = _lb * x._lb;
-            bound_t lu = _lb * x._ub;
-            bound_t ul = _ub * x._lb;
-            bound_t uu = _ub * x._ub;
-            return interval_t{std::min({ll, lu, ul, uu}), std::max({ll, lu, ul, uu})};
-        }
-    }
+    interval_t operator*(const interval_t& x) const;
 
     interval_t& operator*=(const interval_t& x) { return operator=(operator*(x)); }
 
@@ -268,8 +252,8 @@ class interval_t final {
         if (is_bottom()) {
             return false;
         }
-        bound_t b(n);
-        return (_lb <= b) && (b <= _ub);
+        const bound_t b{n};
+        return _lb <= b && b <= _ub;
     }
 
     friend std::ostream& operator<<(std::ostream& o, const interval_t& interval) {
diff --git a/src/crab_utils/num_big.hpp b/src/crab_utils/num_big.hpp
index 538921a2b..3a0d57168 100644
--- a/src/crab_utils/num_big.hpp
+++ b/src/crab_utils/num_big.hpp
@@ -223,6 +223,8 @@ class number_t final {
 
     bool operator>=(const number_t& x) const { return _n >= x._n; }
 
+    number_t abs() const { return _n < 0 ? -_n : _n; }
+
     number_t operator&(const number_t& x) const { return number_t(_n & x._n); }
 
     number_t operator|(const number_t& x) const { return number_t(_n | x._n); }
diff --git a/test-data/muldiv.yaml b/test-data/muldiv.yaml
index f32754b43..f294d667d 100644
--- a/test-data/muldiv.yaml
+++ b/test-data/muldiv.yaml
@@ -56,5 +56,5 @@ code:
 
 post:
   - r0.type=number
-  - r0.svalue=-100532224
-  - r0.uvalue=18446744073609019392
+  - r0.svalue=97779712
+  - r0.uvalue=97779712
diff --git a/test-data/sdivmod.yaml b/test-data/sdivmod.yaml
index 8da36fb5b..3e59e49cd 100644
--- a/test-data/sdivmod.yaml
+++ b/test-data/sdivmod.yaml
@@ -569,3 +569,23 @@ post:
   - r2.type=number
   - r2.svalue=[7, 10]
   - r2.uvalue=[7, 10]
+
+---
+test-case: Inverted range after modulo
+
+# Ranges are and divisor are selected so that the result of modulo has
+# inverted range compared (i.e 0x8000 % 0x808 = 0x788 and 0x7FFFFFFFFFFFFFFF % 0x808 = 0x787)
+# This previously caused an issue in the verifier as the range would be bottom (instead of a valid range).
+# Constants were picked by the libfuzzer tool.
+pre: ["r0.type=number",
+      "r0.svalue=[32768, 9223372036854775807]"]
+
+code:
+  <start>: |
+    r0 s%= 0x808
+
+post:
+  - r0.type=number
+  - r0.svalue=[0, 2055]
+  - r0.uvalue=[0, 2055]
+  - r0.svalue=r0.uvalue
diff --git a/test-data/udivmod.yaml b/test-data/udivmod.yaml
index 5255c026b..85c833b86 100644
--- a/test-data/udivmod.yaml
+++ b/test-data/udivmod.yaml
@@ -573,3 +573,23 @@ post:
   - r2.type=number
   - r2.svalue=[7, 10]
   - r2.uvalue=[7, 10]
+
+---
+test-case: Inverted range after modulo
+
+# Ranges are and divisor are selected so that the result of modulo has
+# inverted range compared (i.e 0x8000 % 0x808 = 0x788 and 0x7FFFFFFFFFFFFFFF % 0x808 = 0x787)
+# This previously caused an issue in the verifier as the range would be bottom (instead of a valid range).
+# Constants were picked by the libfuzzer tool.
+pre: ["r0.type=number",
+      "r0.uvalue=[32768, 9223372036854775807]"]
+
+code:
+  <start>: |
+    r0 %= 0x808
+
+post:
+  - r0.type=number
+  - r0.svalue=[0, 2055]
+  - r0.uvalue=[0, 2055]
+  - r0.svalue=r0.uvalue