Repository: apytypes/apytypes (language: python)

Issue:
Subnormal multiplication not working
```
from apytypes import APyFloat
a = APyFloat(0, 0, 31, 4, 5)
b = APyFloat(0, 0, 5, 4, 5)
```

`a*b` returns `APyFloat(sign=0, exp=0, man=1, exp_bits=4, man_bits=5)`, but this should be zero(?).

`float(a)*float(b)` is 3.695487976074219e-05, while `float(a*b)` is 0.00048828125. The first should not round to the second.


Discussion / hints:

Fix patch:
diff --git a/src/apyfloat.cc b/src/apyfloat.cc
index 8b92c246a..bc7cd5458 100644
--- a/src/apyfloat.cc
+++ b/src/apyfloat.cc
@@ -1020,18 +1020,21 @@ APyFloat APyFloat::operator*(const APyFloat& y) const
             new_man <<= 1;
         } else if (new_man & one_before) {
             // Align with longer result
-            if (tmp_exp == 0) {
-                new_man <<= 1;
-            } else {
-                new_man <<= 2;
-            }
+            new_man <<= 2;
         } else {
-            // One or two of the operands were subnormal
-            new_man <<= std::max(tmp_exp + 1, (std::int64_t)0);
-            tmp_exp = 0;
+            // One or two of the operands were subnormal.
+            // If the exponent is positive, the result is normalized by
+            // left-shifting until the exponent is zero or the mantissa is 1.xx
+            const int leading_zeros = 1 + sum_man_bits - bit_width(new_man);
+            const int shift = std::max(
+                std::min(tmp_exp, (std::int64_t)leading_zeros), (std::int64_t)0
+            );
+            tmp_exp -= shift;
+            // + 2 to align with longer result
+            new_man <<= shift + 2;
         }
 
