Repository: apytypes/apytypes (language: python)

Issue:
Deviations from Berkeley TestFloat
Some deviations and bugs in floating-point arithmetic were found when testing against Berkeley TestFloat:
* [x] (#408) As already known, the sign of zero is incorrect as it doesn't consider the quantization mode. 
* [x] (#407) The result sometimes becomes infinity instead of the maximum normal number for certain quantization modes. This makes sense when reading the definitions again, e.g. for TO_POS (from the 754 standard): "the result shall be the format’s floating-point number (possibly +∞) closest to and no less than the infinitely precise result". The result should therefore never become negative infinity for TO_POS.
* [x] (#428) f16_mul sometimes fails for directed quantization modes, e.g. the following for TO_POS:
```Python
APyFloat(1, 1, 1023, 5, 10) * APyFloat(1, 26, 80, 5, 10)
# Becomes APyFloat(sign=0, exp=13, man=79, exp_bits=5, man_bits=10)
# and not APyFloat(sign=0, exp=13, man=80, exp_bits=5, man_bits=10)
```
* [x] (#411) f32_mul sometimes results in a subnormal number instead of zero, however this isn't seen for f64_mul or f16_mul):
```Python
APyFloat(0, 53, 32895, 8, 23) * APyFloat(1, 1, 1999049, 8, 23)
# Becomes APyFloat(sign=1, exp=0, man=10184, exp_bits=8, man_bits=23)
# and not APyFloat(sign=1, exp=0, man=0, exp_bits=8, man_bits=23)
```
* [x] Added bullet point. The mantissa is sometimes wrong when casting from a larger format to a smaller, similar to the previous bullet. See comment below.

These were found from running
```cmd
python run_berkeley_cases.py -op all -qm all -s 1 -l 1
```
but dedicated test cases should of course be made while fixing these.


Discussion / hints:

Fix patch:
diff --git a/src/apyfloat.cc b/src/apyfloat.cc
index b142b35bc..5ef5eed13 100644
--- a/src/apyfloat.cc
+++ b/src/apyfloat.cc
@@ -256,21 +256,22 @@ void APyFloat::cast_mantissa(std::uint8_t new_man_bits, QuantizationMode quantiz
     // Check if only zeros should be added
     if (man_bits_delta <= 0) {
         if (exp >= max_exp) {
-            exp = max_exp;
-            man = 0;
+            if (do_infinity(quantization, sign)) {
+                exp = max_exp;
+                man = 0;
+            } else {
+                exp = max_exp - 1;
+                man = leading_one() - 1;
+            }
         } else {
             man <<= -man_bits_delta;
         }
         return;
     }
 
-    quantize_mantissa(man, exp, man_bits_delta, sign, leading_one(), quantization);
-
-    // Check for overflow
-    if (exp >= max_exp) {
-        exp = max_exp;
-        man = 0;
-    }
+    quantize_mantissa(
+        man, exp, max_exp, man_bits_delta, sign, leading_one(), quantization
+    );
 }
 
 // Simplified version of cast_mantissa when it is known that new_man_bits is shorter
@@ -282,13 +283,9 @@ void APyFloat::cast_mantissa_shorter(
     auto man_bits_delta = man_bits - new_man_bits;
     man_bits = new_man_bits;
     const auto max_exp = max_exponent();
-    quantize_mantissa(man, exp, man_bits_delta, sign, leading_one(), quantization);
-
-    // Check for overflow
-    if (exp >= max_exp) {
-        exp = max_exp;
-        man = 0;
-    }
+    quantize_mantissa(
+        man, exp, max_exp, man_bits_delta, sign, leading_one(), quantization
+    );
 }
 
 // Simplified version of cast_mantissa with exp = 0
@@ -305,7 +302,10 @@ void APyFloat::cast_mantissa_subnormal(
         return;
     }
 
-    quantize_mantissa(man, exp, man_bits_delta, sign, leading_one(), quantization);
+    // The overflow check with max_exponent() is not needed here, but send it anyway
+    quantize_mantissa(
+        man, exp, max_exponent(), man_bits_delta, sign, leading_one(), quantization
+    );
 }
 
 APyFloat APyFloat::_cast_to_double() const
@@ -464,7 +464,13 @@ APyFloat APyFloat::cast_from_double(
     }
 
     if (new_exp >= res.max_exponent()) {
-        return res.construct_inf();
+        if (do_infinity(get_float_quantization_mode(), res.sign)) {
+            new_exp = res.max_exponent();
+            new_man = 0;
+        } else {
+            new_exp = res.max_exponent() - 1;
+            new_man = res.man_mask();
+        }
     }
 
     res.man = new_man;
@@ -858,7 +864,14 @@ APyFloat APyFloat::operator+(const APyFloat& rhs) const
 
     // Check for overflow
     if (new_exp >= res.max_exponent()) {
-        return res.construct_inf();
+        if (do_infinity(quantization, res.sign)) {
+            res.exp = res.max_exponent();
+            res.man = 0;
+        } else {
+            res.exp = res.max_exponent() - 1;
+            res.man = res.man_mask();
+        }
+        return res;
     }
 
     // Remove leading one
@@ -982,8 +995,15 @@ APyFloat& APyFloat::operator+=(const APyFloat& rhs)
     }
 
     // Check for overflow
-    if (exp >= max_exponent()) {
-        set_to_inf();
+    const exp_t max_exp = max_exponent();
+    if (exp >= max_exp) {
+        if (do_infinity(quantization, sign)) {
+            exp = max_exp;
+            man = 0;
+        } else {
+            exp = max_exp - 1;
+            man = leading_one() - 1;
+        }
         return *this;
     }
 
@@ -1069,18 +1089,17 @@ APyFloat APyFloat::operator*(const APyFloat& y) const
         exp_t res_exp = static_cast<exp_t>(tmp_exp);
         new_man &= (two - 1);
         quantize_mantissa(
-            new_man, res_exp, man_bits_delta, sign, two_res, quantization
+            new_man,
+            res_exp,
+            res.max_exponent(),
+            man_bits_delta,
+            sign,
+            two_res,
+            quantization
         );
 
-        // Check for overflow
-        const auto max_exp = res.max_exponent();
-        if (res_exp >= max_exp) {
-            res.exp = max_exp;
-            res.man = 0;
-        } else {
-            res.man = new_man;
-            res.exp = res_exp;
-        }
+        res.man = new_man;
+        res.exp = res_exp;
         return res;
     } else {
         // Normalize both inputs
@@ -1123,7 +1142,14 @@ APyFloat APyFloat::operator*(const APyFloat& y) const
         }
 
         if (new_exp >= res.max_exponent()) {
-            return res.construct_inf();
+            if (do_infinity(quantization, res.sign)) {
+                res.exp = res.max_exponent();
+                res.man = 0;
+            } else {
+                res.exp = res.max_exponent() - 1;
+                res.man = res.man_mask();
+            }
+            return res;
         }
 
         if (apy_res.positive_greater_than_equal_pow2(0)) { // Remove leading one
@@ -1200,7 +1226,14 @@ APyFloat APyFloat::operator/(const APyFloat& y) const
 
     // Check limits
     if (new_exp >= res.max_exponent()) {
-        return res.construct_inf();
+        if (do_infinity(quantization, res.sign)) {
+            res.exp = res.max_exponent();
+            res.man = 0;
+        } else {
+            res.exp = res.max_exponent() - 1;
+            res.man = res.man_mask();
+        }
+        return res;
     }
 
     // Handle subnormal case
diff --git a/src/apyfloat_util.h b/src/apyfloat_util.h
index 8c0c6414c..a8bfac408 100644
--- a/src/apyfloat_util.h
+++ b/src/apyfloat_util.h
@@ -12,10 +12,28 @@ static constexpr std::size_t _MAN_T_SIZE_BITS = 8 * _MAN_T_SIZE_BYTES;
 static constexpr std::size_t _EXP_T_SIZE_BYTES = sizeof(exp_t);
 static constexpr std::size_t _EXP_T_SIZE_BITS = 8 * _EXP_T_SIZE_BYTES;
 
+//! Check if one should saturate to infinity or maximum normal number
+bool APY_INLINE do_infinity(QuantizationMode mode, bool sign)
+{
+    switch (mode) {
+    case QuantizationMode::TRN_ZERO:     // TO_ZERO
+    case QuantizationMode::JAM:          // JAM
+    case QuantizationMode::JAM_UNBIASED: // JAM_UNBIASED
+        return false;
+    case QuantizationMode::TRN: // TO_NEG
+        return sign;
+    case QuantizationMode::TRN_INF: // TO_POS
+        return !sign;
+    default:
+        return true;
+    }
+}
+
 //! Quantize mantissa
 void APY_INLINE quantize_mantissa(
     man_t& man,
     exp_t& exp,
+    exp_t max_exp,
     std::uint8_t bits_to_quantize,
     bool sign,
     man_t man_msb_constant,
@@ -99,18 +117,31 @@ void APY_INLINE quantize_mantissa(
             "unknown (did you pass `int` as `QuantizationMode`?)"
         );
     }
+
     man = res_man;
     man += B;
     if (man & man_msb_constant) {
         ++exp;
         man = 0;
     }
+
+    // Check for overflow. This must always be checked since other methods depend on it.
+    if (exp >= max_exp) {
+        if (do_infinity(quantization, sign)) {
+            exp = max_exp;
+            man = 0;
+        } else {
+            exp = max_exp - 1;
+            man = man_msb_constant - 1;
+        }
+    }
 }
 
 //! Quantize mantissa
 void APY_INLINE quantize_mantissa(
     man_t& man,
     exp_t& exp,
+    exp_t max_exp,
     std::uint8_t bits_to_quantize,
     bool sign,
     man_t man_msb_constant,
@@ -120,6 +151,7 @@ void APY_INLINE quantize_mantissa(
     quantize_mantissa(
         man,
         exp,
+        max_exp,
         bits_to_quantize,
         sign,
         man_msb_constant,
diff --git a/src/apyfloatarray.cc b/src/apyfloatarray.cc
index e96e37d0a..bda605b6e 100644
--- a/src/apyfloatarray.cc
+++ b/src/apyfloatarray.cc
@@ -218,14 +218,17 @@ APyFloatArray APyFloatArray::operator+(const APyFloatArray& rhs) const
             new_man &= man_mask;
 
             quantize_mantissa(
-                new_man, new_exp, 4, x.sign, final_res_leading_one, 3, 7, quantization
+                new_man,
+                new_exp,
+                res_max_exponent,
+                4,
+                x.sign,
+                final_res_leading_one,
+                3,
+                7,
+                quantization
             );
 
-            // Check for overflow
-            if (new_exp >= res_max_exponent) {
-                new_exp = res_max_exponent;
-                new_man = 0;
-            }
             res.data[i]
                 = { x.sign, static_cast<exp_t>(new_exp), static_cast<man_t>(new_man) };
         }
@@ -405,14 +408,17 @@ APyFloatArray APyFloatArray::operator+(const APyFloat& rhs) const
             new_man &= man_mask;
 
             quantize_mantissa(
-                new_man, new_exp, 4, res_sign, final_res_leading_one, 3, 7, quantization
+                new_man,
+                new_exp,
+                res_max_exponent,
+                4,
+                res_sign,
+                final_res_leading_one,
+                3,
+                7,
+                quantization
             );
 
-            // Check for overflow
-            if (new_exp >= res_max_exponent) {
-                new_exp = res_max_exponent;
-                new_man = 0;
-            }
             res.data[i] = { res_sign,
                             static_cast<exp_t>(new_exp),
                             static_cast<man_t>(new_man) };
@@ -615,6 +621,7 @@ void APyFloatArray::hadamard_multiplication(
             quantize_mantissa(
                 new_man,
                 new_exp,
+                res_max_exponent,
                 man_bits_delta,
                 res_sign,
                 two_res,
@@ -623,11 +630,6 @@ void APyFloatArray::hadamard_multiplication(
                 quantization
             );
 
-            // Check for overflow
-            if (new_exp >= res_max_exponent) {
-                new_exp = res_max_exponent;
-                new_man = 0;
-            }
             res.data[i] = { res_sign,
                             static_cast<exp_t>(new_exp),
                             static_cast<man_t>(new_man) };
@@ -811,6 +813,7 @@ APyFloatArray APyFloatArray::operator*(const APyFloat& rhs) const
             quantize_mantissa(
                 new_man,
                 new_exp,
+                res_max_exponent,
                 man_bits_delta,
                 res_sign,
                 two_res,
@@ -819,11 +822,6 @@ APyFloatArray APyFloatArray::operator*(const APyFloat& rhs) const
                 quantization
             );
 
-            // Check for overflow
-            if (new_exp >= res_max_exponent) {
-                new_exp = res_max_exponent;
-                new_man = 0;
-            }
             res.data[i] = { res_sign,
                             static_cast<exp_t>(new_exp),
                             static_cast<man_t>(new_man) };
@@ -1457,14 +1455,19 @@ APyFloat APyFloatArray::vector_sum(const QuantizationMode quantization) const
             sum_man &= man_mask;
 
             quantize_mantissa(
-                sum_man, sum_exp, 4, sum_sign, final_res_leading_one, 3, 7, quantization
+                sum_man,
+                sum_exp,
+                res_max_exponent,
+                4,
+                sum_sign,
+                final_res_leading_one,
+                3,
+                7,
+                quantization
             );
 
             // Check for overflow
-            if (sum_exp >= res_max_exponent) {
-                // Inf
-                sum_exp = res_max_exponent;
-                sum_man = 0;
+            if (sum_exp == res_max_exponent) {
                 sum_is_max_exponent = true;
             }
         }


Test patch:
diff --git a/lib/test/apyfloat/run_berkeley_cases.py b/lib/test/apyfloat/run_berkeley_cases.py
index 5e6577a2d..d80bdaa4d 100644
--- a/lib/test/apyfloat/run_berkeley_cases.py
+++ b/lib/test/apyfloat/run_berkeley_cases.py
@@ -115,22 +115,23 @@ def run_berkeley_test(
 
     if "to" in operation:
         to_exp_bits, to_man_bits = parse_format(args[2])
-        with open(output_file, "a") as f:
-            for test in test_cases:
-                tests_total += 1
-                lhs = APyFloat.from_bits(test[0], exp_bits, man_bits)
-                ref = APyFloat.from_bits(test[1], to_exp_bits, to_man_bits)
-                try:
-                    res = func_under_test(lhs, to_exp_bits, to_man_bits)
-                except Exception as e:
-                    f.write(f"lhs: {lhs!r}, ref: {ref!r}, res: {res!r}\n")
-                    f.write(f"Exception: {e}")
-                else:
-                    if (ref.is_nan and not res.is_nan) or (
-                        not ref.is_nan and not res.is_identical(ref)
-                    ):
+        with APyFloatQuantizationContext(quantization):
+            with open(output_file, "a") as f:
+                for test in test_cases:
+                    tests_total += 1
+                    lhs = APyFloat.from_bits(test[0], exp_bits, man_bits)
+                    ref = APyFloat.from_bits(test[1], to_exp_bits, to_man_bits)
+                    try:
+                        res = func_under_test(lhs, to_exp_bits, to_man_bits)
+                    except Exception as e:
                         f.write(f"lhs: {lhs!r}, ref: {ref!r}, res: {res!r}\n")
-                        tests_failed += 1
+                        f.write(f"Exception: {e}")
+                    else:
+                        if (ref.is_nan and not res.is_nan) or (
+                            not ref.is_nan and not res.is_identical(ref)
+                        ):
+                            f.write(f"lhs: {lhs!r}, ref: {ref!r}, res: {res!r}\n")
+                            tests_failed += 1
     else:
         with APyFloatQuantizationContext(quantization):
             with open(output_file, "a") as f:
@@ -290,7 +291,7 @@ def print_summary(summary, output_file: str, seed: int, level: int) -> None:
 
     success_rate = (1 - float(failed_tests) / total_tests) * 100
     print(
-        f"Ran {total_tests} tests in total whereof {failed_tests} failed, giving a {success_rate:.2f}% overall success rate."
+        f"Ran {total_tests} tests in total whereof {failed_tests} failed, giving an overall success rate of {success_rate:.2f}%."
     )
 
 
diff --git a/lib/test/apyfloat/test_rounding.py b/lib/test/apyfloat/test_rounding.py
index 876951fe0..2c06efcf4 100644
--- a/lib/test/apyfloat/test_rounding.py
+++ b/lib/test/apyfloat/test_rounding.py
@@ -40,6 +40,14 @@ def test_to_pos(self):
             res = APyFloat(1, 15, 2, 5, 2) + APyFloat(0, 1, 0, 5, 2)
             assert res.is_identical(APyFloat(1, 15, 1, 5, 2))
 
+            # Two big positive numbers should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) + APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 31, 0, 5, 2))
+
+            # Two big negative numbers should saturate
+            res = APyFloat(1, 30, 3, 5, 2) + APyFloat(1, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 30, 3, 5, 2))
+
     def test_to_neg(self):
         with APyFloatQuantizationContext(QuantizationMode.TO_NEG):
             # 1.5 + very small number should quantize to 1.5
@@ -50,6 +58,14 @@ def test_to_neg(self):
             res = APyFloat(1, 15, 2, 5, 2) + APyFloat(0, 1, 0, 5, 2)
             assert res.is_identical(APyFloat(1, 15, 2, 5, 2))
 
+            # Two big positive numbers should saturate
+            res = APyFloat(0, 30, 3, 5, 2) + APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 30, 3, 5, 2))
+
+            # Two big negative numbers should become infinity
+            res = APyFloat(1, 30, 3, 5, 2) + APyFloat(1, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 31, 0, 5, 2))
+
     def test_to_zero(self):
         with APyFloatQuantizationContext(QuantizationMode.TO_ZERO):
             # 1.5 + relatively big number (0.21875) should still quantize to 1.5
@@ -64,6 +80,14 @@ def test_to_zero(self):
             )
             assert res.is_identical(APyFloat(1, 15, 2, 5, 2))
 
+            # Two big positive numbers should saturate
+            res = APyFloat(0, 30, 3, 5, 2) + APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 30, 3, 5, 2))
+
+            # Two big negative numbers should saturate
+            res = APyFloat(1, 30, 3, 5, 2) + APyFloat(1, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 30, 3, 5, 2))
+
     def test_to_ties_even(self):
         with APyFloatQuantizationContext(QuantizationMode.TIES_EVEN):
             # 1.5 + relatively big number (0.21875) should quantize to 1.75
