Repository: apytypes/apytypes (language: python)

Issue:
NumPy mathematical functions
There is a list at https://numpy.org/doc/stable/reference/routines.math.html

- ~Trigonometric functions~
- ~Hyperbolic functions~
- [ ] Rounding (should be straightforward in many cases)
   - [ ] `round`, `around`, `rint`, `fix`, `floor`, `ceil`, `trunc` (should probably return NumPy arrays? Skip `digits`for (a)round? Or just an alias for cast for fixed-point?)
- [ ] Sums, products, differences (most should be possible)
   - [x] `prod` #476
   - [x] `sum` #473
   - [x] `nanprod` #476
   - [x] `nansum` #473
   - [x] `cumprod` #476 
   - [x] `cumsum` #473 
   - [x] `nancumprod` #476 
   - [x] `nancumsum` #473
   - [ ] `diff`
   - [ ] `ediff1d`
   - [ ] `gradient` ? Intermediate word length?
   - [ ] `cross`
   - [ ] `trapz` ? Intermediate word length?
- ~Exponents and logarithms~
- ~Other special functions~
- [ ] Floating point routines (most)
   - [ ] `signbit` (return NumPy array)
   - [ ] `copysign` 
   - [ ] `frexp` (return NumPy arrays of int32 and int64?)
   - [ ] ~`ldexp`~
   - [ ] `nextafter` (use `next_up()` and `next_down()`)
   - [ ] `spacing`?
- [ ] Rational routines (should be possible?)
   - [ ] `lcm`
   - [ ] `gcd` 
- [ ] Arithmetic operations (most can be implemented using +, -, etc)
   - [ ] `add`
   - [ ] `reciprocal` (Word length for fixed-point?)
   - [ ] `positive`
   - [ ] `negative`
   - [ ] `multiply`
   - [ ] `divide`
   - [ ] `power`
   - [ ] `subtract`
   - [ ] `true_divide`
   - [ ] `floor_divide`
   - [ ] ~`float_power`~
   - [ ] `fmod` ?
   - [ ] `mod`
   - [ ] `modf` (Return format?)
   - [ ] `remainder`
   - [ ] `divmod`