-        if (tmp_exp < 0) {
+        if (tmp_exp <= 0) {
             new_man = (new_man >> (-tmp_exp + 1))
                 | ((new_man & ((1 << (-tmp_exp + 1)) - 1)) != 0);
             tmp_exp = 0;
diff --git a/src/apyfloatarray.cc b/src/apyfloatarray.cc
index ac75bcd8d..6439736f6 100644
--- a/src/apyfloatarray.cc
+++ b/src/apyfloatarray.cc
@@ -570,18 +570,21 @@ APyFloatArray APyFloatArray::operator*(const APyFloatArray& rhs) const
                 new_man <<= 1;
             } else if (new_man & one_before) {
                 // Align with longer result
-                if (tmp_exp == 0) {
-                    new_man <<= 1;
-                } else {
-                    new_man <<= 2;
-                }
+                new_man <<= 2;
             } else {
-                // One or two of the operands were subnormal
-                new_man <<= std::max(tmp_exp + 1, (std::int64_t)0);
-                tmp_exp = 0;
+                // One or two of the operands were subnormal.
+                // If the exponent is positive, the result is normalized by
+                // left-shifting until the exponent is zero or the mantissa is 1.xx
+                const int leading_zeros = 1 + sum_man_bits - bit_width(new_man);
+                const int shift = std::max(
+                    std::min(tmp_exp, (std::int64_t)leading_zeros), (std::int64_t)0
+                );
+                tmp_exp -= shift;
+                // + 2 to align with longer result
+                new_man <<= shift + 2;
             }
 
-            if (tmp_exp < 0) {
+            if (tmp_exp <= 0) {
                 // Shift and add sticky bit
                 new_man = (new_man >> (-tmp_exp + 1))
                     | ((new_man & ((1 << (-tmp_exp + 1)) - 1)) != 0);
@@ -762,18 +765,21 @@ APyFloatArray APyFloatArray::operator*(const APyFloat& rhs) const
                 new_man <<= 1;
             } else if (new_man & one_before) {
                 // Align with longer result
-                if (tmp_exp == 0) {
-                    new_man <<= 1;
-                } else {
-                    new_man <<= 2;
-                }
+                new_man <<= 2;
             } else {
-                // One or two of the operands were subnormal
-                new_man <<= std::max(tmp_exp + 1, (std::int64_t)0);
-                tmp_exp = 0;
+                // One or two of the operands were subnormal.
+                // If the exponent is positive, the result is normalized by
+                // left-shifting until the exponent is zero or the mantissa is 1.xx
+                const int leading_zeros = 1 + sum_man_bits - bit_width(new_man);
+                const int shift = std::max(
+                    std::min(tmp_exp, (std::int64_t)leading_zeros), (std::int64_t)0
+                );
+                tmp_exp -= shift;
+                // + 2 to align with longer result
+                new_man <<= shift + 2;
             }
 
-            if (tmp_exp < 0) {
+            if (tmp_exp <= 0) {
                 // Shift and add sticky bit
                 new_man = (new_man >> (-tmp_exp + 1))
                     | ((new_man & ((1 << (-tmp_exp + 1)) - 1)) != 0);


Test patch:
diff --git a/lib/test/apyfloat/test_arithmetic.py b/lib/test/apyfloat/test_arithmetic.py
index 5dce81584..57a7a51a5 100644
--- a/lib/test/apyfloat/test_arithmetic.py
+++ b/lib/test/apyfloat/test_arithmetic.py
@@ -452,6 +452,12 @@ def test_mul_underflow():
     res = APyFloat(0, 1, 1, 5, 2) * APyFloat(0, 1, 3, 5, 3)
     assert res == 0
 
+    res = APyFloat(0, 0, 31, 4, 5) * APyFloat(0, 0, 5, 4, 5)
+    assert res == 0
+
+    res = APyFloat(0, 0, 31, 4, 5) * APyFloat(0, 1, 1, 4, 5)
+    assert res == 0
+
 
 @pytest.mark.float_mul
 def test_mul_with_one():
@@ -515,6 +521,32 @@ def test_mul_short_subnormal():
     res = APyFloat(0, 8, 0, 4, 3) * APyFloat(0, 0, 2, 4, 3)
     assert res.is_identical(APyFloat(0, 0, 4, 4, 3))
 
+    # 4.0 x subn
+    res = APyFloat(0, 9, 0, 4, 3) * APyFloat(0, 0, 1, 4, 3)
+    assert res.is_identical(APyFloat(0, 0, 4, 4, 3))
+
+    # subnormal becoming normal
+    res = APyFloat(0, 10, 0, 4, 3) * APyFloat(0, 0, 1, 4, 3)
+    assert res.is_identical(APyFloat(0, 1, 0, 4, 3))
+
+    # 0.5 x subn
+    res = APyFloat(0, 6, 0, 4, 3) * APyFloat(0, 0, 4, 4, 3)
+    assert res.is_identical(APyFloat(0, 0, 2, 4, 3))
+
+    # 0.25 x subn
+    res = APyFloat(0, 5, 0, 4, 3) * APyFloat(0, 0, 4, 4, 3)
+    assert res.is_identical(APyFloat(0, 0, 1, 4, 3))
+
+    # 0.125 x subn
+    res = APyFloat(0, 4, 0, 4, 3) * APyFloat(0, 0, 4, 4, 3)
+    assert res.is_identical(APyFloat(0, 0, 0, 4, 3))
+
+    # Two normal numbers that generate carry, but with subnormal result
+    a = APyFloat(sign=1, exp=2, man=957, exp_bits=5, man_bits=10)
+    b = APyFloat(sign=0, exp=12, man=1015, exp_bits=5, man_bits=10)
+    res = a * b
+    assert res.is_identical(APyFloat(sign=1, exp=0, man=986, exp_bits=5, man_bits=10))
+
 
 @pytest.mark.float_mul
 def test_long_mul():
@@ -592,12 +624,16 @@ def test_div_overflow():
     """Test that a division can overflow to infinity."""
     assert (APyFloat(0, 30, 1, 5, 2) / APyFloat(0, 1, 3, 5, 2)).is_inf
 
+    assert (APyFloat(0, 30, 1, 5, 2) / APyFloat(0, 0, 3, 5, 2)).is_inf
+
 
 @pytest.mark.float_div
 def test_div_underflow():
-    """Test that a multiplication can underflow to zero."""
+    """Test that a division can underflow to zero."""
     assert (APyFloat(0, 1, 3, 5, 2) / APyFloat(1, 30, 3, 5, 2)) == 0
 
+    assert (APyFloat(0, 0, 3, 5, 2) / APyFloat(1, 30, 3, 5, 2)) == 0
+
 
 @pytest.mark.float_div
 @pytest.mark.parametrize("x", [0, 1, float("inf"), float("nan")])
diff --git a/lib/test/apyfloatarray/test_arithmetic.py b/lib/test/apyfloatarray/test_arithmetic.py
index 3d33121ad..e2b1ee2a5 100644
--- a/lib/test/apyfloatarray/test_arithmetic.py
+++ b/lib/test/apyfloatarray/test_arithmetic.py
@@ -392,6 +392,32 @@ def test_array_mul():
     res = APyFloatArray([0], [8], [0], 4, 3) * APyFloatArray([0], [0], [2], 4, 3)
     assert res.is_identical(APyFloatArray([0], [0], [4], 4, 3))
 
+    # 4.0 x subn
+    res = APyFloatArray([0], [9], [0], 4, 3) * APyFloatArray([0], [0], [1], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [4], 4, 3))
+
+    # subnormal becoming normal
+    res = APyFloatArray([0], [10], [0], 4, 3) * APyFloatArray([0], [0], [1], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [1], [0], 4, 3))
+
+    # 0.5 x subn
+    res = APyFloatArray([0], [6], [0], 4, 3) * APyFloatArray([0], [0], [4], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [2], 4, 3))
+
+    # 0.25 x subn
+    res = APyFloatArray([0], [5], [0], 4, 3) * APyFloatArray([0], [0], [4], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [1], 4, 3))
+
+    # 0.125 x subn
+    res = APyFloatArray([0], [4], [0], 4, 3) * APyFloatArray([0], [0], [4], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [0], 4, 3))
+
+    # Two normal numbers that generate carry, but with subnormal result
+    a = APyFloat(sign=1, exp=2, man=957, exp_bits=5, man_bits=10)
+    b = APyFloatArray([0], [12], [1015], 5, 10)
+    res = a * b
+    assert res.is_identical(APyFloatArray([1], [0], [986], 5, 10))
+
 
 @pytest.mark.float_array
 def test_array_mul_scalar():