@@ -102,6 +126,14 @@ def test_to_ties_even(self):
             res = APyFloat(1, 15, 3, 5, 2) + APyFloat(1, 12, 0, 5, 2)
             assert res.is_identical(APyFloat(1, 16, 0, 5, 2))
 
+            # Two big positive numbers should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) + APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 31, 0, 5, 2))
+
+            # Two big negative numbers should become infinity
+            res = APyFloat(1, 30, 3, 5, 2) + APyFloat(1, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 31, 0, 5, 2))
+
     def test_to_ties_away(self):
         with APyFloatQuantizationContext(QuantizationMode.TIES_EVEN):
             # 1.5 + relatively big number (0.21875) should quantize to 1.75
@@ -140,6 +172,14 @@ def test_to_ties_away(self):
             res = APyFloat(1, 15, 3, 5, 2) + APyFloat(1, 12, 0, 5, 2)
             assert res.is_identical(APyFloat(1, 16, 0, 5, 2))
 
+            # Two big positive numbers should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) + APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 31, 0, 5, 2))
+
+            # Two big negative numbers should become infinity
+            res = APyFloat(1, 30, 3, 5, 2) + APyFloat(1, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 31, 0, 5, 2))
+
     def test_jam(self):
         with APyFloatQuantizationContext(QuantizationMode.JAM):
             # 1.5 + very small number should quantize to 1.75