- Handling complex numbers (depends on #392)
   - [ ] `angle` Return NumPy float array?
   - [ ] `real`
   - [ ] `imag`
   - [ ] `conj` / `conjugate`
- [ ] Extrema Finding (most)
   - [ ] `maximum`
   - [x] `max`/`amax` #479
   - [ ] `fmax`
   - [x] `nanmax` #479
   - [ ] `minimum`
   - [x] `min`/`amin` #479
   - [ ] `fmin`
   - [x] `nanmin` #479
- [ ] Miscellaneous (some)
   - [x] `convolve` #462 
   - [ ] `clip`
   - [ ] ~`sqrt`~ ?
   - [ ] ~`qbrt`~ ?
   - [ ] `square`
   - [ ] `absolute` / `abs` (word length for complex numbers?, sqrt needed then)
   - [ ] `fabs`
   - [ ] `sign` (return NumPy array)
   - [ ] `heaviside`
   - [ ] `nan_to_num`
   - [ ] ~`real_if_close`~ ?
   - [ ] ~`interp`~ ? (Word length?)

At least a trivial version should be possible for many, like `square` => `*this * (*this)` etc.

One may even consider if it is possible to have a common Array parent class for these and some other ones?

I think that there maybe should be both `apytypes.square(a)` and `a.square()`, where possibly the former is just a Python wrapper for the latter.


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index be21dfd16..2925e85a5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -26,6 +26,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   - `reshape()`
   - `flatten() / ravel()`
   - `convolve()`
+  - `sum()`
+  - `cumsum()`
+  - `nansum()`
+  - `nancumsum()`
 - Added `array_utils.h` to contain shared functions between float and fixed arrays.
 
 ### Changed
diff --git a/lib/apytypes/_apytypes.pyi b/lib/apytypes/_apytypes.pyi
index 0b53489b8..c95224a90 100644
--- a/lib/apytypes/_apytypes.pyi
+++ b/lib/apytypes/_apytypes.pyi
@@ -808,7 +808,133 @@ class APyFixedArray:
         """
 
     def squeeze(self, axis: int | tuple | None = None) -> APyFixedArray:
-        """Squeeze the matrix, removing dimensions of size 1."""
+        """
+        Removes axes of size one at the specified axis/axes, if no axis is given removes all dimensions with size one.
+
+        Parameters
+        ----------
+        axis : tuple or int, optional
+            The axes to squeeze, a given axis with a size other than one will result in an error. No given axes  will be remove all dimensions of size one.
+
+        Returns
+        -------
+        :class:`APyFixedArray`
+
+        Raises
+        -------
+        :class:ValueError
+            If given an axis of a size other than one a ValueError will be thrown.
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+        """
+
+    def sum(self, axis: tuple | int | None = None) -> APyFixedArray | APyFixed:
+        """
+        Returns the sum of the elements along specified axis/axes treating Not a Number as 0.
+
+        Parameters
+        ----------
+        axis : tuple or int, optional
+            The axis/axes to summate across. Will summate the whole array if no int or tuple is specified.
+
+        Returns
+        -------
+        :class:`APyFixedArray` or :class:`APyFixed`
+
+        Raises
+        -------
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+
+        Examples
+        -------
+
+        >>> from apytypes import APyFixedArray
+
+        >>> a = APyFixedArray([1, 2, 3, 4, 5, 6], int_bits=10, frac_bits=0)
+
+        >>> a.sum()
+        APyFixed(21, bits=13, int_bits=13)
+
+        -------
+        """
+
+    def cumsum(self, axis: int | None = None) -> APyFixedArray:
+        """
+        Returns the cumulative sum of the elements along a given axes.
+
+        Parameters
+        ----------
+        axis : int, optional
+            The axes to summate across. If not given an axis it will return the cumulative sum of the flattened array.
+
+        Returns
+        -------
+        :class:`APyFixedArray`
+
+        Raises
+        -------
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+
+        Examples
+        -------
+
+        >>> from apytypes import APyFixedArray
+
+        Array `a`, array to summate across.
+
+        >>> a = APyFixedArray([[1, 2, 3], [4, 5, 6]], int_bits=10, frac_bits=0)
+
+        >>> a.cumsum()
+        APyFixedArray([1, 3, 6, 10, 15, 21], shape=(6,), bits=13, int_bits=13)
+
+        >>> a.cumsum(0)
+        APyFixedArray([1, 2, 3, 5, 7, 9], shape=(2, 3), bits=11, int_bits=11)
+
+        >>> a.cumsum(1)
+        APyFixedArray([1, 3, 6, 4, 9, 15], shape=(2, 3), bits=12, int_bits=12)
+
+        -------
+        """
+
+    def nansum(self, axis: tuple | int | None = None) -> APyFixedArray | APyFixed:
+        """
+        Returns the sum of the elements along specified axis/axes treating Not a Number as 0.
+
+        Parameters
+        ----------
+        axis : tuple or int, optional
+            The axis/axes to summate across. Will summate the whole array if no int or tuple is specified.
+
+        Returns
+        -------
+        :class:`APyFixedArray` or :class:`APyFixed`
+
+        Raises
+        -------
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+        """
+
+    def nancumsum(self, axis: int | None = None) -> APyFixedArray:
+        """
+        Returns the cumulative sum of the elements along a given axis treating Not a Number as 0.
+
+        Parameters
+        ----------
+        axis : int, optional
+            The axis to summate across. If not given an axis it will return the cumulative sum of the flattened array.
+
+        Returns
+        -------
+        :class:`APyFixedArray`
+
+        Raises
+        -------
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+        """
 
     @staticmethod
     def from_float(
@@ -2061,7 +2187,173 @@ class APyFloatArray:
         """
 
     def squeeze(self, axis: int | tuple | None = None) -> APyFloatArray:
-        """Squeeze the matrix, removing dimensions of size 1."""
+        """
+        Remove axes of size one at the specified axis/axes, if no axís is given removes all dimensions with size one.
+
+        Parameters
+        ----------
+        axis : tuple or int, optional
+            The axis/axes to squeeze, a given axis with a size other than one will result in an error. No given axes  will be remove all dimensions of size 1.
+
+        Returns
+        -------
+        :class:`APyFloatArray`
+
+        Raises
+        -------
+        :class:ValueError
+            If given an axis of a size other than one, a ValueError will be thrown.
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+        """
+
+    def sum(self, axis: tuple | int | None = None) -> APyFloatArray | APyFloat:
+        """
+        Returns the sum of the elements along specified axis/axes.
+
+        Parameters
+        ----------
+        axis : tuple or int, optional
+            The axis/axes to summate across. Will summate the whole array if no int or tuple is specified.
+
+        Returns
+        -------
+        :class:`APyFloatArray` or :class:`APyFloat`
+
+        Raises
+        -------
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+
+        Examples
+        -------
+
+        >>> from apytypes import APyFloatArray
+
+        >>> a = APyFloatArray.from_float([1, 2, 3, 4, 5, 6], exp_bits=10, man_bits=10)
+
+        >>> a.sum()
+        APyFloat(sign=0, exp=515, man=320, exp_bits=10, man_bits=10)
+
+        -------
+        """
+
+    def cumsum(self, axis: int | None = None) -> APyFloatArray:
+        """
+        Returns the cumulative sum of the elements along a given axis.
+
+        Parameters
+        ----------
+        axis : int, optional
+            The axis to summate across. If not given an axis it will return the cumulative sum of the flattened array.
+
+        Returns
+        -------
+        :class:`APyFloatArray`
+
+        Raises
+        -------
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+
+        Examples
+        -------
+
+        >>> from apytypes import APyFloatArray
+
+        >>> a = APyFloatArray.from_float(
+        ...     [[1, 2, 3], [4, 5, 6]], exp_bits=10, man_bits=10
+        ... )
+
+        >>> a.cumsum()
+        APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 513, 514, 514, 515], [0, 512, 512, 256, 896, 320], shape=(6,), exp_bits=10, man_bits=10, bias=511)
+
+        >>> a.cumsum(0)
+        APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 512, 513, 513, 514], [0, 0, 512, 256, 768, 128], shape=(2, 3), exp_bits=10, man_bits=10, bias=511)
+
+        >>> a.cumsum(1)
+        APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 513, 513, 514, 514], [0, 512, 512, 0, 128, 896], shape=(2, 3), exp_bits=10, man_bits=10, bias=511)
+
+        -------
+        """
+
+    def nansum(self, axis: tuple | int | None = None) -> APyFloatArray | APyFloat:
+        """
+        Returns the sum of the elements along specified axis/axes treating Not a Number as 0.
+
+        Parameters
+        ----------
+        axis : tuple or int, optional
+            The axis/axes to summate across. Will summate the whole array if no int or tuple is specified.
+
+        Returns
+        -------
+        :class:`APyFloatArray` or :class:`APyFloat`
+
+        Raises
+        -------
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+
+        Examples
+        -------
+
+        >>> from apytypes import APyFloatArray
+
+        >>> nan = float("nan")
+
+        >>> a = APyFloatArray.from_float([1, 2, 3, 4, 5, nan], exp_bits=10, man_bits=10)
+
+        >>> a.nansum()
+        APyFloat(sign=0, exp=514, man=896, exp_bits=10, man_bits=10)
+
+
+
+        -------
+        """
+
+    def nancumsum(self, axis: int | None = None) -> APyFloatArray:
+        """
+        Returns the cumulative sum of the elements along a given axis treating Not a Number as 0.
+
+        Parameters
+        ----------
+        axis : int, optional
+            The axis to summate across. If not given an axis it will return the cumulative sum of the flattened array.
+
+        Returns
+        -------
+        :class:`APyFloatArray`
+
+        Raises
+        -------
+        :class:IndexError
+            If a specified axis is outside of the existing number of dimensions for the array.
+
+        Examples
+        -------
+
+
+        >>> from apytypes import APyFloatArray
+
+        >>> nan = float("nan")
+
+        >>> a = APyFloatArray.from_float(
+        ...     [[1, 2, 3], [4, 5, 6]], exp_bits=10, man_bits=10
+        ... )
+
+        >>> a.nancumsum()
+        APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 513, 514, 514, 515], [0, 512, 512, 256, 896, 320], shape=(6,), exp_bits=10, man_bits=10, bias=511)
+
+        >>> a.nancumsum(0)
+        APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 512, 513, 513, 514], [0, 0, 512, 256, 768, 128], shape=(2, 3), exp_bits=10, man_bits=10, bias=511)
+
+        >>> a.nancumsum(1)
+        APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 513, 513, 514, 514], [0, 512, 512, 0, 128, 896], shape=(2, 3), exp_bits=10, man_bits=10, bias=511)
+
+
+        -------
+        """
 
     def __getitem__(self, idx: int) -> APyFloatArray | APyFloat: ...
     def __iter__(self) -> APyFloatArrayIterator: ...
diff --git a/src/apyfixedarray.cc b/src/apyfixedarray.cc
index 2fcbaa90c..cb89f88ea 100644
--- a/src/apyfixedarray.cc
+++ b/src/apyfixedarray.cc
@@ -804,7 +804,7 @@ APyFixedArray::squeeze(std::optional<std::variant<nb::int_, nb::tuple>> axis) co
             int axis_n = int(nanobind::cast<nb::int_>(*ptr));
             if (axis_n >= int(shape.size())) {
                 throw nb::index_error(
-                    "specified axis with larger than number of dimensions in the "
+                    "Specified axis with larger than number of dimensions in the "
                     "APyFixedArray"
                 );
             }
@@ -814,7 +814,7 @@ APyFixedArray::squeeze(std::optional<std::variant<nb::int_, nb::tuple>> axis) co
         auto predicate = [&](std::size_t dim) {
             if (axis_set.find(cnt) != axis_set.end() && dim != 1) {
                 throw nb::value_error(
-                    "cannot select an axis to squeeze that has size other than one"
+                    "Cannot select an axis to squeeze that has size other than one"
                 );
             }
             return axis_set.find(cnt++) != axis_set.end() && dim == 1;
@@ -836,6 +836,191 @@ APyFixedArray::squeeze(std::optional<std::variant<nb::int_, nb::tuple>> axis) co
     return result;
 }
 
+APyFixedArray APyFixedArray::cumsum(std::optional<nb::int_> axis) const
+{
+    auto pos_func = [](std::size_t i,
+                       std::size_t sec_length,
+                       std::size_t elements,
+                       APyFixedArray& src,
+                       APyFixedArray& dst) {
+        std::size_t pos = i - sec_length;
+        if (i % (sec_length * elements) < sec_length) {
+            pos = i;
+        }
+        mpn_add_n_functor<> {}(
+            &dst._data.at(i * src._itemsize),   // dst
+            &dst._data.at(pos * src._itemsize), // src1
+            &src._data.at(i * src._itemsize),   // src2
+            src._itemsize                       // limb vector length
+        );
+    };
+    return cumulative_prod_sum_function(pos_func, axis);
+}
+
+APyFixedArray APyFixedArray::nancumsum(std::optional<nb::int_> axis) const
+{
+    return this->cumsum(axis);
+}
+
+APyFixedArray APyFixedArray::cumulative_prod_sum_function(
+    void (*pos_func)(std::size_t, std::size_t, std::size_t, APyFixedArray&, APyFixedArray&),
+    std::optional<nb::int_> axis
+) const
+{
+    std::size_t _axis = _shape.size();
+    if (axis.has_value()) {
+        _axis = std::size_t(axis.value());
+        if (_axis >= _shape.size()) {
+            throw nb::index_error(
+                "specified axis outside number of dimensions in the APyFixedArray"
+            );
+        }
+    }
+
+    std::size_t elements = _nitems;
+    std::size_t res_nitems = _nitems;
+    std::vector<std::size_t> shape = { _nitems };
+    if (_axis < _shape.size()) {
+        elements = _shape.at(_axis);
+        shape = _shape;
+    }
+
+    int bit_increase = bit_width(elements - 1);
+    const int res_int_bits = int_bits() + bit_increase;
+    const int res_frac_bits = frac_bits();
+    const int res_bits = res_int_bits + res_frac_bits;
+    std::size_t sec_length = _nitems;
+    for (std::size_t i = 0; i < _shape.size(); i++) {
+        sec_length /= _shape.at(i);
+        if (i == _axis) {
+            break;
+        }
+    }
+    // Resulting vector
+    APyFixedArray result(_shape, res_bits, res_int_bits);
+    APyFixedArray source(_shape, bits(), int_bits());
+    std::copy_n(_data.begin(), _data.size(), source._data.begin());
+
+    for (std::size_t i = 0; i < result._nitems; i++) {
+        pos_func(i, sec_length, elements, source, result);
+    }
+
+    APyFixedArray result2(shape, res_bits, res_int_bits);
+    std::copy_n(
+        result._data.begin(), res_nitems * result._itemsize, result2._data.begin()
+    );
+    return result2;
+}
+std::variant<APyFixedArray, APyFixed> APyFixedArray::prod_sum_function(
+    void (*pos_func)(std::size_t, std::size_t, std::size_t, APyFixedArray&, APyFixedArray&),
+    std::optional<std::variant<nb::tuple, nb::int_>> axis
+) const
+{
+    std::set<std::size_t> axes_set;
+    if (axis.has_value()) {
+        std::vector<std::size_t> axes_vector
+            = cpp_shape_from_python_shape_like(axis.value());
+        for (auto i : axes_vector) {
+            if (i >= _shape.size()) {
+                throw nb::index_error(
+                    "specified axis outside number of dimensions in the APyFixedArray"
+                );
+            }
+            axes_set.insert(i);
+        }
+    } else {
+        axes_set.insert(_shape.size());
+    }
+    std::size_t elements = 1;
+    for (std::size_t i = 0; i <= _shape.size(); i++) {
+        if (axes_set.find(i) != axes_set.end()) {
+            if (i == _shape.size()) {
+                elements *= _nitems;
+            } else {
+                elements *= _shape.at(i);
+            }
+        }
+    }
+    int bit_increase = bit_width(elements - 1);
+    const int res_int_bits = int_bits() + bit_increase;
+    const int res_frac_bits = frac_bits();
+    const int res_bits = res_int_bits + res_frac_bits;
+
+    // Resulting vector
+    APyFixedArray result(_shape, res_bits, res_int_bits);
+    APyFixedArray source(_shape, bits(), int_bits());
+    std::copy_n(_data.begin(), _data.size(), source._data.begin());
+
+    std::vector<std::size_t> shape;
+    std::size_t res_elements = source._nitems;
+    std::size_t sec_length = source._nitems;
+    for (std::size_t x = 0; x <= _shape.size(); x++) {
+        std::size_t elements = res_elements;
+        if (axes_set.find(x) == axes_set.end()) {
+            if (x < _shape.size()) {
+                sec_length /= source._shape.at(x);
+                shape.push_back(_shape.at(x));
+            }
+            continue;
+        }
+
+        if (x < _shape.size()) {
+            elements = source._shape.at(x);
+            sec_length /= source._shape.at(x);
+        } else {
+            shape = {};
+        }
+        for (std::size_t i = 0; i < res_elements; i++) {
+            pos_func(i, sec_length, elements, source, result);
+        }
+
+        res_elements /= elements;
+        source = result;
+        result._data.assign(result._data.size(), 0);
+    }
+    if (shape.size() == 0) {
+        APyFixed res(res_bits, res_int_bits);
+        std::copy_n(source._data.begin(), source._itemsize, res._data.begin());
+        return res;
+    } else {
+        APyFixedArray res(shape, res_bits, res_int_bits);
+        std::copy_n(
+            source._data.begin(), res_elements * source._itemsize, res._data.begin()
+        );
+        return res;
+    }
+}
+
+std::variant<APyFixedArray, APyFixed>
+APyFixedArray::sum(std::optional<std::variant<nb::tuple, nb::int_>> axis) const
+{
+    auto pos_func = [](std::size_t i,
+                       std::size_t sec_length,
+                       std::size_t elements,
+                       APyFixedArray& src,
+                       APyFixedArray& dst) {
+        std::size_t pos_in_sec = i % (sec_length);
+        std::size_t sec_pos
+            = (i - i % (elements * sec_length)) / (elements * sec_length);
+        // Perform ripple-carry operation on the limbs
+        auto pos = (pos_in_sec + sec_pos) * src._itemsize;
+        mpn_add_n_functor<> {}(
+            &dst._data.at(pos),               // dst
+            &dst._data.at(pos),               // src1
+            &src._data.at(i * src._itemsize), // src2
+            src._itemsize                     // limb vector length
+        );
+    };
+
+    return prod_sum_function(pos_func, axis);
+}
+
+std::variant<APyFixedArray, APyFixed>
+APyFixedArray::nansum(std::optional<std::variant<nb::tuple, nb::int_>> axis) const
+{
+    return this->sum(axis);
+}
+
 APyFixedArray
 APyFixedArray::broadcast_to_python(const std::variant<nb::tuple, nb::int_> shape) const
 {
diff --git a/src/apyfixedarray.h b/src/apyfixedarray.h
index aa224b6ee..2c2c4e7e1 100644
--- a/src/apyfixedarray.h
+++ b/src/apyfixedarray.h
@@ -18,6 +18,7 @@ namespace nb = nanobind;
 #include <cstddef>  // std::size_t
 #include <limits>   // std::numeric_limits<>::is_iec559
 #include <optional> // std::optional, std::nullopt
+#include <set>      // std::set
 #include <string>   // std::string
 #include <vector>   // std::vector
 
@@ -102,6 +103,16 @@ class APyFixedArray : public APyBuffer<mp_limb_t> {
     template <class ripple_carry_op, class simd_op_const, class simd_shift_op_const>
     inline APyFixedArray _apyfixed_base_add_sub(const APyFixed& rhs) const;
 
+    std::variant<APyFixedArray, APyFixed> prod_sum_function(
+        void (*pos_func)(std::size_t, std::size_t, std::size_t, APyFixedArray&, APyFixedArray&),
+        std::optional<std::variant<nb::tuple, nb::int_>> axis = std::nullopt
+    ) const;
+
+    APyFixedArray cumulative_prod_sum_function(
+        void (*pos_func)(std::size_t, std::size_t, std::size_t, APyFixedArray&, APyFixedArray&),
+        std::optional<nb::int_> axis = std::nullopt
+    ) const;
+
 public:
     APyFixedArray operator+(const APyFixedArray& rhs) const;
     APyFixedArray operator+(const APyFixed& rhs) const;
@@ -148,6 +159,22 @@ class APyFixedArray : public APyBuffer<mp_limb_t> {
     //! Perform a linear convolution with `other` using `mode`
     APyFixedArray convolve(const APyFixedArray& other, const std::string& mode) const;
 
+    //! Returns a copy where the specified axes is summated.
+    std::variant<APyFixedArray, APyFixed>
+    sum(std::optional<std::variant<nb::tuple, nb::int_>> axis = std::nullopt) const;
+
+    //! Returns a copy where the specified axes contains the increasing cumulated
+    //! summation across its own axis.
+    APyFixedArray cumsum(std::optional<nb::int_> axis = std::nullopt) const;
+
+    //! Returns a copy where the specified axes is summated, treating Nan as 0.
+    std::variant<APyFixedArray, APyFixed>
+    nansum(std::optional<std::variant<nb::tuple, nb::int_>> axis = std::nullopt) const;
+
+    //! Returns a copy where the specified axes contains the increasing cumulated
+    //! summation across its own axis, treating Nan as 0.
+    APyFixedArray nancumsum(std::optional<nb::int_> axis = std::nullopt) const;
+
     //! Python `__repr__()` function
     std::string repr() const;
 
diff --git a/src/apyfixedarray_wrapper.cc b/src/apyfixedarray_wrapper.cc
index 34fa2515e..5fa3e89a9 100644
--- a/src/apyfixedarray_wrapper.cc
+++ b/src/apyfixedarray_wrapper.cc
@@ -479,7 +479,158 @@ void bind_fixed_array(nb::module_& m)
             &APyFixedArray::squeeze,
             nb::arg("axis") = nb::none(),
             R"pbdoc(
-            Squeeze the matrix, removing dimensions of size 1.
+            Removes axes of size one at the specified axis/axes, if no axis is given removes all dimensions with size one.
+
+            Parameters
+            ----------
+            axis : tuple or int, optional
+                The axes to squeeze, a given axis with a size other than one will result in an error. No given axes  will be remove all dimensions of size one.
+
+            Returns
+            -------
+            :class:`APyFixedArray`
+
+            Raises
+            -------
+            :class:ValueError
+                If given an axis of a size other than one a ValueError will be thrown.
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
+            )pbdoc"
+        )
+
+        .def(
+            "sum",
+            &APyFixedArray::sum,
+            nb::arg("axis") = nb::none(),
+            R"pbdoc(
+            Returns the sum of the elements along specified axis/axes treating Not a Number as 0.
+
+            Parameters
+            ----------
+            axis : tuple or int, optional
+                The axis/axes to summate across. Will summate the whole array if no int or tuple is specified.
+
+            Returns
+            -------
+            :class:`APyFixedArray` or :class:`APyFixed`
+
+            Raises
+            -------
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
+            Examples
+            -------
+
+            >>> from apytypes import APyFixedArray
+
+            >>> a = APyFixedArray(
+            ...     [1,2,3,4,5,6],
+            ...     int_bits=10,
+            ...     frac_bits=0
+            ... )
+
+            >>> a.sum()
+            APyFixed(21, bits=13, int_bits=13)
+
+            -------
+            )pbdoc"
+        )
+
+        .def(
+            "cumsum",
+            &APyFixedArray::cumsum,
+            nb::arg("axis") = nb::none(),
+            R"pbdoc(
+            Returns the cumulative sum of the elements along a given axes.
+
+            Parameters
+            ----------
+            axis : int, optional
+                The axes to summate across. If not given an axis it will return the cumulative sum of the flattened array.
+
+            Returns
+            -------
+            :class:`APyFixedArray`
+
+            Raises
+            -------
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
+            Examples
+            -------
+
+            >>> from apytypes import APyFixedArray
+
+            Array `a`, array to summate across.
+
+            >>> a = APyFixedArray(
+            ...     [[1,2,3],[4,5,6]],
+            ...     int_bits=10,
+            ...     frac_bits=0
+            ... )
+
+            >>> a.cumsum()
+            APyFixedArray([1, 3, 6, 10, 15, 21], shape=(6,), bits=13, int_bits=13)
+
+            >>> a.cumsum(0)
+            APyFixedArray([1, 2, 3, 5, 7, 9], shape=(2, 3), bits=11, int_bits=11)
+
+            >>> a.cumsum(1)
+            APyFixedArray([1, 3, 6, 4, 9, 15], shape=(2, 3), bits=12, int_bits=12)
+
+            -------
+            )pbdoc"
+        )
+
+        .def(
+            "nansum",
+            &APyFixedArray::nansum,
+            nb::arg("axis") = nb::none(),
+            R"pbdoc(
+            Returns the sum of the elements along specified axis/axes treating Not a Number as 0.
+
+            Parameters
+            ----------
+            axis : tuple or int, optional
+                The axis/axes to summate across. Will summate the whole array if no int or tuple is specified.
+
+            Returns
+            -------
+            :class:`APyFixedArray` or :class:`APyFixed`
+
+            Raises
+            -------
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
+            )pbdoc"
+        )
+
+        .def(
+            "nancumsum",
+            &APyFixedArray::nancumsum,
+            nb::arg("axis") = nb::none(),
+            R"pbdoc(
+            Returns the cumulative sum of the elements along a given axis treating Not a Number as 0.
+
+            Parameters
+            ----------
+            axis : int, optional
+                The axis to summate across. If not given an axis it will return the cumulative sum of the flattened array.
+
+            Returns
+            -------
+            :class:`APyFixedArray`
+
+            Raises
+            -------
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
             )pbdoc"
         )
 
