Repository: apytypes/apytypes (language: python)

Issue:
Bias conversion
The support for mixed-biases are not that well developed, beyond the IEEE-case. I think this is the way to convert general biases.

Two numbers with bias b1 and b2 and exponent bits e1 and e2, converted to e3 exponent bits.

```
b3 = ((b1 + 1)/(2**e1) + (b2 + 1)/(2**e2))*(2**(e3))/2 - 1
``` 
If `e3 >= max(e1, e2)`
```
b3 = ((b1 + 1) << (e3 - e1) + (b2 + 1) << (e3 - e2)) >> 1
```

This should lead to that numbers with IEEE-biases are kept the same and other combinations will get the "average" bias, which is most likely the best we can do. In this case it is rounded downwards, which I think is as good as the other option.
Premature quantization in floating-point accumulator context
The accumulator context for floating-point is implemented by resizing the operands to the accumulator format before the multiplication to avoid quantizing too early. However this only works if the accumulator format can fit the operands. Instead one should perform the multiplication and quantize directly to the accumulator format -- this may require writing a new multiplication method.


Discussion / hints:
Seems like the second expression should be this?
```Python
((((b1 + 1) << (e3 - e1)) + ((b2 + 1) << (e3 - e2))) >> 1) - 1
```

I think your suggestion makes sense but it can lead to some unexpected results. Consider adding a large number with zero as
```Python
x = APyFloat(sign=0, exp=21, man=0, exp_bits=5, man_bits=2, bias=5)
y = APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=25)
# APyFloat(sign=0, exp=31, man=0, exp_bits=5, man_bits=2)
x+y
```
The resulting bias will be $15$ so the exponent of $x$ will be increased by $10$ and thus become infinity.

I think in general one should side with the format representing numbers with greater magnitude, but I'm not sure how that would be formulated. Perhaps something like this:
```Python
b3 = min((b1 + 1) << (e3 - e1), (b2 + 1) << (e3 - e2)) - 1
```
Won't a similar thing happen for small numbers? Adding zero will turn the reault into zero?
Yes that's true.. I think both make sense but in different ways


Fix patch:
diff --git a/src/apyfloat.cc b/src/apyfloat.cc
index a1c642e70..37a8b3628 100644
--- a/src/apyfloat.cc
+++ b/src/apyfloat.cc
@@ -491,6 +491,7 @@ APyFloat APyFloat::cast_from_double(
     return res;
 }
 
+//! Cast to a larger format
 APyFloat APyFloat::cast_no_quant(
     std::uint8_t new_exp_bits, std::uint8_t new_man_bits, std::optional<exp_t> new_bias
 ) const