@@ -154,6 +194,14 @@ def test_jam(self):
             res = APyFloat(0, 15, 0, 5, 2) + APyFloat(0, 15, 0, 5, 2)
             assert res.is_identical(APyFloat(0, 16, 1, 5, 2))
 
+            # Two big positive numbers should saturate
+            res = APyFloat(0, 30, 3, 5, 2) + APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 30, 3, 5, 2))
+
+            # Two big negative numbers should saturate
+            res = APyFloat(1, 30, 3, 5, 2) + APyFloat(1, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 30, 3, 5, 2))
+
     def test_stoch_weighted(self):
         with APyFloatQuantizationContext(QuantizationMode.STOCH_WEIGHTED):
             # 1.5 + 0.125 should quantize to 1.5 or 1.75
@@ -216,6 +264,14 @@ def test_to_pos(self):
             res = APyFloat(1, 0, 1, 4, 3) * APyFloat(0, 0, 1, 4, 3)
             assert res.is_identical(APyFloat(1, 0, 0, 4, 3))
 
+            # Big positive number should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) * APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 31, 0, 5, 2))
+
+            # Big negative number should saturate
+            res = APyFloat(1, 30, 3, 5, 2) * APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 30, 3, 5, 2))
+
     def test_to_neg(self):
         with APyFloatQuantizationContext(QuantizationMode.TO_NEG):
             # Should round down