@@ -411,8 +437,8 @@ def test_array_mul_scalar():
 
     a = APyFloatArray([0, 1], [0, 0], [5, 14], 4, 5)
     b = APyFloat(0, 0, 11, 4, 5)
-    assert (a * b).is_identical(APyFloatArray([0, 1], [0, 0], [0, 1], 4, 5))
-    assert (b * a).is_identical(APyFloatArray([0, 1], [0, 0], [0, 1], 4, 5))
+    assert (a * b).is_identical(APyFloatArray([0, 1], [0, 0], [0, 0], 4, 5))
+    assert (b * a).is_identical(APyFloatArray([0, 1], [0, 0], [0, 0], 4, 5))
 
     a = APyFloatArray.from_float([0.5, 1.5, 5], 10, 35)
     b = APyFloat.from_float(8, 10, 35)
@@ -468,7 +494,7 @@ def test_array_mul_scalar():
 
     a = APyFloatArray([0, 1], [4, 0], [0, 31], 4, 5)
     b = APyFloat(0, 0, 5, 4, 5)
-    ans = APyFloatArray([0, 1], [0, 0], [1, 1], 4, 5)
+    ans = APyFloatArray([0, 1], [0, 0], [1, 0], 4, 5)
     assert (a * b).is_identical(ans)
     assert (b * a).is_identical(ans)
 
@@ -490,6 +516,46 @@ def test_array_mul_scalar():
     assert (a * b).is_identical(ans)
     assert (b * a).is_identical(ans)
 
+    # 1.0 x subn
+    res = APyFloat(0, 7, 0, 4, 3) * APyFloatArray([0], [0], [1], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [1], 4, 3))
+
+    res = APyFloat(0, 7, 0, 4, 3) * APyFloatArray([0], [0], [2], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [2], 4, 3))
+
+    # 2.0 x subn
+    res = APyFloat(0, 8, 0, 4, 3) * APyFloatArray([0], [0], [1], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [2], 4, 3))
+
+    res = APyFloat(0, 8, 0, 4, 3) * APyFloatArray([0], [0], [2], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [4], 4, 3))
+
+    # 4.0 x subn
+    res = APyFloat(0, 9, 0, 4, 3) * APyFloatArray([0], [0], [1], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [4], 4, 3))
+
+    # subnormal becoming normal
+    res = APyFloat(0, 10, 0, 4, 3) * APyFloatArray([0], [0], [1], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [1], [0], 4, 3))
+
+    # 0.5 x subn
+    res = APyFloat(0, 6, 0, 4, 3) * APyFloatArray([0], [0], [4], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [2], 4, 3))
+
+    # 0.25 x subn
+    res = APyFloat(0, 5, 0, 4, 3) * APyFloatArray([0], [0], [4], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [1], 4, 3))
+
+    # 0.125 x subn
+    res = APyFloat(0, 4, 0, 4, 3) * APyFloatArray([0], [0], [4], 4, 3)
+    assert res.is_identical(APyFloatArray([0], [0], [0], 4, 3))
+
+    # Two normal numbers that generate carry, but with subnormal result
+    a = APyFloatArray([1], [2], [957], 5, 10)
+    b = APyFloatArray([0], [12], [1015], 5, 10)
+    res = a * b
+    assert res.is_identical(APyFloatArray([1], [0], [986], 5, 10))
+
 
 @pytest.mark.float_array
 def test_array_mul_int_float():