diff --git a/src/apyfloatarray.cc b/src/apyfloatarray.cc
index a99dab71c..ba59958e1 100644
--- a/src/apyfloatarray.cc
+++ b/src/apyfloatarray.cc
@@ -1007,7 +1007,7 @@ APyFloatArray::squeeze(std::optional<std::variant<nb::int_, nb::tuple>> axis) co
             int axis_n = int(nanobind::cast<nb::int_>(*ptr));
             if (axis_n >= int(_shape.size())) {
                 throw nb::index_error(
-                    "specified axis with larger than number of dimensions in the "
+                    "Specified axis with larger than number of dimensions in the "
                     "APyFloatArray"
                 );
             }
@@ -1018,7 +1018,7 @@ APyFloatArray::squeeze(std::optional<std::variant<nb::int_, nb::tuple>> axis) co
         auto predicate = [&](std::size_t dim) {
             if (axis_set.find(cnt) != axis_set.end() && dim != 1) {
                 throw nb::value_error(
-                    "cannot select an axis to squeeze that has size other than one"
+                    "Cannot select an axis to squeeze that has size other than one"
                 );
             }
             return axis_set.find(cnt++) != axis_set.end() && dim == 1;
@@ -1123,6 +1123,231 @@ APyFloatArray::convolve(const APyFloatArray& other, const std::string& mode) con
     return result;
 }
 