@@ -241,6 +297,14 @@ def test_to_neg(self):
             res = APyFloat(1, 0, 1, 4, 3) * APyFloat(0, 0, 1, 4, 3)
             assert res.is_identical(APyFloat(1, 0, 1, 4, 3))
 
+            # Big positive number should saturate
+            res = APyFloat(0, 30, 3, 5, 2) * APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 30, 3, 5, 2))
+
+            # Big negative number should become infinity
+            res = APyFloat(1, 30, 3, 5, 2) * APyFloat(0, 30, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 31, 0, 5, 2))
+
 
 @pytest.mark.float_div
 class TestAPyFloatQuantizationDiv:
@@ -258,6 +322,14 @@ def test_to_pos(self):
             res = APyFloat(1, 15, 2, 5, 2) / APyFloat(0, 15, 1, 5, 2)
             assert res.is_identical(APyFloat(1, 15, 0, 5, 2))
 
+            # Should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(0, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 31, 0, 5, 2))
+
+            # Should saturate
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(1, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 30, 3, 5, 2))
+
     def test_to_neg(self):
         with APyFloatQuantizationContext(QuantizationMode.TO_NEG):
             # 1.5 / 1.25 (=1.2) should quantize to 1.0
@@ -268,6 +340,14 @@ def test_to_neg(self):
             res = APyFloat(1, 15, 2, 5, 2) / APyFloat(0, 15, 1, 5, 2)
             assert res.is_identical(APyFloat(1, 15, 1, 5, 2))
 