@@ -732,24 +733,22 @@ APyFloat APyFloat::operator+(const APyFloat& rhs) const
     } else {
         std::uint8_t res_exp_bits = std::max(exp_bits, rhs.exp_bits);
         std::uint8_t res_man_bits = std::max(man_bits, rhs.man_bits);
-        exp_t res_bias = APyFloat::ieee_bias(res_exp_bits);
-        // Cast once to resulting word length to get faster comparisons later
-        x = cast_no_quant(res_exp_bits, res_man_bits, res_bias);
-        y = rhs.cast_no_quant(res_exp_bits, res_man_bits, res_bias);
+        exp_t res_bias
+            = calc_bias(res_exp_bits, exp_bits, bias, rhs.exp_bits, rhs.bias);
         res = APyFloat(res_exp_bits, res_man_bits, res_bias);
 
         if (is_zero()) {
-            if (y.is_zero()) {
-                const bool new_sign = (sign == y.sign)
+            if (rhs.is_zero()) {
+                const bool new_sign = (sign == rhs.sign)
                     ? sign
                     : get_float_quantization_mode() == QuantizationMode::TRN;
                 res.set_to_zero(new_sign);
                 return res;
             }
-            return y;
+            return rhs.cast(res_exp_bits, res_man_bits, res_bias);
         }
         if (rhs.is_zero()) {
-            return x;
+            return cast(res_exp_bits, res_man_bits, res_bias);
         }
         // Handle the NaN and inf cases
         if (is_max_exponent() || rhs.is_max_exponent()) {
@@ -762,9 +761,13 @@ APyFloat APyFloat::operator+(const APyFloat& rhs) const
             res.set_to_inf();
             return res;
         }
+
+        x = cast_no_quant(res_exp_bits, res_man_bits, bias);
+        y = rhs.cast_no_quant(res_exp_bits, res_man_bits, rhs.bias);
+
         if (x.abs() < y.abs()) {
             std::swap(x, y);
-        } else if (sign != rhs.sign && x.exp == y.exp && x.man == y.man) {
+        } else if (sign != rhs.sign && x.true_exp() == y.true_exp() && x.man == y.man) {
             // +0 for all quantization modes except TO_NEG
             res.set_to_zero(get_float_quantization_mode() == QuantizationMode::TRN);
             return res;
@@ -1053,9 +1056,10 @@ APyFloat& APyFloat::operator+=(const APyFloat& rhs)
 
 APyFloat APyFloat::operator*(const APyFloat& y) const
 {
-    auto res_exp_bits = std::max(exp_bits, y.exp_bits);
-    auto res_man_bits = std::max(man_bits, y.man_bits);
-    APyFloat res(res_exp_bits, res_man_bits);
+    const auto res_exp_bits = std::max(exp_bits, y.exp_bits);
+    const auto res_man_bits = std::max(man_bits, y.man_bits);
+    const auto res_bias = calc_bias(res_exp_bits, exp_bits, bias, y.exp_bits, y.bias);
+    APyFloat res(res_exp_bits, res_man_bits, res_bias);
 
     // Calculate sign
     res.sign = sign ^ y.sign;
@@ -1209,7 +1213,10 @@ APyFloat APyFloat::operator*(const APyFloat& y) const
 
 APyFloat APyFloat::operator/(const APyFloat& y) const
 {
-    APyFloat res(std::max(exp_bits, y.exp_bits), std::max(man_bits, y.man_bits));
+    const auto res_exp_bits = std::max(exp_bits, y.exp_bits);
+    const auto res_man_bits = std::max(man_bits, y.man_bits);
+    const auto res_bias = calc_bias(res_exp_bits, exp_bits, bias, y.exp_bits, y.bias);
+    APyFloat res(res_exp_bits, res_man_bits, res_bias);
 
     // Calculate sign
     res.sign = sign ^ y.sign;
@@ -1753,7 +1760,7 @@ void APyFloat::set_to_nan(std::optional<bool> new_sign, man_t payload /*= 1*/)
 
 APyFloat APyFloat::construct_zero(std::optional<bool> new_sign) const
 {
-    return APyFloat(new_sign.value_or(sign), 0, 0, exp_bits, man_bits);
+    return APyFloat(new_sign.value_or(sign), 0, 0, exp_bits, man_bits, bias);
 }
 
 APyFloat APyFloat::construct_inf(std::optional<bool> new_sign) const
@@ -1765,7 +1772,7 @@ APyFloat
 APyFloat::construct_nan(std::optional<bool> new_sign, man_t payload /*= 1*/) const
 {
     return APyFloat(
-        new_sign.value_or(sign), max_exponent(), payload, exp_bits, man_bits
+        new_sign.value_or(sign), max_exponent(), payload, exp_bits, man_bits, bias
     );
 }
 
diff --git a/src/apyfloat_util.cc b/src/apyfloat_util.cc
index 94c77cc80..f3e9361be 100644
--- a/src/apyfloat_util.cc
+++ b/src/apyfloat_util.cc
@@ -1,4 +1,5 @@
 #include "apyfloat_util.h"
+#include <math.h>
 
 /* Helper functions */
 man_t ipow(man_t base, unsigned int n)
@@ -23,6 +24,16 @@ man_t ipow(man_t base, unsigned int n)
     return result;
 }
 
+exp_t calc_bias_general(
+    int new_exp_bits, int exp_bits1, exp_t bias1, int exp_bits2, exp_t bias2
+)
+{
+    return ((double)(bias1 + 1) / (1 << exp_bits1)
+            + (double)(bias2 + 1) / (1 << exp_bits2))
+        * (1 << (new_exp_bits - 1))
+        - 1;
+}
+
 void quantize_apymantissa(
     APyFixed& apyman, bool sign, int bits, QuantizationMode quantization
 )
diff --git a/src/apyfloat_util.h b/src/apyfloat_util.h
index 3e74c5356..3ee353bdf 100644
--- a/src/apyfloat_util.h
+++ b/src/apyfloat_util.h
@@ -29,6 +29,22 @@ bool APY_INLINE do_infinity(QuantizationMode mode, bool sign)
     }
 }
 
+//! Calculate new bias. Assumes new_exp_bits is larger than exp_bits1 and exp_bits2.
+exp_t APY_INLINE
+calc_bias(int new_exp_bits, int exp_bits1, exp_t bias1, int exp_bits2, exp_t bias2)
+{
+    return ((((bias1 + 1) << (new_exp_bits - exp_bits1))
+             + ((bias2 + 1) << (new_exp_bits - exp_bits2)))
+            >> 1)
+        - 1;
+}
+
+//! General calculation of new bias. Should only be used if new_exp_bits can be strictly
+//! less than exp_bits1 or exp_bits2.
+exp_t calc_bias_general(
+    int new_exp_bits, int exp_bits1, exp_t bias1, int exp_bits2, exp_t bias2
+);
+
 //! Quantize mantissa
 void APY_INLINE quantize_mantissa(
     man_t& man,
diff --git a/src/apyfloatarray.cc b/src/apyfloatarray.cc
index b96547ddd..8e77b5ae1 100644
--- a/src/apyfloatarray.cc
+++ b/src/apyfloatarray.cc
@@ -25,12 +25,7 @@ APyFloatArray::APyFloatArray(
     : exp_bits(exp_bits)
     , man_bits(man_bits)
 {
-    this->bias = APyFloat::ieee_bias(exp_bits);
-    if (bias.value_or(this->bias) != this->bias) {
-        throw NotImplementedException(
-            "Not implemented: APyFloatArray with non IEEE-like bias"
-        );
-    }
+    this->bias = bias.value_or(APyFloat::ieee_bias(exp_bits));
 
     const auto signs_shape = python_sequence_extract_shape(sign_seq);
     const auto exps_shape = python_sequence_extract_shape(exp_seq);
@@ -99,7 +94,7 @@ APyFloatArray APyFloatArray::operator+(const APyFloatArray& rhs) const
     if (same_type_as(rhs) && (max_man_bits <= _MAN_T_SIZE_BITS)
         && (quantization != QuantizationMode::STOCH_WEIGHTED)) {
         // Result array
-        APyFloatArray res(shape, exp_bits, man_bits);
+        APyFloatArray res(shape, exp_bits, man_bits, bias);
 
         APyFloatData x, y;
         const exp_t res_max_exponent = ((1ULL << exp_bits) - 1);
@@ -240,9 +235,11 @@ APyFloatArray APyFloatArray::operator+(const APyFloatArray& rhs) const
         return res;
     }
     // Calculate new format
-    APyFloatArray res(
-        shape, std::max(exp_bits, rhs.exp_bits), std::max(man_bits, rhs.man_bits)
-    );
+    const auto res_exp_bits = std::max(exp_bits, rhs.exp_bits);
+    const auto res_man_bits = std::max(man_bits, rhs.man_bits);
+    const auto res_bias
+        = calc_bias(res_exp_bits, exp_bits, bias, rhs.exp_bits, rhs.bias);
+    APyFloatArray res(shape, res_exp_bits, res_man_bits, res_bias);
 
     APyFloat lhs_scalar(exp_bits, man_bits, bias);
     APyFloat rhs_scalar(rhs.exp_bits, rhs.man_bits, rhs.bias);
@@ -270,7 +267,7 @@ APyFloatArray APyFloatArray::operator+(const APyFloat& rhs) const
         const exp_t res_max_exponent = ((1ULL << exp_bits) - 1);
         APyFloatData x;
         // Result array
-        APyFloatArray res(shape, exp_bits, man_bits);
+        APyFloatArray res(shape, exp_bits, man_bits, bias);
         APyFloatData y = rhs.get_data();
         bool y_is_max_exponent = y.exp == res_max_exponent;
 
@@ -434,9 +431,11 @@ APyFloatArray APyFloatArray::operator+(const APyFloat& rhs) const
     }
     auto new_exp_bits = std::max(exp_bits, rhs.get_exp_bits());
     auto new_man_bits = std::max(man_bits, rhs.get_man_bits());
+    const auto res_bias
+        = calc_bias(new_exp_bits, exp_bits, bias, rhs.get_exp_bits(), rhs.get_bias());
 
     // Calculate new format
-    APyFloatArray res(shape, new_exp_bits, new_man_bits);
+    APyFloatArray res(shape, new_exp_bits, new_man_bits, res_bias);
 
     APyFloat lhs_scalar(exp_bits, man_bits, bias);
     // Perform operations
@@ -498,7 +497,8 @@ APyFloatArray APyFloatArray::operator*(const APyFloatArray& rhs) const
     // Calculate new format
     const uint8_t res_exp_bits = std::max(exp_bits, rhs.exp_bits);
     const uint8_t res_man_bits = std::max(man_bits, rhs.man_bits);
-    const auto res_bias = APyFloat::ieee_bias(res_exp_bits);
+    const auto res_bias
+        = calc_bias(res_exp_bits, exp_bits, bias, rhs.exp_bits, rhs.bias);
     APyFloatArray res(shape, res_exp_bits, res_man_bits, res_bias);
     const auto quantization = get_float_quantization_mode();
     hadamard_multiplication(
@@ -660,7 +660,8 @@ APyFloatArray APyFloatArray::operator*(const APyFloat& rhs) const
     // Calculate new format
     const auto res_exp_bits = std::max(exp_bits, rhs.get_exp_bits());
     const auto res_man_bits = std::max(man_bits, rhs.get_man_bits());
-    const auto res_bias = APyFloat::ieee_bias(res_exp_bits);
+    const auto res_bias
+        = calc_bias(res_exp_bits, exp_bits, bias, rhs.get_exp_bits(), rhs.get_bias());
     APyFloatArray res(shape, res_exp_bits, res_man_bits, res_bias);
 
     const int sum_man_bits = man_bits + rhs.get_man_bits();
@@ -867,9 +868,11 @@ APyFloatArray APyFloatArray::operator/(const APyFloatArray& rhs) const
     }
 
     // Calculate new format
-    APyFloatArray res(
-        shape, std::max(exp_bits, rhs.exp_bits), std::max(man_bits, rhs.man_bits)
-    );
+    const auto res_exp_bits = std::max(exp_bits, rhs.exp_bits);
+    const auto res_man_bits = std::max(man_bits, rhs.man_bits);
+    const auto res_bias
+        = calc_bias(res_exp_bits, exp_bits, bias, rhs.exp_bits, rhs.bias);
+    APyFloatArray res(shape, res_exp_bits, res_man_bits, res_bias);
 
     APyFloat lhs_scalar(exp_bits, man_bits, bias);
     APyFloat rhs_scalar(rhs.exp_bits, rhs.man_bits, rhs.bias);
@@ -887,11 +890,11 @@ APyFloatArray APyFloatArray::operator/(const APyFloatArray& rhs) const
 APyFloatArray APyFloatArray::operator/(const APyFloat& rhs) const
 {
     // Calculate new format
-    APyFloatArray res(
-        shape,
-        std::max(exp_bits, rhs.get_exp_bits()),
-        std::max(man_bits, rhs.get_man_bits())
-    );
+    const auto res_exp_bits = std::max(exp_bits, rhs.get_exp_bits());
+    const auto res_man_bits = std::max(man_bits, rhs.get_man_bits());
+    const auto res_bias
+        = calc_bias(res_exp_bits, exp_bits, bias, rhs.get_exp_bits(), rhs.get_bias());
+    APyFloatArray res(shape, res_exp_bits, res_man_bits, res_bias);
 
     APyFloat lhs_scalar(exp_bits, man_bits, bias);
     // Perform operations
@@ -1070,12 +1073,6 @@ APyFloatArray APyFloatArray::from_double(
         return from_array(ndarray, exp_bits, man_bits, bias);
     }
 
-    if (bias.has_value() && bias.value() != APyFloat::ieee_bias(exp_bits)) {
-        throw NotImplementedException(
-            "Not implemented: APyFloatArray with non IEEE-like bias"
-        );
-    }
-
     APyFloatArray result(
         python_sequence_extract_shape(double_seq), exp_bits, man_bits, bias
     );
@@ -1109,12 +1106,6 @@ APyFloatArray APyFloatArray::from_array(
     std::optional<exp_t> bias
 )
 {
-    if (bias.has_value() && bias.value() != APyFloat::ieee_bias(exp_bits)) {
-        throw NotImplementedException(
-            "Not implemented: APyFloatArray with non IEEE-like bias"
-        );
-    }
-
     std::size_t ndim = ndarray.ndim();
     assert(ndim > 0);
     std::vector<std::size_t> shape(ndim, 0);
@@ -1296,10 +1287,14 @@ APyFloat APyFloatArray::checked_inner_product(
         const auto acc_option = accumulator_mode.value();
         const auto tmp_exp_bits = acc_option.exp_bits;
         const auto tmp_man_bits = acc_option.man_bits;
-        const auto tmp_bias = APyFloat::ieee_bias(tmp_exp_bits);
+        const auto tmp_bias = acc_option.bias.value_or(
+            calc_bias(tmp_exp_bits, exp_bits, bias, rhs.exp_bits, rhs.bias)
+        );
+        ;
         // If an accumulator is used, the operands must be resized before the
         // multiplication. This is because the products would otherwise get quantized
         // too early.
+        // NOTE: This assumes that the format of the accumulator is larger
 
         APyFloatArray hadamard(shape, tmp_exp_bits, tmp_man_bits, tmp_bias);
         // Hadamard product of `*this` and `rhs`
@@ -1320,7 +1315,7 @@ APyFloat APyFloatArray::checked_inner_product(
         sum = sum._cast(
             max_exp_bits,
             max_man_bits,
-            APyFloat::ieee_bias(max_exp_bits),
+            calc_bias_general(max_exp_bits, exp_bits, bias, rhs.exp_bits, rhs.bias),
             acc_option.quantization
         );
         return sum;
@@ -1338,7 +1333,7 @@ APyFloat APyFloatArray::vector_sum(const QuantizationMode quantization) const
 {
     // +5 to give room for leading one, carry, and 3 guard bits
     const unsigned int max_man_bits = man_bits + 5;
-    APyFloat ret(0, 0, 0, exp_bits, man_bits);
+    APyFloat ret(0, 0, 0, exp_bits, man_bits, bias);
     if ((max_man_bits <= _MAN_T_SIZE_BITS)
         && (quantization != QuantizationMode::STOCH_WEIGHTED)) {
         APyFloatData x;
@@ -1496,7 +1491,7 @@ APyFloat APyFloatArray::vector_sum(const QuantizationMode quantization) const
         ret.set_data({ sum_sign, sum_exp, sum_man });
         return ret;
     }
-    APyFloat tmp(0, 0, 0, exp_bits, man_bits);
+    APyFloat tmp(0, 0, 0, exp_bits, man_bits, bias);
 
     for (std::size_t i = 0; i < data.size(); i++) {
         tmp.set_data(data[i]);
@@ -1515,7 +1510,8 @@ APyFloatArray APyFloatArray::checked_2d_matmul(const APyFloatArray& rhs) const
         : std::vector<std::size_t> { shape[0] };              // rhs is 1-D
     const std::uint8_t max_exp_bits = std::max(exp_bits, rhs.exp_bits);
     const std::uint8_t max_man_bits = std::max(man_bits, rhs.man_bits);
-    const auto res_bias = APyFloat::ieee_bias(max_exp_bits);
+    const auto res_bias
+        = calc_bias(max_exp_bits, exp_bits, bias, rhs.exp_bits, rhs.bias);
     const auto res_cols = rhs.shape.size() > 1 ? rhs.shape[1] : 1;
 
     auto accumulator_mode = get_accumulator_mode_float();
@@ -1533,7 +1529,10 @@ APyFloatArray APyFloatArray::checked_2d_matmul(const APyFloatArray& rhs) const
         const auto acc_option = accumulator_mode.value();
         const auto tmp_exp_bits = acc_option.exp_bits;
         const auto tmp_man_bits = acc_option.man_bits;
-        const auto tmp_bias = APyFloat::ieee_bias(tmp_exp_bits);
+
+        const auto tmp_bias = acc_option.bias.value_or(
+            calc_bias(tmp_exp_bits, exp_bits, bias, rhs.exp_bits, rhs.bias)
+        );
 
         APyFloatArray hadamard({ shape[1] }, tmp_exp_bits, tmp_man_bits, tmp_bias);
         APyFloatArray casted_this
diff --git a/src/apytypes_common.cc b/src/apytypes_common.cc
index 3302ab4b2..7cd8631b0 100644
--- a/src/apytypes_common.cc
+++ b/src/apytypes_common.cc
@@ -191,7 +191,7 @@ APyFloatAccumulatorContext::APyFloatAccumulatorContext(
 
     new_mode.exp_bits = exp_bits.value();
     new_mode.man_bits = man_bits.value();
-    new_mode.bias = bias.value_or(APyFloat::ieee_bias(new_mode.exp_bits));
+    new_mode.bias = bias;
     new_mode.quantization = quantization.value_or(get_float_quantization_mode());
 
     // Set the current mode
diff --git a/src/apytypes_common.h b/src/apytypes_common.h
index 7e6b31a39..02f45b434 100644
--- a/src/apytypes_common.h
+++ b/src/apytypes_common.h
@@ -177,7 +177,7 @@ std::optional<APyFixedAccumulatorOption> get_accumulator_mode_fixed();
 struct APyFloatAccumulatorOption {
     std::uint8_t exp_bits;
     std::uint8_t man_bits;
-    exp_t bias;
+    std::optional<exp_t> bias;
     QuantizationMode quantization;
 };
 


Test patch:
diff --git a/lib/test/apyfloat/test_arithmetic.py b/lib/test/apyfloat/test_arithmetic.py
index 9a1943e3b..cd76de347 100644
--- a/lib/test/apyfloat/test_arithmetic.py
+++ b/lib/test/apyfloat/test_arithmetic.py
@@ -132,6 +132,53 @@ def test_add_inf_nan():
     assert (APyFloat(0, 0x1F, 1, 5, 7) + APyFloat(0, 0x1F, 1, 5, 7)).is_nan
 
 
+@pytest.mark.float_sub
+def test_sub_mixed_bias():
+    # Test that the implementation doesn't do "x.cast() - y.cast()"
+    x = APyFloat(sign=0, exp=30, man=0, exp_bits=5, man_bits=2, bias=14)
+    y = APyFloat(sign=0, exp=30, man=0, exp_bits=5, man_bits=2, bias=16)
+    assert (_ := x - y).is_identical(
+        APyFloat(sign=0, exp=30, man=2, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Test subtraction when two numbers are equal but with different formats
+    x = APyFloat(sign=0, exp=14, man=0, exp_bits=5, man_bits=2, bias=14)
+    y = APyFloat(sign=0, exp=21, man=0, exp_bits=5, man_bits=2, bias=21)
+    assert (_ := x - y).is_identical(
+        APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=17)
+    )
+
+
+@pytest.mark.float_add
+def test_add_mixed_bias_overflow():
+    """Test that a result can overflow to infinity due to a change in bias."""
+    x = APyFloat(sign=0, exp=21, man=0, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=25)
+
+    # Add with zero
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=31, man=0, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Add with zero but with larger bias difference
+    y = APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=27)
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=31, man=0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Add with small normal number
+    y = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=31, man=0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Add with subnormal number
+    y = APyFloat(sign=0, exp=0, man=1, exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=31, man=0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+
 @pytest.mark.float_add
 def test_add_double():
     """Some tests for sanity checking addition with binary64. This uses the faster path."""
@@ -445,6 +492,72 @@ def test_long_mul_with_one():
             assert res.is_identical(b)
 
 
+@pytest.mark.float_mul
+def test_mul_mixed_bias():
+    # Test that the implementation doesn't do "x.cast() * y.cast()"
+    x = APyFloat(sign=0, exp=30, man=0, exp_bits=5, man_bits=2, bias=14)
+    y = APyFloat(sign=0, exp=15, man=0, exp_bits=5, man_bits=2, bias=16)
+    assert (_ := x * y).is_identical(
+        APyFloat(sign=0, exp=30, man=0, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Multiply two one's but with different formats
+    x = APyFloat(sign=0, exp=14, man=0, exp_bits=5, man_bits=2, bias=14)
+    y = APyFloat(sign=0, exp=21, man=0, exp_bits=5, man_bits=2, bias=21)
+    assert (_ := x * y).is_identical(
+        APyFloat(sign=0, exp=17, man=0, exp_bits=5, man_bits=2, bias=17)
+    )
+
+
+@pytest.mark.float_mul
+def test_mul_mixed_bias_overflow():
+    """Test that a result can overflow to infinity due to a change in bias."""
+    x = APyFloat(sign=0, exp=21, man=0, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=25, man=0, exp_bits=5, man_bits=2, bias=25)
+
+    # Multiply with one
+    assert (_ := x * y).is_identical(
+        APyFloat(sign=0, exp=31, man=0, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Multiply with one but with larger bias difference
+    y = APyFloat(sign=0, exp=27, man=0, exp_bits=5, man_bits=2, bias=27)
+    assert (_ := x * y).is_identical(
+        APyFloat(sign=0, exp=31, man=0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Should become one
+    x = APyFloat(sign=0, exp=30, man=0, exp_bits=5, man_bits=2, bias=3)
+    y = APyFloat(sign=0, exp=1, man=0, exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x * y).is_identical(
+        APyFloat(sign=0, exp=15, man=0, exp_bits=5, man_bits=2)
+    )
+
+
+@pytest.mark.float_add
+def test_add_mixed_bias_underflow():
+    """Test that a result can become zero due to a change in bias."""
+    x = APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=25)
+
+    # Add with zero
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Add with zero but with larger bias difference
+    y = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=27)
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Add with subnormal number
+    y = APyFloat(sign=0, exp=0, man=1, exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+
 @pytest.mark.float_mul
 @pytest.mark.parametrize("x", [0.0, 1.0, float("inf"), float("nan")])
 @pytest.mark.parametrize("y", [0.0, 1.0, float("inf"), float("nan")])
@@ -725,6 +838,52 @@ def test_long_div():
     assert res.is_identical(APyFloat(sign=0, exp=0, man=2, exp_bits=14, man_bits=59))
 
 
+@pytest.mark.float_div
+def test_div_mixed_bias():
+    # Test that the implementation doesn't do "x.cast() / y.cast()"
+    x = APyFloat(sign=0, exp=30, man=0, exp_bits=5, man_bits=2, bias=4)
+    y = APyFloat(sign=0, exp=30, man=0, exp_bits=5, man_bits=2, bias=8)
+    assert (_ := x / y).is_identical(
+        APyFloat(sign=0, exp=10, man=0, exp_bits=5, man_bits=2, bias=6)
+    )
+
+
+@pytest.mark.float_mul
+def test_div_mixed_bias_overflow():
+    """Test that a result can overflow to infinity due to a change in bias."""
+    x = APyFloat(sign=0, exp=21, man=0, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=25, man=0, exp_bits=5, man_bits=2, bias=25)
+
+    # Divide with one
+    assert (_ := x / y).is_identical(
+        APyFloat(sign=0, exp=31, man=0, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Divide with one but with larger bias difference
+    y = APyFloat(sign=0, exp=27, man=0, exp_bits=5, man_bits=2, bias=27)
+    assert (_ := x / y).is_identical(
+        APyFloat(sign=0, exp=31, man=0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+
+@pytest.mark.float_add
+def test_div_mixed_bias_underflow():
+    """Test that a result can become zero due to a change in bias."""
+    x = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=25)
+    y = APyFloat(sign=0, exp=5, man=0, exp_bits=5, man_bits=2, bias=5)
+
+    # Divide with one
+    assert (_ := x / y).is_identical(
+        APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Divide with one but with larger bias difference
+    x = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=30)
+    assert (_ := x / y).is_identical(
+        APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=17)
+    )
+
+
 # Power
 @pytest.mark.float_pow
 def test_issue_253():
diff --git a/lib/test/apyfloat/test_comparisons.py b/lib/test/apyfloat/test_comparisons.py
index a4fca86c9..c71bd9e56 100644
--- a/lib/test/apyfloat/test_comparisons.py
+++ b/lib/test/apyfloat/test_comparisons.py
@@ -208,3 +208,98 @@ def test_signed_zero_comparison():
     assert b >= a
     assert not (b < a)
     assert b <= a
+
+
+@pytest.mark.float_comp
+def test_mixed_bias_comparisons():
+    """Comparisons where all fields are different."""
+    # x and y are equal
+    x = APyFloat(sign=0, exp=5, man=1, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=25, man=1, exp_bits=5, man_bits=2, bias=25)
+    assert x == y
+    assert x >= y
+    assert x <= y
+    assert not x > y
+    assert not x < y
+    assert not x != y
+
+    # Two normal numbers but y is smaller
+    y = APyFloat(sign=0, exp=5, man=1, exp_bits=5, man_bits=2, bias=6)
+    assert not x == y
+    assert x >= y
+    assert not x <= y
+    assert x > y
+    assert not x < y
+    assert x != y
+
+    # Two equal subnormal numbers
+    x = APyFloat(sign=0, exp=0, man=1, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=0, man=2, exp_bits=5, man_bits=2, bias=6)
+    assert x == y
+    assert x >= y
+    assert x <= y
+    assert not x > y
+    assert not x < y
+    assert not x != y
+
+    # Two subnormal numbers but y is slightly larger
+    y = APyFloat(sign=0, exp=0, man=3, exp_bits=5, man_bits=2, bias=6)
+    assert not x == y
+    assert not x >= y
+    assert x <= y
+    assert not x > y
+    assert x < y
+    assert x != y
+
+    # Subnormal number equal to normal number
+    x = APyFloat(sign=0, exp=0, man=1, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=1, man=0, exp_bits=5, man_bits=2, bias=7)
+    assert x == y
+    assert x >= y
+    assert x <= y
+    assert not x > y
+    assert not x < y
+    assert not x != y
+
+
+@pytest.mark.float_comp
+def test_mixed_format_comparisons():
+    """Comparisons with mixed formats."""
+    # x and y are equal
+    x = APyFloat(sign=0, exp=5, man=2, exp_bits=4, man_bits=3, bias=5)
+    y = APyFloat(sign=0, exp=25, man=1, exp_bits=5, man_bits=2, bias=25)
+    assert x == y
+    assert x >= y
+    assert x <= y
+    assert not x > y
+    assert not x < y
+    assert not x != y
+
+    # Two normal numbers but y is smaller
+    y = APyFloat(sign=0, exp=7, man=3, exp_bits=5, man_bits=5, bias=8)
+    assert not x == y
+    assert x >= y
+    assert not x <= y
+    assert x > y
+    assert not x < y
+    assert x != y
+
+    # Two equal subnormal numbers
+    x = APyFloat(sign=0, exp=0, man=1, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=0, man=4, exp_bits=4, man_bits=3, bias=6)
+    assert x == y
+    assert x >= y
+    assert x <= y
+    assert not x > y
+    assert not x < y
+    assert not x != y
+
+    # Two subnormal numbers but y is slightly larger
+    x = APyFloat(sign=0, exp=0, man=1, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=0, man=5, exp_bits=4, man_bits=3, bias=6)
+    assert not x == y
+    assert not x >= y
+    assert x <= y
+    assert not x > y
+    assert x < y
+    assert x != y
diff --git a/lib/test/apyfloatarray/test_arithmetic.py b/lib/test/apyfloatarray/test_arithmetic.py
index 4b0dfcb22..f27ab0cba 100644
--- a/lib/test/apyfloatarray/test_arithmetic.py
+++ b/lib/test/apyfloatarray/test_arithmetic.py
@@ -869,6 +869,93 @@ def test_array_infinity_saturation_mul():
         assert res.is_identical(APyFloatArray([1], [30], [3], 5, 2))
 
 
+@pytest.mark.float_mul
+@pytest.mark.parametrize("with_scalar", [False, True])
+def test_array_mul_mixed_bias(with_scalar):
+    # Test that the implementation doesn't do "x.cast() * y.cast()"
+    x = APyFloatArray([0], [30], [0], exp_bits=5, man_bits=2, bias=14)
+
+    if with_scalar:
+        y = APyFloat(sign=0, exp=15, man=0, exp_bits=5, man_bits=2, bias=16)
+    else:
+        y = APyFloatArray([0], [15], [0], exp_bits=5, man_bits=2, bias=16)
+    assert (_ := x * y).is_identical(
+        APyFloatArray([0], [30], [0], exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Multiply two one's but with different formats
+    x = APyFloatArray([0], [14], [0], exp_bits=5, man_bits=2, bias=14)
+
+    if with_scalar:
+        y = APyFloat(sign=0, exp=21, man=0, exp_bits=5, man_bits=2, bias=21)
+    else:
+        y = APyFloatArray([0], [21], [0], exp_bits=5, man_bits=2, bias=21)
+    assert (_ := x * y).is_identical(
+        APyFloatArray([0], [17], [0], exp_bits=5, man_bits=2, bias=17)
+    )
+
+
+@pytest.mark.float_mul
+@pytest.mark.parametrize("with_scalar", [False, True])
+def test_array_mul_mixed_bias_overflow(with_scalar):
+    """Test that a result can overflow to infinity due to a change in bias."""
+    x = APyFloatArray([0], [21], [0], exp_bits=5, man_bits=2, bias=5)
+
+    if with_scalar:
+        y = APyFloat(sign=0, exp=25, man=0, exp_bits=5, man_bits=2, bias=25)
+    else:
+        y = APyFloatArray([0], [25], [0], exp_bits=5, man_bits=2, bias=25)
+
+    # Multiply with one
+    assert (_ := x * y).is_identical(
+        APyFloatArray([0], [31], [0], exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Multiply with one but with larger bias difference
+    if with_scalar:
+        y = APyFloat(sign=0, exp=27, man=0, exp_bits=5, man_bits=2, bias=27)
+    else:
+        y = APyFloatArray([0], [27], [0], exp_bits=5, man_bits=2, bias=27)
+    assert (_ := x * y).is_identical(
+        APyFloatArray([0], [31], [0], exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Should become one
+    x = APyFloatArray([0], [30], [0], exp_bits=5, man_bits=2, bias=3)
+
+    if with_scalar:
+        y = APyFloat(sign=0, exp=1, man=0, exp_bits=5, man_bits=2, bias=28)
+    else:
+        y = APyFloatArray([0], [1], [0], exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x * y).is_identical(
+        APyFloatArray([0], [15], [0], exp_bits=5, man_bits=2)
+    )
+
+
+@pytest.mark.float_add
+def test_add_mixed_bias_underflow():
+    """Test that a result can become zero due to a change in bias."""
+    x = APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=5)
+    y = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=25)
+
+    # Add with zero
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Add with zero but with larger bias difference
+    y = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=27)
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Add with subnormal number
+    y = APyFloat(sign=0, exp=0, man=1, exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x + y).is_identical(
+        APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+
 @pytest.mark.float_add
 @pytest.mark.float_sub
 @pytest.mark.float_array
@@ -1012,3 +1099,110 @@ def test_array_infinity_saturation_div():
         # 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))
+
+
+@pytest.mark.float_add
+@pytest.mark.parametrize("with_scalar", [False, True])
+def test_array_add_mixed_bias_overflow(with_scalar):
+    """Test that a result can overflow to infinity due to a change in bias."""
+    x = APyFloatArray([0], [21], [0], exp_bits=5, man_bits=2, bias=5)
+
+    if with_scalar:
+        y = APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=25)
+    else:
+        y = APyFloatArray([0], [0], [0], exp_bits=5, man_bits=2, bias=25)
+
+    # Add with zero
+    assert (_ := x + y).is_identical(
+        APyFloatArray([0], [31], [0], exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Add with zero but with larger bias difference
+    if with_scalar:
+        y = APyFloat(sign=0, exp=0, man=0, exp_bits=5, man_bits=2, bias=27)
+    else:
+        y = APyFloatArray([0], [0], [0], exp_bits=5, man_bits=2, bias=27)
+    assert (_ := x + y).is_identical(
+        APyFloatArray([0], [31], [0], exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Add with small normal number
+    if with_scalar:
+        y = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=28)
+    else:
+        y = APyFloatArray([0], [1], [1], exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x + y).is_identical(
+        APyFloatArray([0], [31], [0], exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Add with subnormal number
+    if with_scalar:
+        y = APyFloat(sign=0, exp=0, man=1, exp_bits=5, man_bits=2, bias=28)
+    else:
+        y = APyFloatArray([0], [0], [1], exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x + y).is_identical(
+        APyFloatArray([0], [31], [0], exp_bits=5, man_bits=2, bias=16)
+    )
+
+
+@pytest.mark.float_add
+@pytest.mark.parametrize("with_scalar", [False, True])
+def test_array_add_mixed_bias_underflow(with_scalar):
+    """Test that a result can become zero due to a change in bias."""
+    x = APyFloatArray([0], [0], [0], exp_bits=5, man_bits=2, bias=5)
+
+    if with_scalar:
+        y = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=25)
+    else:
+        y = APyFloatArray([0], [1], [1], exp_bits=5, man_bits=2, bias=25)
+
+    # Add with zero
+    assert (_ := x + y).is_identical(
+        APyFloatArray([0], [0], [0], exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Add with zero but with larger bias difference
+    if with_scalar:
+        y = APyFloat(sign=0, exp=1, man=1, exp_bits=5, man_bits=2, bias=27)
+    else:
+        y = APyFloatArray([0], [1], [1], exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x + y).is_identical(
+        APyFloatArray([0], [0], [0], exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Add with subnormal number
+    if with_scalar:
+        y = APyFloat(sign=0, exp=0, man=1, exp_bits=5, man_bits=2, bias=28)
+    else:
+        y = APyFloatArray([0], [0], [1], exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x + y).is_identical(
+        APyFloatArray([0], [0], [0], exp_bits=5, man_bits=2, bias=16)
+    )
+
+
+@pytest.mark.float_sub
+@pytest.mark.parametrize("with_scalar", [False, True])
+def test_array_sub_mixed_bias(with_scalar):
+    # Test that the implementation doesn't do "x.cast() - y.cast()"
+    x = APyFloatArray([0], [30], [0], exp_bits=5, man_bits=2, bias=14)
+
+    if with_scalar:
+        y = APyFloat(sign=0, exp=30, man=0, exp_bits=5, man_bits=2, bias=16)
+    else:
+        y = APyFloatArray([0], [30], [0], exp_bits=5, man_bits=2, bias=16)
+
+    assert (_ := x - y).is_identical(
+        APyFloatArray([0], [30], [2], exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Test subtraction when two numbers are equal but with different formats
+    x = APyFloatArray([0], [14], [0], exp_bits=5, man_bits=2, bias=14)
+
+    if with_scalar:
+        y = APyFloat(sign=0, exp=21, man=0, exp_bits=5, man_bits=2, bias=21)
+    else:
+        y = APyFloatArray([0], [21], [0], exp_bits=5, man_bits=2, bias=21)
+
+    assert (_ := x - y).is_identical(
+        APyFloatArray([0], [0], [0], exp_bits=5, man_bits=2, bias=17)
+    )
diff --git a/lib/test/apyfloatarray/test_constructors.py b/lib/test/apyfloatarray/test_constructors.py
index 74e3e3b87..ebd1eb071 100644
--- a/lib/test/apyfloatarray/test_constructors.py
+++ b/lib/test/apyfloatarray/test_constructors.py
@@ -8,8 +8,6 @@ def test_constructor_raises():
         APyFloatArray([1], [5, 2], [4], 10, 10)
     with pytest.raises(ValueError, match="Inhomogeneous sequence"):
         APyFloatArray([1, 2], [5, 2], [4, "str"], 10, 10)
-    with pytest.raises(ValueError, match="Not.*implemented.*bias"):
-        APyFloatArray([1], [5], [4], 10, 10, 12)
     with pytest.raises(ValueError, match="python_sequence_extract_shape"):
         APyFloatArray(["foo"], [5], [4], 10, 10)
     with pytest.raises(ValueError, match="python_sequence_extract_shape"):
@@ -29,8 +27,6 @@ def test_constructor_raises():
         match="Non <type>/sequence found when walking: '<class 'apytypes._apytypes.APyFloatArray'>'",
     ):
         APyFloatArray([True], [4], [APyFloatArray], 10, 10)
-    with pytest.raises(ValueError, match="Not.*implemented.*bias"):
-        APyFloatArray.from_float([1, 5], 10, 10, 12)
 
 
 @pytest.mark.float_array
diff --git a/lib/test/apyfloatarray/test_matmul.py b/lib/test/apyfloatarray/test_matmul.py
index 33cbc320b..0ba5413c6 100644
--- a/lib/test/apyfloatarray/test_matmul.py
+++ b/lib/test/apyfloatarray/test_matmul.py
@@ -272,3 +272,337 @@ def test_inner_product_with_cancellation_long():
     assert (a @ b).is_identical(b @ a)
     assert (b @ a).is_identical(APyFloat(0, 488, 1152921504598458368, 10, 60))
     assert (a @ b).is_identical(sum(a * b))
+
+
+@pytest.mark.float_array
+def test_array_matmul_mixed_bias():
+    """Test that the implementation doesn't do "x.cast() @ y.cast()"."""
+    x = APyFloatArray([0], [30], [0], exp_bits=5, man_bits=2, bias=14)
+    y = APyFloatArray([0], [15], [0], exp_bits=5, man_bits=2, bias=16)
+    assert (_ := x @ y).is_identical(
+        APyFloat(0, 30, 0, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Multiply two one's but with different formats
+    x = APyFloatArray([0], [14], [0], exp_bits=5, man_bits=2, bias=14)
+    y = APyFloatArray([0], [21], [0], exp_bits=5, man_bits=2, bias=21)
+    assert (_ := x @ y).is_identical(
+        APyFloat(0, 17, 0, exp_bits=5, man_bits=2, bias=17)
+    )
+
+    # Same tests again but with (2, 2) matrices
+    x = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[30, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=14,
+    )
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[15, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=16,
+    )
+    assert (_ := x @ y.T).is_identical(
+        APyFloatArray(
+            [[0, 0], [0, 0]],
+            [[30, 0], [0, 0]],
+            [[0, 0], [0, 0]],
+            exp_bits=5,
+            man_bits=2,
+            bias=15,
+        )
+    )
+
+    x = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[14, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=14,
+    )
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[21, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=21,
+    )
+    assert (_ := x @ y.T).is_identical(
+        APyFloatArray(
+            [[0, 0], [0, 0]],
+            [[17, 0], [0, 0]],
+            [[0, 0], [0, 0]],
+            exp_bits=5,
+            man_bits=2,
+            bias=17,
+        )
+    )
+
+
+@pytest.mark.float_array
+def test_array_matmul_mixed_bias_overflow():
+    """Test that a result can overflow to infinity due to a change in bias."""
+    x = APyFloatArray([0], [21], [0], exp_bits=5, man_bits=2, bias=5)
+    y = APyFloatArray([0], [25], [0], exp_bits=5, man_bits=2, bias=25)
+
+    # Multiply with one
+    assert (_ := x @ y).is_identical(
+        APyFloat(0, 31, 0, exp_bits=5, man_bits=2, bias=15)
+    )
+
+    # Multiply with one but with larger bias difference
+    y = APyFloatArray([0], [27], [0], exp_bits=5, man_bits=2, bias=27)
+    assert (_ := x @ y).is_identical(
+        APyFloat(0, 31, 0, exp_bits=5, man_bits=2, bias=16)
+    )
+
+    # Should become one
+    x = APyFloatArray([0], [30], [0], exp_bits=5, man_bits=2, bias=3)
+    y = APyFloatArray([0], [1], [0], exp_bits=5, man_bits=2, bias=28)
+    assert (_ := x @ y).is_identical(APyFloat(0, 15, 0, exp_bits=5, man_bits=2))
+
+    # Same tests again but now (2, 2) matrices
+    x = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[21, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=5,
+    )
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[25, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=25,
+    )
+    assert (_ := x @ y.T).is_identical(
+        APyFloatArray(
+            [[0, 0], [0, 0]],
+            [[31, 0], [0, 0]],
+            [[0, 0], [0, 0]],
+            exp_bits=5,
+            man_bits=2,
+            bias=15,
+        )
+    )
+
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[27, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=27,
+    )
+    assert (_ := x @ y.T).is_identical(
+        APyFloatArray(
+            [[0, 0], [0, 0]],
+            [[31, 0], [0, 0]],
+            [[0, 0], [0, 0]],
+            exp_bits=5,
+            man_bits=2,
+            bias=16,
+        )
+    )
+
+    x = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[30, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=3,
+    )
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[1, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=28,
+    )
+    assert (_ := x @ y.T).is_identical(
+        APyFloatArray(
+            [[0, 0], [0, 0]],
+            [[15, 0], [0, 0]],
+            [[0, 0], [0, 0]],
+            exp_bits=5,
+            man_bits=2,
+        )
+    )
+
+
+@pytest.mark.xfail()
+@pytest.mark.float_array
+def test_matrix_mul_acc_context_mixed_bias():
+    """Test that the accumulator doesn't cast prematurely. All test cases here use an accumulator that shouldn't have any effect."""
+    x = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[30, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=14,
+    )
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[15, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=16,
+    )
+
+    with APyFloatAccumulatorContext(exp_bits=5, man_bits=2, bias=15):
+        assert (_ := x @ y.T).is_identical(
+            APyFloatArray(
+                [[0, 0], [0, 0]],
+                [[30, 0], [0, 0]],
+                [[0, 0], [0, 0]],
+                exp_bits=5,
+                man_bits=2,
+                bias=15,
+            )
+        )
+
+    x = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[14, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=14,
+    )
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[21, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=21,
+    )
+
+    with APyFloatAccumulatorContext(exp_bits=5, man_bits=2, bias=17):
+        assert (_ := x @ y.T).is_identical(
+            APyFloatArray(
+                [[0, 0], [0, 0]],
+                [[17, 0], [0, 0]],
+                [[0, 0], [0, 0]],
+                exp_bits=5,
+                man_bits=2,
+                bias=17,
+            )
+        )
+
+
+@pytest.mark.xfail()
+@pytest.mark.float_array
+def test_matrix_mul_acc_context_mixed_bias_overflow():
+    """Test that the accumulator doesn't cast prematurely. All test cases here use an accumulator that shouldn't have any effect."""
+    x = APyFloatArray([0], [21], [0], exp_bits=5, man_bits=2, bias=5)
+    y = APyFloatArray([0], [25], [0], exp_bits=5, man_bits=2, bias=25)
+
+    # Multiply with one
+    with APyFloatAccumulatorContext(exp_bits=5, man_bits=2, bias=15):
+        assert (_ := x @ y).is_identical(
+            APyFloat(0, 31, 0, exp_bits=5, man_bits=2, bias=15)
+        )
+
+    # Multiply with one but with larger bias difference
+    y = APyFloatArray([0], [27], [0], exp_bits=5, man_bits=2, bias=27)
+    with APyFloatAccumulatorContext(exp_bits=5, man_bits=2, bias=16):
+        assert (_ := x @ y).is_identical(
+            APyFloat(0, 31, 0, exp_bits=5, man_bits=2, bias=16)
+        )
+
+    # Should become one
+    x = APyFloatArray([0], [30], [0], exp_bits=5, man_bits=2, bias=3)
+    y = APyFloatArray([0], [1], [0], exp_bits=5, man_bits=2, bias=28)
+    with APyFloatAccumulatorContext(exp_bits=5, man_bits=2):
+        assert (_ := x @ y).is_identical(APyFloat(0, 15, 0, exp_bits=5, man_bits=2))
+
+    # Same tests again but now (2, 2) matrices
+    x = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[21, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=5,
+    )
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[25, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=25,
+    )
+    with APyFloatAccumulatorContext(exp_bits=5, man_bits=2):
+        assert (_ := x @ y.T).is_identical(
+            APyFloatArray(
+                [[0, 0], [0, 0]],
+                [[31, 0], [0, 0]],
+                [[0, 0], [0, 0]],
+                exp_bits=5,
+                man_bits=2,
+                bias=15,
+            )
+        )
+
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[27, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=27,
+    )
+    with APyFloatAccumulatorContext(exp_bits=5, man_bits=2, bias=16):
+        assert (_ := x @ y.T).is_identical(
+            APyFloatArray(
+                [[0, 0], [0, 0]],
+                [[31, 0], [0, 0]],
+                [[0, 0], [0, 0]],
+                exp_bits=5,
+                man_bits=2,
+                bias=16,
+            )
+        )
+
+    x = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[30, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=3,
+    )
+    y = APyFloatArray(
+        [[0, 0], [0, 0]],
+        [[1, 0], [0, 0]],
+        [[0, 0], [0, 0]],
+        exp_bits=5,
+        man_bits=2,
+        bias=28,
+    )
+    with APyFloatAccumulatorContext(exp_bits=5, man_bits=2):
+        assert (_ := x @ y.T).is_identical(
+            APyFloatArray(
+                [[0, 0], [0, 0]],
+                [[15, 0], [0, 0]],
+                [[0, 0], [0, 0]],
+                exp_bits=5,
+                man_bits=2,
+            )
+        )