+std::variant<APyFloatArray, APyFloat> APyFloatArray::prod_sum_function(
+    void (*pos_func)(std::size_t, std::size_t, std::size_t, APyFloatArray&, APyFloatArray&, APyFloat&, APyFloat&),
+    std::optional<std::variant<nb::tuple, nb::int_>> axis
+) const
+{
+    std::set<std::size_t> axes_set;
+    if (axis.has_value()) {
+        std::vector<std::size_t> axes_vector
+            = cpp_shape_from_python_shape_like(axis.value());
+        for (auto i : axes_vector) {
+            if (i >= shape.size()) {
+                throw nb::index_error(
+                    "specified axis outside number of dimensions in the APyFixedArray"
+                );
+            }
+            axes_set.insert(i);
+        }
+    } else {
+        axes_set.insert(shape.size());
+    }
+
+    // Resulting vector
+    APyFloatArray result(shape, exp_bits, man_bits);
+    APyFloatArray source(shape, exp_bits, man_bits);
+    std::copy_n(data.begin(), data.size(), source.data.begin());
+
+    APyFloat lhs_scalar(exp_bits, man_bits, bias);
+    APyFloat rhs_scalar(exp_bits, man_bits, bias);
+
+    std::vector<std::size_t> shape_;
+    std::size_t res_elements = source.data.size();
+    std::size_t sec_length = source.data.size();
+    for (std::size_t x = 0; x <= shape.size(); x++) {
+        std::size_t elements = res_elements;
+        if (axes_set.find(x) == axes_set.end()) {
+            if (x < shape.size()) {
+                sec_length /= source.shape.at(x);
+                shape_.push_back(shape.at(x));
+            }
+            continue;
+        }
+
+        if (x < shape.size()) {
+            elements = source.shape.at(x);
+            sec_length /= source.shape.at(x);
+        } else {
+            shape_ = {};
+        }
+        for (std::size_t i = 0; i < res_elements; i++) {
+            pos_func(i, sec_length, elements, source, result, lhs_scalar, rhs_scalar);
+        }
+
+        res_elements /= elements;
+        source = result;
+        result.data.assign(result.data.size(), { 0, 0, 0 });
+    }
+    if (shape_.size() == 0) {
+        APyFloat res(exp_bits, man_bits);
+        res.set_data(source.data.at(0));
+        return res;
+    } else {
+        APyFloatArray res(shape_, exp_bits, man_bits);
+        std::copy_n(source.data.begin(), res_elements, res.data.begin());
+        return res;
+    }
+}
+
+APyFloatArray APyFloatArray::cumulative_prod_sum_function(
+    void (*pos_func)(std::size_t, std::size_t, std::size_t, APyFloatArray&, APyFloatArray&, APyFloat&, APyFloat&),
+    std::optional<nb::int_> axis
+) const
+{
+    std::size_t _axis = shape.size();
+    if (axis.has_value()) {
+        _axis = std::size_t(axis.value());
+        if (_axis >= shape.size()) {
+            throw nb::index_error(
+                "specified axis outside number of dimensions in the APyFloatArray"
+            );
+        }
+    }
+
+    std::size_t elements = data.size();
+    std::size_t res_nitems = data.size();
+    std::vector<std::size_t> shape_ = { data.size() };
+    if (_axis < shape.size()) {
+        elements = shape.at(_axis);
+        shape_ = shape;
+    }
+
+    std::size_t sec_length = data.size();
+    for (std::size_t i = 0; i < shape.size(); i++) {
+        sec_length /= shape.at(i);
+        if (i == _axis) {
+            break;
+        }
+    }
+    // Resulting vector
+    APyFloatArray result(shape, exp_bits, man_bits);
+    APyFloatArray source(shape, exp_bits, man_bits);
+    std::copy_n(data.begin(), data.size(), source.data.begin());
+
+    APyFloat lhs_scalar(exp_bits, man_bits, bias);
+    APyFloat rhs_scalar(exp_bits, man_bits, bias);
+
+    for (std::size_t i = 0; i < result.data.size(); i++) {
+        pos_func(i, sec_length, elements, source, result, lhs_scalar, rhs_scalar);
+    }
+
+    APyFloatArray result2(shape_, exp_bits, man_bits);
+    std::copy_n(result.data.begin(), res_nitems, result2.data.begin());
+    return result2;
+}
+
+std::variant<APyFloatArray, APyFloat>
+APyFloatArray::sum(std::optional<std::variant<nb::tuple, nb::int_>> axis) const
+{
+    auto pos_func = [](std::size_t i,
+                       std::size_t sec_length,
+                       std::size_t elements,
+                       APyFloatArray& src,
+                       APyFloatArray& dst,
+                       APyFloat& lhs_scalar,
+                       APyFloat& rhs_scalar) {
+        // calculate new position
+        std::size_t pos_in_sec = i % (sec_length);
+        std::size_t sec_pos
+            = (i - i % (elements * sec_length)) / (elements * sec_length);
+        auto pos = (pos_in_sec + sec_pos);
+
+        // perform addition
+        lhs_scalar.set_data(dst.data.at(pos));
+        rhs_scalar.set_data(src.data.at(i));
+        dst.data[pos] = (lhs_scalar + rhs_scalar).get_data();
+    };
+
+    return prod_sum_function(pos_func, axis);
+}
+
+APyFloatArray APyFloatArray::cumsum(std::optional<nb::int_> axis) const
+{
+    auto pos_func = [](std::size_t i,
+                       std::size_t sec_length,
+                       std::size_t elements,
+                       APyFloatArray& src,
+                       APyFloatArray& dst,
+                       APyFloat& lhs_scalar,
+                       APyFloat& rhs_scalar) {
+        // calculate new position
+        std::size_t pos = i - sec_length;
+        if (i % (sec_length * elements) < sec_length) {
+            pos = i;
+        }
+
+        // perform addition
+        lhs_scalar.set_data(dst.data.at(pos));
+        rhs_scalar.set_data(src.data.at(i));
+        dst.data[i] = (lhs_scalar + rhs_scalar).get_data();
+    };
+
+    return cumulative_prod_sum_function(pos_func, axis);
+}
+
+std::variant<APyFloatArray, APyFloat>
+APyFloatArray::nansum(std::optional<std::variant<nb::tuple, nb::int_>> axis) const
+{
+    auto pos_func = [](std::size_t i,
+                       std::size_t sec_length,
+                       std::size_t elements,
+                       APyFloatArray& src,
+                       APyFloatArray& dst,
+                       APyFloat& lhs_scalar,
+                       APyFloat& rhs_scalar) {
+        // calculate new position
+        std::size_t pos_in_sec = i % (sec_length);
+        std::size_t sec_pos
+            = (i - i % (elements * sec_length)) / (elements * sec_length);
+        auto pos = (pos_in_sec + sec_pos);
+
+        // perform addition
+        lhs_scalar.set_data(dst.data.at(pos));
+        rhs_scalar.set_data(src.data.at(i));
+
+        if (lhs_scalar.is_nan()) {
+            lhs_scalar.set_data({ 0, 0, 0 });
+        }
+        if (rhs_scalar.is_nan()) {
+            rhs_scalar.set_data({ 0, 0, 0 });
+        }
+        dst.data[pos] = (lhs_scalar + rhs_scalar).get_data();
+    };
+
+    return prod_sum_function(pos_func, axis);
+}
+
+APyFloatArray APyFloatArray::nancumsum(std::optional<nb::int_> axis) const
+{
+    auto pos_func = [](std::size_t i,
+                       std::size_t sec_length,
+                       std::size_t elements,
+                       APyFloatArray& src,
+                       APyFloatArray& dst,
+                       APyFloat& lhs_scalar,
+                       APyFloat& rhs_scalar) {
+        // calculate new position
+        std::size_t pos = i - sec_length;
+        if (i % (sec_length * elements) < sec_length) {
+            pos = i;
+        }
+
+        // perform addition
+        lhs_scalar.set_data(dst.data.at(pos));
+        rhs_scalar.set_data(src.data.at(i));
+        if (lhs_scalar.is_nan()) {
+            lhs_scalar.set_data({ 0, 0, 0 });
+        }
+        if (rhs_scalar.is_nan()) {
+            rhs_scalar.set_data({ 0, 0, 0 });
+        }
+        dst.data[i] = (lhs_scalar + rhs_scalar).get_data();
+    };
+
+    return cumulative_prod_sum_function(pos_func, axis);
+}
+
 std::string APyFloatArray::repr() const
 {
     std::stringstream ss {};
diff --git a/src/apyfloatarray.h b/src/apyfloatarray.h
index 5afbaa5b3..00602ed98 100644
--- a/src/apyfloatarray.h
+++ b/src/apyfloatarray.h
@@ -111,12 +111,29 @@ class APyFloatArray {
     //! Broadcast to a new shape
     APyFloatArray broadcast_to(const std::vector<std::size_t> shape) const;
 
+    //! Returns a copy where all axes of size 1 is removed.
     APyFloatArray
     squeeze(std::optional<std::variant<nb::int_, nb::tuple>> axis = std::nullopt) const;
 
     //! Perform a linear convolution with `other` using `mode`
     APyFloatArray convolve(const APyFloatArray& other, const std::string& mode) const;
 
+    //! Returns a copy where the specified axes is summated.
+    std::variant<APyFloatArray, APyFloat>
+    sum(std::optional<std::variant<nb::tuple, nb::int_>> axis = std::nullopt) const;
+
+    //! Returns a copy where the specified axes contains the increasing cumulated
+    //! summation across its own axis.
+    APyFloatArray cumsum(std::optional<nb::int_> axis = std::nullopt) const;
+
+    //! Returns a copy where the specified axes is summated, treating Nan as 0.
+    std::variant<APyFloatArray, APyFloat>
+    nansum(std::optional<std::variant<nb::tuple, nb::int_>> axis = std::nullopt) const;
+
+    //! Returns a copy where the specified axes contains the increasing cumulated
+    //! summation across its own axis, treating Nan as 0.
+    APyFloatArray nancumsum(std::optional<nb::int_> axis = std::nullopt) const;
+
     //! Python-exposed `broadcast_to`
     APyFloatArray broadcast_to_python(const std::variant<nb::tuple, nb::int_> shape
     ) const;
@@ -289,6 +306,16 @@ class APyFloatArray {
         const APYFLOAT_TYPE& x, // Floating point src1
         const APYFLOAT_TYPE& y  // Floating point src2
     );
+
+    std::variant<APyFloatArray, APyFloat> prod_sum_function(
+        void (*pos_func)(std::size_t, std::size_t, std::size_t, APyFloatArray&, APyFloatArray&, APyFloat&, APyFloat&),
+        std::optional<std::variant<nb::tuple, nb::int_>> axis = std::nullopt
+    ) const;
+
+    APyFloatArray cumulative_prod_sum_function(
+        void (*pos_func)(std::size_t, std::size_t, std::size_t, APyFloatArray&, APyFloatArray&, APyFloat&, APyFloat&),
+        std::optional<nb::int_> axis = std::nullopt
+    ) const;
 };
 
 #endif
diff --git a/src/apyfloatarray_wrapper.cc b/src/apyfloatarray_wrapper.cc
index 593e03734..496f6c8a0 100644
--- a/src/apyfloatarray_wrapper.cc
+++ b/src/apyfloatarray_wrapper.cc
@@ -606,7 +606,200 @@ void bind_float_array(nb::module_& m)
             &APyFloatArray::squeeze,
             nb::arg("axis") = nb::none(),
             R"pbdoc(
-            Squeeze the matrix, removing dimensions of size 1.
+            Remove axes of size one at the specified axis/axes, if no axís is given removes all dimensions with size one.
+
+            Parameters
+            ----------
+            axis : tuple or int, optional
+                The axis/axes to squeeze, a given axis with a size other than one will result in an error. No given axes  will be remove all dimensions of size 1.
+
+            Returns
+            -------
+            :class:`APyFloatArray`
+
+            Raises
+            -------
+            :class:ValueError
+                If given an axis of a size other than one, a ValueError will be thrown.
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
+            )pbdoc"
+        )
+
+        .def(
+            "sum",
+            &APyFloatArray::sum,
+            nb::arg("axis") = nb::none(),
+            R"pbdoc(
+            Returns the sum of the elements along specified axis/axes.
+
+            Parameters
+            ----------
+            axis : tuple or int, optional
+                The axis/axes to summate across. Will summate the whole array if no int or tuple is specified.
+
+            Returns
+            -------
+            :class:`APyFloatArray` or :class:`APyFloat`
+
+            Raises
+            -------
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
+            Examples
+            -------
+
+            >>> from apytypes import APyFloatArray
+
+            >>> a = APyFloatArray.from_float(
+            ...     [1,2,3,4,5,6],
+            ...     exp_bits=10,
+            ...     man_bits=10
+            ... )
+
+            >>> a.sum()
+            APyFloat(sign=0, exp=515, man=320, exp_bits=10, man_bits=10)
+
+            -------
+            )pbdoc"
+        )
+
+        .def(
+            "cumsum",
+            &APyFloatArray::cumsum,
+            nb::arg("axis") = nb::none(),
+            R"pbdoc(
+            Returns the cumulative sum of the elements along a given axis.
+
+            Parameters
+            ----------
+            axis : int, optional
+                The axis to summate across. If not given an axis it will return the cumulative sum of the flattened array.
+
+            Returns
+            -------
+            :class:`APyFloatArray`
+
+            Raises
+            -------
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
+            Examples
+            -------
+
+            >>> from apytypes import APyFloatArray
+
+            >>> a = APyFloatArray.from_float(
+            ...     [[1,2,3],[4,5,6]],
+            ...     exp_bits=10,
+            ...     man_bits=10
+            ... )
+
+            >>> a.cumsum()
+            APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 513, 514, 514, 515], [0, 512, 512, 256, 896, 320], shape=(6,), exp_bits=10, man_bits=10, bias=511)
+
+            >>> a.cumsum(0)
+            APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 512, 513, 513, 514], [0, 0, 512, 256, 768, 128], shape=(2, 3), exp_bits=10, man_bits=10, bias=511)
+
+            >>> a.cumsum(1)
+            APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 513, 513, 514, 514], [0, 512, 512, 0, 128, 896], shape=(2, 3), exp_bits=10, man_bits=10, bias=511)
+
+            -------
+            )pbdoc"
+        )
+
+        .def(
+            "nansum",
+            &APyFloatArray::nansum,
+            nb::arg("axis") = nb::none(),
+            R"pbdoc(
+            Returns the sum of the elements along specified axis/axes treating Not a Number as 0.
+
+            Parameters
+            ----------
+            axis : tuple or int, optional
+                The axis/axes to summate across. Will summate the whole array if no int or tuple is specified.
+
+            Returns
+            -------
+            :class:`APyFloatArray` or :class:`APyFloat`
+
+            Raises
+            -------
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
+            Examples
+            -------
+
+            >>> from apytypes import APyFloatArray
+
+            >>> nan = float("nan")
+
+            >>> a = APyFloatArray.from_float(
+            ...     [1,2,3,4,5,nan],
+            ...     exp_bits=10,
+            ...     man_bits=10
+            ... )
+
+            >>> a.nansum()
+            APyFloat(sign=0, exp=514, man=896, exp_bits=10, man_bits=10)
+
+
+
+            -------
+            )pbdoc"
+        )
+
+        .def(
+            "nancumsum",
+            &APyFloatArray::nancumsum,
+            nb::arg("axis") = nb::none(),
+            R"pbdoc(
+            Returns the cumulative sum of the elements along a given axis treating Not a Number as 0.
+
+            Parameters
+            ----------
+            axis : int, optional
+                The axis to summate across. If not given an axis it will return the cumulative sum of the flattened array.
+
+            Returns
+            -------
+            :class:`APyFloatArray`
+
+            Raises
+            -------
+            :class:IndexError
+                If a specified axis is outside of the existing number of dimensions for the array.
+
+            Examples
+            -------
+
+
+            >>> from apytypes import APyFloatArray
+
+            >>> nan = float("nan")
+
+            >>> a = APyFloatArray.from_float(
+            ...     [[1,2,3],[4,5,6]],
+            ...     exp_bits=10,
+            ...     man_bits=10
+            ... )
+
+            >>> a.nancumsum()
+            APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 513, 514, 514, 515], [0, 512, 512, 256, 896, 320], shape=(6,), exp_bits=10, man_bits=10, bias=511)
+
+            >>> a.nancumsum(0)
+            APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 512, 513, 513, 514], [0, 0, 512, 256, 768, 128], shape=(2, 3), exp_bits=10, man_bits=10, bias=511)
+
+            >>> a.nancumsum(1)
+            APyFloatArray([0, 0, 0, 0, 0, 0], [511, 512, 513, 513, 514, 514], [0, 512, 512, 0, 128, 896], shape=(2, 3), exp_bits=10, man_bits=10, bias=511)
+
+
+            -------
             )pbdoc"
         )
 