+            # Should saturate
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(0, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 30, 3, 5, 2))
+
+            # Should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(1, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 31, 0, 5, 2))
+
     def test_to_zero(self):
         with APyFloatQuantizationContext(QuantizationMode.TO_ZERO):
             # 1.25 / 1.5 should quantize to 0.75
@@ -278,6 +358,14 @@ def test_to_zero(self):
             res = APyFloat(1, 15, 1, 5, 2) / APyFloat(0, 15, 2, 5, 2)
             assert res.is_identical(APyFloat(1, 14, 2, 5, 2))
 
+            # Should saturate
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(0, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 30, 3, 5, 2))
+
+            # Should saturate
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(1, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 30, 3, 5, 2))
+
     def test_to_ties_even(self):
         with APyFloatQuantizationContext(QuantizationMode.TIES_EVEN):
             # 1.25 / 1.5 should quantize to closes which is 0.875
@@ -312,6 +400,14 @@ def test_to_ties_even(self):
             res = APyFloat(1, 0, 3, 5, 2) / APyFloat(0, 16, 0, 5, 2)
             assert res.is_identical(APyFloat(1, 0, 2, 5, 2))
 
+            # Should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(0, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 31, 0, 5, 2))
+
+            # Should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(1, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 31, 0, 5, 2))
+
     def test_to_ties_away(self):
         with APyFloatQuantizationContext(QuantizationMode.TIES_AWAY):
             # 1.25 / 1.5 should quantize to closes which is 0.875
@@ -346,6 +442,14 @@ def test_to_ties_away(self):
             res = APyFloat(1, 0, 3, 5, 2) / APyFloat(0, 16, 0, 5, 2)
             assert res.is_identical(APyFloat(1, 0, 2, 5, 2))
 
+            # Should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(0, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 31, 0, 5, 2))
+
+            # Should become infinity
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(1, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 31, 0, 5, 2))
+
     def test_jam(self):
         with APyFloatQuantizationContext(QuantizationMode.JAM):
             # 1.5 / 1.25 (=1.2) should quantize to 1.25
@@ -360,6 +464,14 @@ def test_jam(self):
             res = APyFloat(0, 17, 0, 5, 2) / APyFloat(0, 16, 0, 5, 2)
             assert res.is_identical(APyFloat(0, 16, 1, 5, 2))
 
