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 1adbe8082..ea9a54824 100644
--- a/src/apyfloat.cc
+++ b/src/apyfloat.cc
@@ -670,7 +670,14 @@ APyFloat APyFloat::operator+(const APyFloat& rhs) const
     APyFloat res, x, y;
     // Handle the zero cases, other special cases are further down
     if (same_type_as(rhs)) {
+        // Sign is currently not handled correctly if both operands are zero
         if (is_zero()) {
+            if (rhs.is_zero()) {
+                const bool new_sign = (sign == rhs.sign)
+                    ? sign
+                    : get_float_quantization_mode() == QuantizationMode::TRN;
+                return construct_zero(new_sign);
+            }
             return rhs;
         }
         if (rhs.is_zero()) {
@@ -694,7 +701,8 @@ APyFloat APyFloat::operator+(const APyFloat& rhs) const
         if (exp < rhs.exp || (exp == rhs.exp && man < rhs.man)) {
             std::swap(x, y);
         } else if (sign != rhs.sign && exp == rhs.exp && man == rhs.man) {
-            res.set_to_zero(true);
+            // +0 for all quantization modes except TO_NEG
+            res.set_to_zero(get_float_quantization_mode() == QuantizationMode::TRN);
             return res;
         }
     } else {
@@ -704,13 +712,20 @@ APyFloat APyFloat::operator+(const APyFloat& rhs) const
         // 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);
+        res = APyFloat(res_exp_bits, res_man_bits, res_bias);
+
         if (is_zero()) {
+            if (y.is_zero()) {
+                const bool new_sign = (sign == y.sign)
+                    ? sign
+                    : get_float_quantization_mode() == QuantizationMode::TRN;
+                return res.construct_zero(new_sign);
+            }
             return y;
         }
         if (rhs.is_zero()) {
             return x;
         }
-        res = APyFloat(res_exp_bits, res_man_bits, res_bias);
         // Handle the NaN and inf cases
         if (is_max_exponent() || rhs.is_max_exponent()) {
             if (is_nan() || rhs.is_nan()
@@ -725,7 +740,8 @@ APyFloat APyFloat::operator+(const APyFloat& rhs) const
         if (x.abs() < y.abs()) {
             std::swap(x, y);
         } else if (sign != rhs.sign && x.exp == y.exp && x.man == y.man) {
-            res.set_to_zero(true);
+            // +0 for all quantization modes except TO_NEG
+            res.set_to_zero(get_float_quantization_mode() == QuantizationMode::TRN);
             return res;
         }
     }
@@ -889,9 +905,17 @@ APyFloat& APyFloat::operator+=(const APyFloat& rhs)
 
     // Handle the zero cases, other special cases are further down
     if (is_zero()) {
-        sign = rhs.sign;
-        man = rhs.man;
-        exp = rhs.exp;
+        if (rhs.is_zero()) {
+            sign = (sign == rhs.sign)
+                ? sign
+                : get_float_quantization_mode() == QuantizationMode::TRN;
+            exp = 0;
+            man = 0;
+        } else {
+            sign = rhs.sign;
+            man = rhs.man;
+            exp = rhs.exp;
+        }
         return *this;
     }
     if (rhs.is_zero()) {
@@ -916,17 +940,16 @@ APyFloat& APyFloat::operator+=(const APyFloat& rhs)
     bool same_sign = sign == rhs.sign;
     const APyFloat* x = &*this;
     const APyFloat* y = &rhs;
+    const auto quantization = get_float_quantization_mode();
     // Compute sign and swap operands if need to make sure |x| >= |y|
     if (exp < rhs.exp || (exp == rhs.exp && man < rhs.man)) {
         sign = rhs.sign;
         std::swap(x, y);
     } else if (!same_sign && exp == rhs.exp && man == rhs.man) {
-        set_to_zero(true);
+        set_to_zero(quantization == QuantizationMode::TRN);
         return *this;
     }
 
-    const auto quantization = get_float_quantization_mode();
-
     // Compute smaller exponent so that one can overwrite exp later
     exp_t smaller_exp = y->exp + y->is_subnormal();
 
diff --git a/src/apyfloatarray.cc b/src/apyfloatarray.cc
index 3546ec3a6..27d1000ae 100644
--- a/src/apyfloatarray.cc
+++ b/src/apyfloatarray.cc
@@ -93,6 +93,7 @@ APyFloatArray APyFloatArray::operator+(const APyFloatArray& rhs) const
     }
 
     const auto quantization = get_float_quantization_mode();
+    const bool is_to_neg = quantization == QuantizationMode::TRN;
     // +5 to give room for leading one, carry, and 3 guard bits
     const unsigned int max_man_bits = man_bits + 5;
     if (same_type_as(rhs) && (max_man_bits <= _MAN_T_SIZE_BITS)
@@ -114,12 +115,15 @@ APyFloatArray APyFloatArray::operator+(const APyFloatArray& rhs) const
             x = data[i];
             y = rhs.data[i];
             bool x_is_zero_exponent = (x.exp == 0);
+            bool y_is_zero_exponent = (y.exp == 0);
             // Handle zero cases
             if (x_is_zero_exponent && x.man == 0) {
+                if (y_is_zero_exponent && y.man == 0) {
+                    y.sign = (x.sign == y.sign) ? x.sign : is_to_neg;
+                }
                 res.data[i] = y;
                 continue;
             }
-            bool y_is_zero_exponent = (y.exp == 0);
             if (y_is_zero_exponent && y.man == 0) {
                 res.data[i] = x;
                 continue;
@@ -162,7 +166,8 @@ APyFloatArray APyFloatArray::operator+(const APyFloatArray& rhs) const
                 std::swap(x_is_zero_exponent, y_is_zero_exponent);
             } else if (x.sign != y.sign && x.exp == y.exp && x.man == y.man) {
                 // Set to zero
-                res.data[i] = { true, static_cast<exp_t>(0), static_cast<man_t>(0) };
+                res.data[i]
+                    = { is_to_neg, static_cast<exp_t>(0), static_cast<man_t>(0) };
                 continue;
             }
 
@@ -258,13 +263,10 @@ APyFloatArray APyFloatArray::operator+(const APyFloat& rhs) const
     // +5 to give room for leading one, carry, and 3 guard bits
     const unsigned int max_man_bits = man_bits + 5;
     const bool same_types = same_type_as(rhs);
+    const bool is_to_neg = quantization == QuantizationMode::TRN;
     if (same_types && (max_man_bits <= _MAN_T_SIZE_BITS)
         && (quantization != QuantizationMode::STOCH_WEIGHTED)) {
 
-        // Handle zero
-        if (rhs.is_zero()) {
-            return *this;
-        }
         const exp_t res_max_exponent = ((1ULL << exp_bits) - 1);
         APyFloatData x;
         // Result array
@@ -308,7 +310,11 @@ APyFloatArray APyFloatArray::operator+(const APyFloat& rhs) const
             x = data[i];
             // Handle zero case
             if (x.exp == 0 && x.man == 0) {
-                res.data[i] = y;
+                if (y_is_zero_exponent && y.man == 0) {
+                    res.data[i].sign = (x.sign == y.sign) ? x.sign : is_to_neg;
+                } else {
+                    res.data[i] = y;
+                }
                 continue;
             }
             bool x_is_max_exponent = x.exp == res_max_exponent;
@@ -331,7 +337,8 @@ APyFloatArray APyFloatArray::operator+(const APyFloat& rhs) const
             // Check if sum is zero
             if (x.sign != y.sign && x.exp == y.exp && x.man == y.man) {
                 // Set to zero
-                res.data[i] = { false, static_cast<exp_t>(0), static_cast<man_t>(0) };
+                res.data[i]
+                    = { is_to_neg, static_cast<exp_t>(0), static_cast<man_t>(0) };
                 continue;
             }
 
@@ -427,12 +434,6 @@ 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());
-    if (rhs.is_zero()) {
-        if (same_types) {
-            return *this;
-        }
-        return cast_no_quant(new_exp_bits, new_man_bits);
-    }
 
     // Calculate new format
     APyFloatArray res(shape, new_exp_bits, new_man_bits);


Test patch:
diff --git a/lib/test/apyfloat/run_berkeley_cases.py b/lib/test/apyfloat/run_berkeley_cases.py
index cb38d299f..26607b13d 100644
--- a/lib/test/apyfloat/run_berkeley_cases.py
+++ b/lib/test/apyfloat/run_berkeley_cases.py
@@ -14,6 +14,7 @@
     APyFloatQuantizationContext,
 )
 import itertools
+import math
 import os
 import random
 import subprocess
@@ -324,9 +325,10 @@ def print_summary(summary, output_file: str, seed: int, level: int) -> None:
                 )
                 failed_tests += summary[op][quant][0]
 
-    success_rate = (1 - float(failed_tests) / total_tests) * 100
+    # Round percentage down so it doesn't show 100% success even if something failed
+    success_rate = math.floor((1 - float(failed_tests) / total_tests) * 10000) / 100.0
     print(
-        f"Ran {total_tests} tests in total whereof {failed_tests} failed, giving an overall success rate of {success_rate:.2f}%."
+        f"Ran {total_tests} tests in total whereof {failed_tests} failed, giving an overall success rate of {success_rate}%."
     )
 
 
diff --git a/lib/test/apyfloat/test_arithmetic.py b/lib/test/apyfloat/test_arithmetic.py
index 7a6fd33be..9a1943e3b 100644
--- a/lib/test/apyfloat/test_arithmetic.py
+++ b/lib/test/apyfloat/test_arithmetic.py
@@ -4,70 +4,6 @@
 from apytypes import APyFloat
 
 
-# Exhaustive tests, long running
-# These should be replaced by testing against a third-party tool
-# @pytest.mark.slow
-# @pytest.mark.float_add
-# @pytest.mark.parametrize("x_bits", range(0, 2**8))
-# @pytest.mark.parametrize("y_bits", range(0, 2**8))
-# def test_add_all(x_bits, y_bits):
-#    y = APyFloat.from_bits(y_bits, 4, 3)
-#    x = APyFloat.from_bits(x_bits, 4, 3)
-#    ans = x + y
-#    ref = APyFloat.from_float(float(x) + float(y), 4, 3)
-#    print(f'{float(x)} + {float(y)} = {float(ref)} != {float(ans)}')
-#    print(f'{x!r} + {y!r} = {ref!r} != {ans!r}')
-#    assert ans == ref or (ans.is_nan and ref.is_nan) or (ans == 0 and ref == 0)
-
-
-# @pytest.mark.slow
-# @pytest.mark.float_sub
-# @pytest.mark.parametrize("x_bits", range(0, 2**8))
-# @pytest.mark.parametrize("y_bits", range(0, 2**8))
-# def test_sub_all(x_bits, y_bits):
-#    y = APyFloat.from_bits(y_bits, 4, 3)
-#    x = APyFloat.from_bits(x_bits, 4, 3)
-#    ans = x - y
-#    ref = APyFloat.from_float(float(x) - float(y), 4, 3)
-#    print(f'{float(x)} - {float(y)} = {float(ref)} != {float(ans)}')
-#    print(f'{x!r} - {y!r} = {ref!r} != {ans!r}')
-#    assert ans == ref or (ans.is_nan and ref.is_nan) or (ans == 0 and ref == 0)
-
-
-# @pytest.mark.slow
-# @pytest.mark.float_mul
-# @pytest.mark.parametrize("x_bits", range(0, 2**8))
-# @pytest.mark.parametrize("y_bits", range(0, 2**8))
-# def test_mul_all(x_bits, y_bits):
-#    y = APyFloat.from_bits(y_bits, 4, 3)
-#    x = APyFloat.from_bits(x_bits, 4, 3)
-#    ans = x * y
-#    prod = float(x) * float(y)
-#    ref = APyFloat.from_float(prod, 4, 3)
-#    print(f'{float(x)} * {float(y)} = {float(ref)}({prod}) != {float(ans)}')
-#    print(f'{x!r} * {y!r} = {ref!r} != {ans!r}')
-#    assert ans.is_identical(ref) or (ans.is_nan and ref.is_nan) or (ans == 0 and ref == 0)
-
-
-# @pytest.mark.slow
-# @pytest.mark.float_div
-# @pytest.mark.parametrize("x_bits", range(0, 2 * 8))
-# @pytest.mark.parametrize("y_bits", range(0, 2**8))
-# def test_div_all(x_bits, y_bits):
-#     y = APyFloat.from_bits(y_bits, 4, 3)
-#     if y != 0:
-#         x = APyFloat.from_bits(x_bits, 4, 3)
-#         ans = x / y
-#         ref = APyFloat.from_float(float(x) / float(y), 4, 3)
-#         print(f"{float(x)} / {float(y)} = {float(ref)} != {float(ans)}")
-#         print(f"{x!r} / {y!r} = (ref){ref!r} != (ans){ans!r}")
-#         assert (
-#             ans.is_identical(ref)
-#             or (ans.is_nan and ref.is_nan)
-#             or (ans == 0 and ref == 0)
-#         )
-
-
 # Negation
 @pytest.mark.parametrize("sign", ["-", ""])
 @pytest.mark.parametrize(
@@ -301,7 +237,7 @@ def test_long_add():
     x = APyFloat(0, 1046, 1232143, 11, 60)
     y = APyFloat(1, 1046, 1232143, 11, 60)
     res = x + y
-    assert res.is_identical(APyFloat(1, 0, 0, 11, 60))
+    assert res.is_identical(APyFloat(0, 0, 0, 11, 60))
 
 
 # Subtraction
@@ -355,11 +291,11 @@ def test_sub_diff_sign():
 def test_sub_res_zero():
     # 12 - 12, same format
     res = APyFloat(0, 258, 8, 9, 4) - APyFloat(0, 258, 8, 9, 4)
-    assert res.is_identical(APyFloat(1, 0, 0, 9, 4))
+    assert res.is_identical(APyFloat(0, 0, 0, 9, 4))
 
     # 12 - 12, different format
     res = APyFloat(0, 258, 4, 9, 3) - APyFloat(0, 258, 8, 9, 4)
-    assert res.is_identical(APyFloat(1, 0, 0, 9, 4))
+    assert res.is_identical(APyFloat(0, 0, 0, 9, 4))
 
 
 @pytest.mark.float_sub
diff --git a/lib/test/apyfloat/test_rounding.py b/lib/test/apyfloat/test_rounding.py
index 2c06efcf4..df144ac8f 100644
--- a/lib/test/apyfloat/test_rounding.py
+++ b/lib/test/apyfloat/test_rounding.py
@@ -23,6 +23,47 @@ def test_issue_245():
         assert res.is_identical(APyFloat(0, 0, 0, 5, 2))
 
 
+@pytest.mark.float_add
+@pytest.mark.float_sub
+@pytest.mark.parametrize("man", [10, 60])
+def test_add_sub_zero_sign(man):
+    pos_zero = APyFloat(0, 0, 0, 5, 10)
+    neg_zero = APyFloat(1, 0, 0, 5, man)
+    non_zero = APyFloat.from_float(1.0, 5, man)
+    modes = (
+        QuantizationMode.TIES_EVEN,
+        QuantizationMode.TO_ZERO,
+        QuantizationMode.TO_POS,
+    )
+
+    for mode in modes:
+        with APyFloatQuantizationContext(mode):
+            assert (_ := pos_zero + pos_zero).sign == False
+            assert (_ := pos_zero + neg_zero).sign == False
+            assert (_ := neg_zero + pos_zero).sign == False
+            assert (_ := neg_zero + neg_zero).sign == True
+
+            assert (_ := pos_zero - pos_zero).sign == False
+            assert (_ := pos_zero - neg_zero).sign == False
+            assert (_ := neg_zero - pos_zero).sign == True
+            assert (_ := neg_zero - neg_zero).sign == False
+
+            assert (_ := non_zero + (-non_zero)).sign == False
+
+    with APyFloatQuantizationContext(QuantizationMode.TO_NEG):
+        assert (_ := pos_zero + pos_zero).sign == False
+        assert (_ := pos_zero + neg_zero).sign == True
+        assert (_ := neg_zero + pos_zero).sign == True
+        assert (_ := neg_zero + neg_zero).sign == True
+
+        assert (_ := pos_zero - pos_zero).sign == True
+        assert (_ := pos_zero - neg_zero).sign == False
+        assert (_ := neg_zero - pos_zero).sign == True
+        assert (_ := neg_zero - neg_zero).sign == True
+
+        assert (_ := non_zero - non_zero).sign == True
+
+
 @pytest.mark.float_add
 class TestAPyFloatQuantizationAddSub:
     """
diff --git a/lib/test/apyfloatarray/test_arithmetic.py b/lib/test/apyfloatarray/test_arithmetic.py
index 0be4c8803..4b0dfcb22 100644
--- a/lib/test/apyfloatarray/test_arithmetic.py
+++ b/lib/test/apyfloatarray/test_arithmetic.py
@@ -48,7 +48,7 @@ def test_array_add():
     ans = APyFloatArray.from_float(
         [
             [7.0, 4.5, 10.5],
-            [-8.0, 15.0, -0.0],
+            [-8.0, 15.0, 0.0],
         ],
         exp_bits=8,
         man_bits=8,
@@ -246,7 +246,7 @@ def test_array_sub():
     ans = APyFloatArray.from_float(
         [
             [-5, -0.5, -4.5],
-            [-0.0, -25.0, -12.0],
+            [0.0, -25.0, -12.0],
         ],
         exp_bits=8,
         man_bits=8,
@@ -869,6 +869,73 @@ def test_array_infinity_saturation_mul():
         assert res.is_identical(APyFloatArray([1], [30], [3], 5, 2))
 
 
+@pytest.mark.float_add
+@pytest.mark.float_sub
+@pytest.mark.float_array
+@pytest.mark.parametrize("man", [10, 60])
+@pytest.mark.parametrize("with_scalar", [False, True])
+def test_array_add_sub_zero_sign(man, with_scalar):
+    pos_zero = APyFloatArray([0], [0], [0], 5, 10)
+    neg_zero = APyFloatArray([1], [0], [0], 5, man)
+    non_zero = APyFloatArray.from_float([1.0], 5, man)
+
+    # Scalar equivalents
+    pos_zero_scal = APyFloat(0, 0, 0, 5, 10)
+    neg_zero_scal = APyFloat(1, 0, 0, 5, man)
+    non_zero_scal = APyFloat.from_float(1.0, 5, man)
+
+    modes = (
+        QuantizationMode.TIES_EVEN,
+        QuantizationMode.TO_ZERO,
+        QuantizationMode.TO_POS,
+    )
+
+    for mode in modes:
+        with APyFloatQuantizationContext(mode):
+            rhs = pos_zero_scal if with_scalar else pos_zero
+            assert (_ := pos_zero + rhs)[0].sign == False
+            rhs = neg_zero_scal if with_scalar else neg_zero
+            assert (_ := pos_zero + rhs)[0].sign == False
+            rhs = pos_zero_scal if with_scalar else pos_zero
+            assert (_ := neg_zero + rhs)[0].sign == False
+            rhs = neg_zero_scal if with_scalar else neg_zero
+            assert (_ := neg_zero + rhs)[0].sign == True
+
+            rhs = pos_zero_scal if with_scalar else pos_zero
+            assert (_ := pos_zero - rhs)[0].sign == False
+            rhs = neg_zero_scal if with_scalar else neg_zero
+            assert (_ := pos_zero - rhs)[0].sign == False
+            rhs = pos_zero_scal if with_scalar else pos_zero
+            assert (_ := neg_zero - rhs)[0].sign == True
+            rhs = neg_zero_scal if with_scalar else neg_zero
+            assert (_ := neg_zero - rhs)[0].sign == False
+
+            rhs = non_zero_scal if with_scalar else non_zero
+            assert (_ := non_zero + (-rhs))[0].sign == False
+
+    with APyFloatQuantizationContext(QuantizationMode.TO_NEG):
+        rhs = pos_zero_scal if with_scalar else pos_zero
+        assert (_ := pos_zero + rhs)[0].sign == False
+        rhs = neg_zero_scal if with_scalar else neg_zero
+        assert (_ := pos_zero + rhs)[0].sign == True
+        rhs = pos_zero_scal if with_scalar else pos_zero
+        assert (_ := neg_zero + rhs)[0].sign == True
+        rhs = neg_zero_scal if with_scalar else neg_zero
+        assert (_ := neg_zero + rhs)[0].sign == True
+
+        rhs = pos_zero_scal if with_scalar else pos_zero
+        assert (_ := pos_zero - pos_zero)[0].sign == True
+        rhs = neg_zero_scal if with_scalar else neg_zero
+        assert (_ := pos_zero - rhs)[0].sign == False
+        rhs = pos_zero_scal if with_scalar else pos_zero
+        assert (_ := neg_zero - rhs)[0].sign == True
+        rhs = neg_zero_scal if with_scalar else neg_zero
+        assert (_ := neg_zero - rhs)[0].sign == True
+
+        rhs = non_zero_scal if with_scalar else non_zero
+        assert (_ := non_zero - rhs)[0].sign == True
+
+
 @pytest.mark.float_add
 def test_array_infinity_saturation_add():
     # Array x Scalar