diff --git a/src/apytypes_util.h b/src/apytypes_util.h
index 030c1d1d0..c5ddb4bbd 100644
--- a/src/apytypes_util.h
+++ b/src/apytypes_util.h
@@ -110,8 +110,7 @@ template <typename INT_TYPE>
     return leading_zeros(~n);
 }
 
-//! Compute bit-width (`1 + ceil(log2(n))` for unsigned `n`) of an integer `n`. If `n`
-//! is zero, return zero.
+//! Compute bit-width (`ceil(log2(1 + n))` for unsigned `n`) of an integer `n`.
 template <typename INT_TYPE>
 [[maybe_unused, nodiscard]] static APY_INLINE std::size_t bit_width(INT_TYPE n)
 {


Test patch:
diff --git a/lib/test/apyfixedarray/test_methods.py b/lib/test/apyfixedarray/test_methods.py
index 386d89d7e..b36717832 100644
--- a/lib/test/apyfixedarray/test_methods.py
+++ b/lib/test/apyfixedarray/test_methods.py
@@ -40,11 +40,17 @@ def test_squeeze():
     assert k.is_identical(i)
     assert z.is_identical(i)
     m = APyFixedArray([[1], [2], [3]], bits=2, int_bits=2)
-    with pytest.raises(ValueError):
+    with pytest.raises(
+        ValueError,
+        match="Cannot select an axis to squeeze that has size other than one",
+    ):
         _ = m.squeeze(axis=0)
     m1 = m.squeeze(axis=1)
     assert m1.is_identical(APyFixedArray([1, 2, 3], bits=2, int_bits=2))
-    with pytest.raises(IndexError):
+    with pytest.raises(
+        IndexError,
+        match="Specified axis with larger than number of dimensions in the APyFixedArray",
+    ):
         _ = m.squeeze(axis=2)
 
     n = APyFixedArray([[[[[[[[2]]]]]]]], bits=2, int_bits=2)
@@ -53,12 +59,120 @@ def test_squeeze():
     p = n.squeeze((0, 1, 3))
     assert p.is_identical(APyFixedArray([[[[[2]]]]], bits=2, int_bits=2))
     q = APyFixedArray([[[1]], [[2]], [[3]], [[4]]], bits=2, int_bits=2)
-    with pytest.raises(ValueError):
+    with pytest.raises(
+        ValueError,
+        match="Cannot select an axis to squeeze that has size other than one",
+    ):
         _ = q.squeeze((0, 1, 2))
-    with pytest.raises(IndexError):
+    with pytest.raises(
+        IndexError,
+        match="Specified axis with larger than number of dimensions in the APyFixedArray",
+    ):
         _ = m.squeeze((1, 4))
 
 
+@pytest.mark.parametrize("sum_func", ["sum", "nansum"])
+def test_sum(sum_func):
+    a = APyFixedArray([[1, 2], [3, 4], [5, 6], [7, 8]], bits=5, int_bits=5)
+    b = a.__getattribute__(sum_func)()
+    assert b.is_identical(APyFixed(36, bits=8, int_bits=8))
+    c = APyFixedArray([[0, 1, 2], [3, 4, 5]], frac_bits=0, int_bits=5)
+    d = c.__getattribute__(sum_func)((0, 1))
+    e = c.__getattribute__(sum_func)(0)
+    f = c.__getattribute__(sum_func)(1)
+    assert d.is_identical(APyFixed(15, bits=8, int_bits=8))
+    assert e.is_identical(APyFixedArray([3, 5, 7], bits=6, int_bits=6))
+    assert f.is_identical(APyFixedArray([3, 12], bits=7, int_bits=7))
+
+    # test for size larger than 32 and 64 when number over multiple limbs
+    g = APyFixedArray([[0, 1, 2], [3, 4, 5]], frac_bits=0, int_bits=33)
+    h = g.__getattribute__(sum_func)(0)
+    assert h.is_identical(APyFixedArray([3, 5, 7], frac_bits=0, int_bits=34))
+    j = APyFixedArray([[0, 1, 2], [3, 4, 5]], frac_bits=0, int_bits=65)
+    k = j.__getattribute__(sum_func)(0)
+    assert k.is_identical(APyFixedArray([3, 5, 7], frac_bits=0, int_bits=66))
+
+    # test some float and negative summation
+    j = APyFixedArray.from_float([0.25, 1.375, 3.25], frac_bits=3, int_bits=5)
+    k = j.__getattribute__(sum_func)()
+    assert k.is_identical(APyFixed.from_float(4.875, bits=10, int_bits=7))
+    m = APyFixedArray.from_float([0.25, 1.25, 3.25], frac_bits=3, int_bits=5)
+    n = m.__getattribute__(sum_func)()
+    assert n.is_identical(APyFixed.from_float(4.75, frac_bits=3, int_bits=7))
+
+    o = APyFixedArray([[-1, -2], [-3, -4]], frac_bits=0, int_bits=5)
+    p = o.__getattribute__(sum_func)(1)
+    assert p.is_identical(APyFixedArray([-3, -7], frac_bits=0, int_bits=6))
+
+    q = APyFixedArray([[-1, -2], [1, 2]], frac_bits=0, int_bits=5)
+    r = q.__getattribute__(sum_func)(0)
+    assert r.is_identical(APyFixedArray([0, 0], frac_bits=0, int_bits=6))
+
+    m = APyFixedArray([1, 2, 3], bits=2, int_bits=2)
+    with pytest.raises(
+        IndexError,
+        match="specified axis outside number of dimensions in the APyFixedArray",
+    ):
+        _ = m.__getattribute__(sum_func)(1)
+
+
+@pytest.mark.parametrize("sum_func", ["cumsum", "nancumsum"])
+def test_cumsum(sum_func):
+    a = APyFixedArray([[1, 2, 3], [4, 5, 6]], frac_bits=0, int_bits=5)
+    b = a.__getattribute__(sum_func)()
+    assert b.is_identical(APyFixedArray([1, 3, 6, 10, 15, 21], frac_bits=0, int_bits=8))
+    c = a.__getattribute__(sum_func)(0)
+    assert c.is_identical(
+        APyFixedArray([[1, 2, 3], [5, 7, 9]], frac_bits=0, int_bits=6)
+    )
+    d = a.__getattribute__(sum_func)(1)
+    assert d.is_identical(
+        APyFixedArray([[1, 3, 6], [4, 9, 15]], frac_bits=0, int_bits=7)
+    )
+    e = APyFixedArray([[[1, 2], [3, 4]], [[5, 6], [7, 8]]], frac_bits=0, int_bits=8)
+    f = e.__getattribute__(sum_func)()
+    g = e.__getattribute__(sum_func)(0)
+    h = e.__getattribute__(sum_func)(1)
+    i = e.__getattribute__(sum_func)(2)
+    assert f.is_identical(
+        APyFixedArray([1, 3, 6, 10, 15, 21, 28, 36], frac_bits=0, int_bits=11)
+    )
+    assert g.is_identical(
+        APyFixedArray([[[1, 2], [3, 4]], [[6, 8], [10, 12]]], frac_bits=0, int_bits=9)
+    )
+    assert h.is_identical(
+        APyFixedArray([[[1, 2], [4, 6]], [[5, 6], [12, 14]]], frac_bits=0, int_bits=9)
+    )
+    assert i.is_identical(
+        APyFixedArray([[[1, 3], [3, 7]], [[5, 11], [7, 15]]], frac_bits=0, int_bits=9)
+    )
+    with pytest.raises(
+        IndexError,
+        match="specified axis outside number of dimensions in the APyFixedArray",
+    ):
+        _ = e.__getattribute__(sum_func)(4)
+
+    k = APyFixedArray.from_float(
+        [[0.25, 0.25], [0.25, 0.25]], frac_bits=10, int_bits=10
+    )
+    m = k.__getattribute__(sum_func)()
+    assert m.is_identical(
+        APyFixedArray.from_float([0.25, 0.5, 0.75, 1], frac_bits=10, int_bits=12)
+    )
+
+    # test for size larger than 32 and 64 when number over multiple limbs
+    g = APyFixedArray([[0, 1, 2], [3, 4, 5]], frac_bits=0, int_bits=33)
+    h = g.__getattribute__(sum_func)(0)
+    assert h.is_identical(
+        APyFixedArray([[0, 1, 2], [3, 5, 7]], frac_bits=0, int_bits=34)
+    )
+    j = APyFixedArray([[0, 1, 2], [3, 4, 5]], frac_bits=0, int_bits=65)
+    k = j.__getattribute__(sum_func)(0)
+    assert k.is_identical(
+        APyFixedArray([[0, 1, 2], [3, 5, 7]], frac_bits=0, int_bits=66)
+    )
+
+
 def test_to_numpy():
     # Skip this test if `NumPy` is not present on the machine
     np = pytest.importorskip("numpy")
@@ -309,4 +423,3 @@ def test_ravel(shape):
     reshaped = arr.reshape(shape)
     if not APyFixedArray.is_identical(reshaped.ravel(), arr):
         pytest.fail(f"ravel didn't return to original 1d list after reshape {shape}")
-    pass
diff --git a/lib/test/apyfloatarray/test_methods.py b/lib/test/apyfloatarray/test_methods.py
index fbad4feb7..66c11258f 100644
--- a/lib/test/apyfloatarray/test_methods.py
+++ b/lib/test/apyfloatarray/test_methods.py
@@ -304,6 +304,195 @@ def test_squeeze():
         _ = m.squeeze((1, 4))
 
 
+def test_sum():
+    a = APyFloatArray.from_float(
+        [[1, 2], [3, 4], [5, 6], [7, 8]], exp_bits=10, man_bits=10
+    )
+    b = a.sum()
+    assert b.is_identical(APyFloat.from_float(36, exp_bits=10, man_bits=10))
+    c = APyFloatArray.from_float([[0, 1, 2], [3, 4, 5]], exp_bits=10, man_bits=10)
+    d = c.sum((0, 1))
+    e = c.sum(0)
+    f = c.sum(1)
+    assert d.is_identical(APyFloat.from_float(15, exp_bits=10, man_bits=10))
+    assert e.is_identical(APyFloatArray.from_float([3, 5, 7], exp_bits=10, man_bits=10))
+    assert f.is_identical(APyFloatArray.from_float([3, 12], exp_bits=10, man_bits=10))
+
+    # test for size larger than 32 and 64 when number over multiple limbs
+    g = APyFloatArray.from_float([[0, 1, 2], [3, 4, 5]], exp_bits=10, man_bits=10)
+    h = g.sum(0)
+    assert h.is_identical(APyFloatArray.from_float([3, 5, 7], exp_bits=10, man_bits=10))
+    j = APyFloatArray.from_float([[0, 1, 2], [3, 4, 5]], exp_bits=10, man_bits=10)
+    k = j.sum(0)
+    assert k.is_identical(APyFloatArray.from_float([3, 5, 7], exp_bits=10, man_bits=10))
+
+    # test some float and negative summation
+    j = APyFloatArray.from_float([0.25, 1.375, 3.25], exp_bits=10, man_bits=10)
+    k = j.sum()
+    assert k.is_identical(APyFloat.from_float(4.875, exp_bits=10, man_bits=10))
+    m = APyFloatArray.from_float([0.25, 1.25, 3.25], exp_bits=10, man_bits=10)
+    n = m.sum()
+    assert n.is_identical(APyFloat.from_float(4.75, exp_bits=10, man_bits=10))
+
+    o = APyFloatArray.from_float([[-1, -2], [-3, -4]], exp_bits=10, man_bits=10)
+    p = o.sum(1)
+    assert p.is_identical(APyFloatArray.from_float([-3, -7], exp_bits=10, man_bits=10))
+
+    q = APyFloatArray.from_float([[-1, -2], [1, 2]], exp_bits=10, man_bits=10)
+    r = q.sum(0)
+    assert r.is_identical(APyFloatArray.from_float([0, 0], exp_bits=10, man_bits=10))
+
+    m = APyFloatArray.from_float([1, 2, 3], exp_bits=10, man_bits=10)
+    with pytest.raises(IndexError):
+        _ = m.sum(1)
+
+
+def test_cumsum():
+    a = APyFloatArray.from_float([[1, 2, 3], [4, 5, 6]], exp_bits=10, man_bits=10)
+    b = a.cumsum()
+    assert b.is_identical(
+        APyFloatArray.from_float([1, 3, 6, 10, 15, 21], exp_bits=10, man_bits=10)
+    )
+    c = a.cumsum(0)
+    assert c.is_identical(
+        APyFloatArray.from_float([[1, 2, 3], [5, 7, 9]], exp_bits=10, man_bits=10)
+    )
+    d = a.cumsum(1)
+    assert d.is_identical(
+        APyFloatArray.from_float([[1, 3, 6], [4, 9, 15]], exp_bits=10, man_bits=10)
+    )
+    e = APyFloatArray.from_float(
+        [[[1, 2], [3, 4]], [[5, 6], [7, 8]]], exp_bits=10, man_bits=10
+    )
+    f = e.cumsum()
+    g = e.cumsum(0)
+    h = e.cumsum(1)
+    i = e.cumsum(2)
+    assert f.is_identical(
+        APyFloatArray.from_float(
+            [1, 3, 6, 10, 15, 21, 28, 36], exp_bits=10, man_bits=10
+        )
+    )
+    assert g.is_identical(
+        APyFloatArray.from_float(
+            [[[1, 2], [3, 4]], [[6, 8], [10, 12]]], exp_bits=10, man_bits=10
+        )
+    )
+    assert h.is_identical(
+        APyFloatArray.from_float(
+            [[[1, 2], [4, 6]], [[5, 6], [12, 14]]], exp_bits=10, man_bits=10
+        )
+    )
+    assert i.is_identical(
+        APyFloatArray.from_float(
+            [[[1, 3], [3, 7]], [[5, 11], [7, 15]]], exp_bits=10, man_bits=10
+        )
+    )
+    with pytest.raises(IndexError):
+        _ = e.cumsum(4)
+
+    k = APyFloatArray.from_float([[0.25, 0.25], [0.25, 0.25]], exp_bits=10, man_bits=10)
+    m = k.cumsum()
+    assert m.is_identical(
+        APyFloatArray.from_float([0.25, 0.5, 0.75, 1], exp_bits=10, man_bits=10)
+    )
+
+
+def test_nansum():
+    nan = float("nan")
+    a = APyFloatArray.from_float(
+        [[nan, 2], [3, 4], [5, 6], [7, 8]], exp_bits=10, man_bits=10
+    )
+    b = a.nansum()
+    assert b.is_identical(APyFloat.from_float(35, exp_bits=10, man_bits=10))
+    c = APyFloatArray.from_float([[0, 1, nan], [3, 4, 5]], exp_bits=10, man_bits=10)
+    d = c.nansum((0, 1))
+    e = c.nansum(0)
+    f = c.nansum(1)
+    assert d.is_identical(APyFloat.from_float(13, exp_bits=10, man_bits=10))
+    assert e.is_identical(APyFloatArray.from_float([3, 5, 5], exp_bits=10, man_bits=10))
+    assert f.is_identical(APyFloatArray.from_float([1, 12], exp_bits=10, man_bits=10))
+
+    g = APyFloatArray.from_float([[0, 1, 2], [nan, 4, 5]], exp_bits=10, man_bits=10)
+    h = g.nansum(0)
+    assert h.is_identical(APyFloatArray.from_float([0, 5, 7], exp_bits=10, man_bits=10))
+    j = APyFloatArray.from_float([[0, 1, 2], [nan, nan, nan]], exp_bits=10, man_bits=10)
+    k = j.nansum(0)
+    assert k.is_identical(APyFloatArray.from_float([0, 1, 2], exp_bits=10, man_bits=10))
+
+    # test some float and negative summation
+    j = APyFloatArray.from_float([0.25, 1.375, 3.25], exp_bits=10, man_bits=10)
+    k = j.sum()
+    assert k.is_identical(APyFloat.from_float(4.875, exp_bits=10, man_bits=10))
+    m = APyFloatArray.from_float([0.25, 1.25, 3.25], exp_bits=10, man_bits=10)
+    n = m.sum()
+    assert n.is_identical(APyFloat.from_float(4.75, exp_bits=10, man_bits=10))
+
+    o = APyFloatArray.from_float([[-1, -2], [-3, nan]], exp_bits=10, man_bits=10)
+    p = o.nansum(1)
+    assert p.is_identical(APyFloatArray.from_float([-3, -3], exp_bits=10, man_bits=10))
+
+    q = APyFloatArray.from_float([[-1, -2], [1, nan]], exp_bits=10, man_bits=10)
+    r = q.nansum(0)
+    assert r.is_identical(APyFloatArray.from_float([0, -2], exp_bits=10, man_bits=10))
+
+    m = APyFloatArray.from_float([1, 2, 3], exp_bits=10, man_bits=10)
+    with pytest.raises(IndexError):
+        _ = m.nansum(1)
+
+
+def test_nancumsum():
+    nan = float("nan")
+    a = APyFloatArray.from_float([[1, 2, 3], [4, 5, nan]], exp_bits=10, man_bits=10)
+    b = a.nancumsum()
+    assert b.is_identical(
+        APyFloatArray.from_float([1, 3, 6, 10, 15, 15], exp_bits=10, man_bits=10)
+    )
+    c = a.nancumsum(0)
+    assert c.is_identical(
+        APyFloatArray.from_float([[1, 2, 3], [5, 7, 3]], exp_bits=10, man_bits=10)
+    )
+    d = a.nancumsum(1)
+    assert d.is_identical(
+        APyFloatArray.from_float([[1, 3, 6], [4, 9, 9]], exp_bits=10, man_bits=10)
+    )
+    e = APyFloatArray.from_float(
+        [[[1, 2], [3, 4]], [[nan, nan], [7, 8]]], exp_bits=10, man_bits=10
+    )
+    f = e.nancumsum()
+    g = e.nancumsum(0)
+    h = e.nancumsum(1)
+    i = e.nancumsum(2)
+    assert f.is_identical(
+        APyFloatArray.from_float(
+            [1, 3, 6, 10, 10, 10, 17, 25], exp_bits=10, man_bits=10
+        )
+    )
+    assert g.is_identical(
+        APyFloatArray.from_float(
+            [[[1, 2], [3, 4]], [[1, 2], [10, 12]]], exp_bits=10, man_bits=10
+        )
+    )
+    assert h.is_identical(
+        APyFloatArray.from_float(
+            [[[1, 2], [4, 6]], [[0, 0], [7, 8]]], exp_bits=10, man_bits=10
+        )
+    )
+    assert i.is_identical(
+        APyFloatArray.from_float(
+            [[[1, 3], [3, 7]], [[0, 0], [7, 15]]], exp_bits=10, man_bits=10
+        )
+    )
+    with pytest.raises(IndexError):
+        _ = e.nancumsum(4)
+
+    k = APyFloatArray.from_float([[0.25, 0.25], [0.25, nan]], exp_bits=10, man_bits=10)
+    m = k.nancumsum()
+    assert m.is_identical(
+        APyFloatArray.from_float([0.25, 0.5, 0.75, 0.75], exp_bits=10, man_bits=10)
+    )
+
+
 def test_convenience_cast():
     a = APyFloatArray.from_float([0.893820, 3e20, -float("inf"), float("nan")], 10, 50)
     assert a.cast_to_double().is_identical(