+            # Should saturate
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(0, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(0, 30, 3, 5, 2))
+
+            # Should saturate
+            res = APyFloat(0, 30, 3, 5, 2) / APyFloat(1, 0, 3, 5, 2)
+            assert res.is_identical(APyFloat(1, 30, 3, 5, 2))
+
     def test_stoch_weighted(self):
         with APyFloatQuantizationContext(QuantizationMode.STOCH_WEIGHTED):
             # 1.5 / 1.25 (=1.2) should quantize to 1.0 or 1.25
@@ -519,6 +631,15 @@ def test_quantization_to_pos(self):
             .is_identical(APyFloat(1, 5, 0b101, 5, 3))
         )  # Round down
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 15, 0, 4, 3))
+        )  # Should become infinity
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 14, 7, 4, 3))
+        )  # Should saturate
+
     def test_quantization_to_neg(self):
         apytypes.set_float_quantization_mode(QuantizationMode.TO_NEG)
         # Quantization from 0.xx
@@ -609,6 +730,15 @@ def test_quantization_to_neg(self):
             .is_identical(APyFloat(1, 5, 0b110, 5, 3))
         )  # Round up
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 14, 7, 4, 3))
+        )  # Should saturate
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 15, 0, 4, 3))
+        )  # Should become infinity
+
     def test_quantization_to_away(self):
         apytypes.set_float_quantization_mode(QuantizationMode.TO_AWAY)
         # Quantization from 0.xx
@@ -699,6 +829,15 @@ def test_quantization_to_away(self):
             .is_identical(APyFloat(1, 5, 0b110, 5, 3))
         )  # Round down
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 15, 0, 4, 3))
+        )  # Should become infinity
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 15, 0, 4, 3))
+        )  # Should become infinity
+
     def test_quantization_to_zero(self):
         apytypes.set_float_quantization_mode(QuantizationMode.TO_ZERO)
         # Quantization from 0.xx
@@ -789,6 +928,15 @@ def test_quantization_to_zero(self):
             .is_identical(APyFloat(1, 5, 0b101, 5, 3))
         )  # Round down
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 14, 7, 4, 3))
+        )  # Should saturate
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 14, 7, 4, 3))
+        )  # Should saturate
+
     def test_quantization_ties_even(self):
         apytypes.set_float_quantization_mode(QuantizationMode.TIES_EVEN)
         # Quantization from 0.xx
@@ -879,6 +1027,15 @@ def test_quantization_ties_even(self):
             .is_identical(APyFloat(1, 5, 0b110, 5, 3))
         )  # Round up
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 15, 0, 4, 3))
+        )  # Should become infinity
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 15, 0, 4, 3))
+        )  # Should become infinity
+
     def test_quantization_ties_odd(self):
         apytypes.set_float_quantization_mode(QuantizationMode.TIES_ODD)
         # Quantization from 0.xx
@@ -969,6 +1126,15 @@ def test_quantization_ties_odd(self):
             .is_identical(APyFloat(1, 5, 0b110, 5, 3))
         )  # Round up
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 15, 0, 4, 3))
+        )  # Should become infinity
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 15, 0, 4, 3))
+        )  # Should become infinity
+
     def test_quantization_ties_pos(self):
         apytypes.set_float_quantization_mode(QuantizationMode.TIES_POS)
         # Quantization from 0.xx
@@ -1059,6 +1225,15 @@ def test_quantization_ties_pos(self):
             .is_identical(APyFloat(1, 5, 0b110, 5, 3))
         )  # Round up
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 15, 0, 4, 3))
+        )  # Should become infinity
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 15, 0, 4, 3))
+        )  # Should become infinity
+
     def test_quantization_ties_neg(self):
         apytypes.set_float_quantization_mode(QuantizationMode.TIES_NEG)
         # Quantization from 0.xx
@@ -1149,6 +1324,15 @@ def test_quantization_ties_neg(self):
             .is_identical(APyFloat(1, 5, 0b110, 5, 3))
         )  # Round up
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 15, 0, 4, 3))
+        )  # Should become infinity
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 15, 0, 4, 3))
+        )  # Should become infinity
+
     def test_quantization_ties_away(self):
         apytypes.set_float_quantization_mode(QuantizationMode.TIES_AWAY)
         # Quantization from 0.xx
@@ -1239,6 +1423,15 @@ def test_quantization_ties_away(self):
             .is_identical(APyFloat(1, 5, 0b110, 5, 3))
         )  # Round up
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 15, 0, 4, 3))
+        )  # Should become infinity
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 15, 0, 4, 3))
+        )  # Should become infinity
+
     def test_quantization_ties_to_zero(self):
         apytypes.set_float_quantization_mode(QuantizationMode.TIES_ZERO)
         # Quantization from 0.xx
@@ -1329,6 +1522,15 @@ def test_quantization_ties_to_zero(self):
             .is_identical(APyFloat(1, 5, 0b110, 5, 3))
         )  # Round up
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 15, 0, 4, 3))
+        )  # Should become infinity
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 15, 0, 4, 3))
+        )  # Should become infinity
+
     def test_quantization_jamming(self):
         apytypes.set_float_quantization_mode(QuantizationMode.JAM)
         # Quantization from 0.xx
@@ -1419,6 +1621,15 @@ def test_quantization_jamming(self):
             .is_identical(APyFloat(1, 5, 0b101, 5, 3))
         )
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 14, 7, 4, 3))
+        )  # Should saturate
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 14, 7, 4, 3))
+        )  # Should saturate
+
     def test_quantization_unbiased_jamming(self):
         apytypes.set_float_quantization_mode(QuantizationMode.JAM_UNBIASED)
         # Quantization from 0.xx
@@ -1445,6 +1656,15 @@ def test_quantization_unbiased_jamming(self):
             .is_identical(APyFloat(0, 5, 0b101, 5, 3))
         )
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 14, 7, 4, 3))
+        )  # Should saturate
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 14, 7, 4, 3))
+        )  # Should saturate
+
     def test_quantization_magnitude_truncation(self):
         # Shouldn't be used for floating-point
         apytypes.set_float_quantization_mode(QuantizationMode.TRN_MAG)
@@ -1536,6 +1756,15 @@ def test_quantization_magnitude_truncation(self):
             .is_identical(APyFloat(1, 5, 0b110, 5, 3))
         )
 
+        # Test handling of infinities
+        assert (
+            APyFloat(0, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(0, 15, 0, 4, 3))
+        )  # Should become infinity
+
+        assert (
+            APyFloat(1, 30, 0, 5, 5).cast(4, 3).is_identical(APyFloat(1, 15, 0, 4, 3))
+        )  # Should become infinity
+
     def test_quantization_stochastic_weighted(self):
         """
         A bit naive, but test that a value can be rounded both up and down in 1000 tries.
diff --git a/lib/test/apyfloatarray/test_arithmetic.py b/lib/test/apyfloatarray/test_arithmetic.py
index 3ad58badc..0be4c8803 100644
--- a/lib/test/apyfloatarray/test_arithmetic.py
+++ b/lib/test/apyfloatarray/test_arithmetic.py
@@ -1,5 +1,10 @@
 import pytest
-from apytypes import APyFloat, APyFloatArray
+from apytypes import (
+    APyFloat,
+    APyFloatArray,
+    APyFloatQuantizationContext,
+    QuantizationMode,
+)
 
 
 @pytest.mark.float_array
@@ -804,3 +809,139 @@ def test_scalar_array_mul_with_one():
             res = a * b
             assert res.is_identical(b)
             assert (b * a).is_identical(b)
+
+
+def test_array_infinity_saturation_cast():
+    # Big positive number should become infinity
+    res = APyFloatArray([0], [30], [2], 5, 2).cast(4, 2)
+    assert res.is_identical(APyFloatArray([0], [15], [0], 4, 2))
+
+    # Big negative number should become infinity
+    res = APyFloatArray([1], [30], [2], 5, 2).cast(4, 2)
+    assert res.is_identical(APyFloatArray([1], [15], [0], 4, 2))
+
+    with APyFloatQuantizationContext(QuantizationMode.TO_ZERO):
+        # Big positive number should saturate
+        res = APyFloatArray([0], [30], [2], 5, 2).cast(4, 2)
+        assert res.is_identical(APyFloatArray([0], [14], [3], 4, 2))
+
+        # Big negative number should saturate
+        res = APyFloatArray([1], [30], [2], 5, 2).cast(4, 2)
+        assert res.is_identical(APyFloatArray([1], [14], [3], 4, 2))
+
+
+@pytest.mark.float_mul
+def test_array_infinity_saturation_mul():
+    # Array x Scalar
+    # Big positive number should become infinity
+    res = APyFloatArray([0], [30], [2], 5, 2) * APyFloat(0, 30, 3, 5, 2)
+    assert res.is_identical(APyFloatArray([0], [31], [0], 5, 2))
+
+    # Big negative number should become infinity
+    res = APyFloatArray([1], [30], [2], 5, 2) * APyFloat(0, 30, 3, 5, 2)
+    assert res.is_identical(APyFloatArray([1], [31], [0], 5, 2))
+
+    with APyFloatQuantizationContext(QuantizationMode.TO_ZERO):
+        # Big positive number should saturate
+        res = APyFloatArray([0], [30], [2], 5, 2) * APyFloat(0, 30, 3, 5, 2)
+        assert res.is_identical(APyFloatArray([0], [30], [3], 5, 2))
+
+        # Big negative number should saturate
+        res = APyFloatArray([1], [30], [2], 5, 2) * APyFloat(0, 30, 3, 5, 2)
+        assert res.is_identical(APyFloatArray([1], [30], [3], 5, 2))
+
+    # Array x Array
+    # Big positive number should become infinity
+    res = APyFloatArray([0], [30], [2], 5, 2) * APyFloatArray([0], [30], [2], 5, 2)
+    assert res.is_identical(APyFloatArray([0], [31], [0], 5, 2))
+
+    # Big negative number should become infinity
+    res = APyFloatArray([1], [30], [2], 5, 2) * APyFloatArray([0], [30], [2], 5, 2)
+    assert res.is_identical(APyFloatArray([1], [31], [0], 5, 2))
+
+    with APyFloatQuantizationContext(QuantizationMode.TO_ZERO):
+        # Big positive number should saturate
+        res = APyFloatArray([0], [30], [2], 5, 2) * APyFloatArray([0], [30], [2], 5, 2)
+        assert res.is_identical(APyFloatArray([0], [30], [3], 5, 2))
+
+        # Big negative number should saturate
+        res = APyFloatArray([1], [30], [2], 5, 2) * APyFloatArray([0], [30], [2], 5, 2)
+        assert res.is_identical(APyFloatArray([1], [30], [3], 5, 2))
+
+
+@pytest.mark.float_add
+def test_array_infinity_saturation_add():
+    # Array x Scalar
+    # Big positive number should become infinity
+    res = APyFloatArray([0], [30], [2], 5, 2) + APyFloat(0, 30, 3, 5, 2)
+    assert res.is_identical(APyFloatArray([0], [31], [0], 5, 2))
+
+    # Big negative number should become infinity
+    res = APyFloatArray([1], [30], [2], 5, 2) + APyFloat(1, 30, 3, 5, 2)
+    assert res.is_identical(APyFloatArray([1], [31], [0], 5, 2))
+
+    with APyFloatQuantizationContext(QuantizationMode.TO_ZERO):
+        # Big positive number should saturate
+        res = APyFloatArray([0], [30], [2], 5, 2) + APyFloat(0, 30, 3, 5, 2)
+        assert res.is_identical(APyFloatArray([0], [30], [3], 5, 2))
+
+        # Big negative number should saturate
+        res = APyFloatArray([1], [30], [2], 5, 2) + APyFloat(1, 30, 3, 5, 2)
+        assert res.is_identical(APyFloatArray([1], [30], [3], 5, 2))
+
+    # Array x Array
+    # Big positive number should become infinity
+    res = APyFloatArray([0], [30], [2], 5, 2) + APyFloatArray([0], [30], [2], 5, 2)
+    assert res.is_identical(APyFloatArray([0], [31], [0], 5, 2))
+
+    # Big negative number should become infinity
+    res = APyFloatArray([1], [30], [2], 5, 2) + APyFloatArray([1], [30], [2], 5, 2)
+    assert res.is_identical(APyFloatArray([1], [31], [0], 5, 2))
+
+    with APyFloatQuantizationContext(QuantizationMode.TO_ZERO):
+        # Big positive number should saturate
+        res = APyFloatArray([0], [30], [2], 5, 2) + APyFloatArray([0], [30], [2], 5, 2)
+        assert res.is_identical(APyFloatArray([0], [30], [3], 5, 2))
+
+        # Big negative number should saturate
+        res = APyFloatArray([1], [30], [2], 5, 2) + APyFloatArray([1], [30], [2], 5, 2)
+        assert res.is_identical(APyFloatArray([1], [30], [3], 5, 2))
+
+
+@pytest.mark.float_div
+def test_array_infinity_saturation_div():
+    # Array x Scalar
+    # Big positive number should become infinity
+    res = APyFloatArray([0], [30], [2], 5, 2) / APyFloat(0, 1, 3, 5, 2)
+    assert res.is_identical(APyFloatArray([0], [31], [0], 5, 2))
+
+    # Big negative number should become infinity
+    res = APyFloatArray([1], [30], [2], 5, 2) / APyFloat(0, 1, 3, 5, 2)
+    assert res.is_identical(APyFloatArray([1], [31], [0], 5, 2))
+
+    with APyFloatQuantizationContext(QuantizationMode.TO_ZERO):
+        # Big positive number should saturate
+        res = APyFloatArray([0], [30], [2], 5, 2) / APyFloat(0, 1, 3, 5, 2)
+        assert res.is_identical(APyFloatArray([0], [30], [3], 5, 2))
+
+        # Big negative number should saturate
+        res = APyFloatArray([1], [30], [2], 5, 2) / APyFloat(0, 1, 3, 5, 2)
+        assert res.is_identical(APyFloatArray([1], [30], [3], 5, 2))
+
+    # Array x Array
+    # Big positive number should become infinity
+    res = APyFloatArray([0], [30], [2], 5, 2) / APyFloatArray([0], [1], [2], 5, 2)
+    assert res.is_identical(APyFloatArray([0], [31], [0], 5, 2))
+
+    # Big negative number should become infinity
+    res = APyFloatArray([1], [30], [2], 5, 2) / APyFloatArray([0], [1], [2], 5, 2)
+    assert res.is_identical(APyFloatArray([1], [31], [0], 5, 2))
+
+    with APyFloatQuantizationContext(QuantizationMode.TO_ZERO):
+        # Big positive number should saturate
+        res = APyFloatArray([0], [30], [2], 5, 2) / APyFloatArray([0], [1], [2], 5, 2)
+        assert res.is_identical(APyFloatArray([0], [30], [3], 5, 2))
+
+        # Big negative number should saturate
+        res = APyFloatArray([1], [30], [2], 5, 2) / APyFloatArray([0], [1], [2], 5, 2)
+        assert res.is_identical(APyFloatArray([1], [30], [3], 5, 2))