Repository: apytypes/apytypes (language: python)

Issue:
Broadcasting
https://numpy.org/doc/stable/user/basics.broadcasting.html#general-broadcasting-rules (this part we support).

https://numpy.org/doc/stable/user/basics.broadcasting.html#broadcastable-arrays can be added.


Discussion / hints:
A key step to solve this is to implement `broadcast_to` https://numpy.org/doc/stable/reference/generated/numpy.broadcast_to.html

Once that is done, one can simply broadcast one argument to the other size and then use the standard operations. Although this will often be inefficient (as one will have to allocate a new object and copy stuff), it will at least make things work and one can, as a second step, introduce broadcasting rules "locally" in the arithmetic operations to avoid copying etc.
* [`numpy.broadcast()`](https://numpy.org/doc/stable/reference/generated/numpy.broadcast.html#numpy.broadcast)
* [`numpy.broadcast_to()`](https://numpy.org/doc/stable/reference/generated/numpy.broadcast_to.html)
* [`numpy.broadcast_arrays()`](https://numpy.org/doc/stable/reference/generated/numpy.broadcast_arrays.html)
### NumPy error on faulty broadcast:
```python
>>> import numpy as np
>>> np.broadcast_to(np.zeros([4,3,2]), (6,2,6,4))
ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (4,3,2)  and requested shape (6,2,6,4)
```

Fix patch:
diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md
index 08948c332..7d812155e 100644
--- a/.github/PULL_REQUEST_TEMPLATE.md
+++ b/.github/PULL_REQUEST_TEMPLATE.md
@@ -14,3 +14,4 @@ issues and PRs.
 - [ ] "closes #0000" is in the body of the PR description to [link the related issue](https://docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue)
 - [ ] new and changed code is tested
 - [ ] relevant changes added to [CHANGELOG.md](https://github.com/apytypes/apytypes/blob/main/CHANGELOG.md)
+- [ ] new functionality is documented
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 26ffa4d07..437a22b58 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,13 +11,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Methods `next_up()` and `next_down()` for floating-point,
   based on the IEEE-754 standard.
+- Broadcasting support [as NumPy](https://numpy.org/doc/stable/user/basics.broadcasting.html)
+  - `APyFixedArray.broadcast_to()`
+  - `APyFloatArray.broadcast_to()`
+  - Automatic broadcasting on array arithmetic
 
 ### Changed
 
 ### Fixed
 
 - Smaller fixes in quantization for floating-point scalar multiplication.
-- Fix in floating-point 'cast'-method when result became subnormal.
+- Fix in floating-point `cast`-method when result became subnormal.
+- Fix bug in `APyFixedArray.__truediv__()` where program crashes on long-limb zero
+  denominators.
+- Fix bug where single-dimensional tuple representations missed out on a comma
 
 ### Removed
 
diff --git a/docs/api/apyfixedarray.rst b/docs/api/apyfixedarray.rst
index 3e1e620ee..df5b77842 100644
--- a/docs/api/apyfixedarray.rst
+++ b/docs/api/apyfixedarray.rst
@@ -30,11 +30,16 @@
 
    .. automethod:: to_numpy
 
-   Tranposition
-   ------------
+   Transposition
+   -------------
 
    .. automethod:: transpose
 
+   Broadcasting
+   ------------
+
+   .. automethod:: broadcast_to
+
    Properties
    ----------
 
@@ -54,7 +59,7 @@
 
    .. autoproperty:: shape
 
-   Tranposition
-   ^^^^^^^^^^^^
+   Transposition
+   ^^^^^^^^^^^^^
 
    .. autoproperty:: T
diff --git a/docs/api/apyfloatarray.rst b/docs/api/apyfloatarray.rst
index bfaf42ae6..cf34fa889 100644
--- a/docs/api/apyfloatarray.rst
+++ b/docs/api/apyfloatarray.rst
@@ -35,6 +35,11 @@
 
    .. automethod:: transpose
 
+   Broadcasting
+   ------------
+
+   .. automethod:: broadcast_to
+
    Convenience casting methods
    ---------------------------
 
diff --git a/docs/cpp_api/broadcast.rst b/docs/cpp_api/broadcast.rst
new file mode 100644
index 000000000..228d0ec4c
--- /dev/null
+++ b/docs/cpp_api/broadcast.rst
@@ -0,0 +1,5 @@
+``broadcast.h``
+===============
+
+.. doxygenfile:: broadcast.h
+    :project: APyTypes
diff --git a/docs/cpp_api/index.rst b/docs/cpp_api/index.rst
index a8cc0d30a..d926d465c 100644
--- a/docs/cpp_api/index.rst
+++ b/docs/cpp_api/index.rst
@@ -16,5 +16,6 @@ APyTypes C++ API
    apytypes_scratch_vector
    apytypes_simd
    apytypes_util
+   broadcast
    ieee754
    python_util
diff --git a/lib/apytypes/_apytypes.pyi b/lib/apytypes/_apytypes.pyi
index 6bf9dc7ce..fd8e63f32 100644
--- a/lib/apytypes/_apytypes.pyi
+++ b/lib/apytypes/_apytypes.pyi
@@ -574,6 +574,7 @@ class APyFixedArray:
         :class:`APyFixedArray`
         """
 
+    def broadcast_to(self, shape: tuple | int) -> APyFixedArray: ...
     def is_identical(self, other: APyFixedArray) -> bool:
         """
         Test if two :class:`APyFixedArray` objects are identical.
@@ -1156,6 +1157,7 @@ class APyFloatArray:
             see :func:`get_float_quantization_mode`, is used.
         """
 
+    def broadcast_to(self, shape: tuple | int) -> APyFloatArray: ...
     def cast_to_double(
         self, quantization: QuantizationMode | None = None
     ) -> APyFloatArray:
diff --git a/src/apybuffer.h b/src/apybuffer.h
index c389132a5..c30d97f64 100644
--- a/src/apybuffer.h
+++ b/src/apybuffer.h
@@ -18,28 +18,6 @@
 
 #include "apytypes_util.h"
 
-//! Fold a shape under multiplication
-[[maybe_unused]] static APY_INLINE std::size_t
-fold_shape(const std::vector<std::size_t>& shape)
-{
-    return std::accumulate(shape.cbegin(), shape.cend(), 1, std::multiplies {});
-}
-
-//! Retrieve the byte-strides from a shape
-template <typename T>
-[[maybe_unused]] static APY_INLINE std::vector<std::size_t>
-byte_strides_from_shape(const std::vector<std::size_t>& shape, std::size_t itemsize = 1)
-{
-    std::size_t n_bytes = sizeof(T) * itemsize;
-    std::vector<std::size_t> strides(shape.size(), 0);
-    for (std::size_t i = 0; i < shape.size(); i++) {
-        strides[shape.size() - 1 - i] = std::accumulate(
-            shape.crbegin(), shape.crbegin() + i, n_bytes, std::multiplies {}
-        );
-    }
-    return strides;
-}
-
 template <typename T, typename Allocator = std::allocator<T>> class APyBuffer {
 
     //! APyBuffers are to be inherited from. All fields and constructors are protected.
@@ -58,7 +36,7 @@ template <typename T, typename Allocator = std::allocator<T>> class APyBuffer {
     //! Retrieve a Python Buffer structure compatible with the Buffer Protocol
     Py_buffer get_py_buffer()
     {
-        _strides = byte_strides_from_shape<T>(_shape, _itemsize);
+        _strides = strides_from_shape(_shape, _itemsize * sizeof(T));
         return Py_buffer {
             (void*)&_data[0],                // void       *buf
             nullptr,                         // PyObject   *obj
diff --git a/src/apyfixedarray.cc b/src/apyfixedarray.cc
index cc2cb18c7..6eae28db4 100644
--- a/src/apyfixedarray.cc
+++ b/src/apyfixedarray.cc
@@ -2,6 +2,7 @@
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
 #include <nanobind/stl/variant.h> // std::variant (with nanobind support)
+#include <variant>
 namespace nb = nanobind;
 
 // Python details. These should be included before standard header files:
@@ -27,6 +28,7 @@ namespace nb = nanobind;
 #include "apytypes_common.h"
 #include "apytypes_simd.h"
 #include "apytypes_util.h"
+#include "broadcast.h"
 #include "python_util.h"
 
 #include <fmt/format.h>
@@ -245,13 +247,16 @@ inline APyFixedArray APyFixedArray::_apyfixed_base_add_sub(const APyFixed& rhs)
 
 APyFixedArray APyFixedArray::operator+(const APyFixedArray& rhs) const
 {
-    // Make sure `_shape` of `*this` and `rhs` are the same
     if (_shape != rhs._shape) {
-        throw std::length_error(fmt::format(
-            "APyFixedArray.__add__: shape mismatch, lhs.shape={}, rhs.shape={}",
-            string_from_vec(_shape),
-            string_from_vec(rhs._shape)
-        ));
+        auto broadcast_shape = smallest_broadcastable_shape(_shape, rhs._shape);
+        if (broadcast_shape.size() == 0) {
+            throw std::length_error(fmt::format(
+                "APyFixedArray.__add__: shape mismatch, lhs.shape={}, rhs.shape={}",
+                tuple_string_from_vec(_shape),
+                tuple_string_from_vec(rhs._shape)
+            ));
+        }
+        return broadcast_to(broadcast_shape) + rhs.broadcast_to(broadcast_shape);
     }
 
     return _apyfixedarray_base_add_sub<
@@ -270,13 +275,16 @@ APyFixedArray APyFixedArray::operator+(const APyFixed& rhs) const
 
 APyFixedArray APyFixedArray::operator-(const APyFixedArray& rhs) const
 {
-    // Make sure `_shape` of `*this` and `rhs` are the same
     if (_shape != rhs._shape) {
-        throw std::length_error(fmt::format(
-            "APyFixedArray.__sub__: shape mismatch, lhs.shape={}, rhs.shape={}",
-            string_from_vec(_shape),
-            string_from_vec(rhs._shape)
-        ));
+        auto broadcast_shape = smallest_broadcastable_shape(_shape, rhs._shape);
+        if (broadcast_shape.size() == 0) {
+            throw std::length_error(fmt::format(
+                "APyFixedArray.__sub__: shape mismatch, lhs.shape={}, rhs.shape={}",
+                tuple_string_from_vec(_shape),
+                tuple_string_from_vec(rhs._shape)
+            ));
+        }
+        return broadcast_to(broadcast_shape) - rhs.broadcast_to(broadcast_shape);
     }
 
     return _apyfixedarray_base_add_sub<
@@ -331,13 +339,16 @@ APyFixedArray APyFixedArray::rsub(const APyFixed& lhs) const
 
 APyFixedArray APyFixedArray::operator*(const APyFixedArray& rhs) const
 {
-    // Make sure `_shape` of `*this` and `rhs` are the same
     if (_shape != rhs._shape) {
-        throw std::length_error(fmt::format(
-            "APyFixedArray.__mul__: shape mismatch, lhs.shape=({}), rhs.shape=({})",
-            string_from_vec(_shape),
-            string_from_vec(rhs._shape)
-        ));
+        auto broadcast_shape = smallest_broadcastable_shape(_shape, rhs._shape);
+        if (broadcast_shape.size() == 0) {
+            throw std::length_error(fmt::format(
+                "APyFixedArray.__mul__: shape mismatch, lhs.shape={}, rhs.shape={}",
+                tuple_string_from_vec(_shape),
+                tuple_string_from_vec(rhs._shape)
+            ));
+        }
+        return broadcast_to(broadcast_shape) * rhs.broadcast_to(broadcast_shape);
     }
 
     const int res_int_bits = int_bits() + rhs.int_bits();
@@ -449,11 +460,15 @@ APyFixedArray APyFixedArray::operator/(const APyFixedArray& rhs) const
 {
     // Make sure `_shape` of `*this` and `rhs` are the same
     if (_shape != rhs._shape) {
-        throw std::length_error(fmt::format(
-            "APyFixedArray.__div__: shape mismatch, lhs.shape=({}), rhs.shape=({})",
-            string_from_vec(_shape),
-            string_from_vec(rhs._shape)
-        ));
+        auto broadcast_shape = smallest_broadcastable_shape(_shape, rhs._shape);
+        if (broadcast_shape.size() == 0) {
+            throw std::length_error(fmt::format(
+                "APyFixedArray.__div__: shape mismatch, lhs.shape={}, rhs.shape={}",
+                tuple_string_from_vec(_shape),
+                tuple_string_from_vec(rhs._shape)
+            ));
+        }
+        return broadcast_to(broadcast_shape) / rhs.broadcast_to(broadcast_shape);
     }
 
     const int res_int_bits = int_bits() + rhs.frac_bits() + 1;
@@ -483,6 +498,12 @@ APyFixedArray APyFixedArray::operator/(const APyFixedArray& rhs) const
 
     for (std::size_t i = 0; i < _nitems; i++) {
         std::fill(std::begin(abs_num), std::end(abs_num), 0);
+        if (limb_vector_is_zero(
+                std::begin(rhs._data) + (i + 0) * rhs._itemsize,
+                std::begin(rhs._data) + (i + 1) * rhs._itemsize
+            )) {
+            continue;
+        }
         bool den_sign = limb_vector_abs(
             std::begin(rhs._data) + (i + 0) * rhs._itemsize,
             std::begin(rhs._data) + (i + 1) * rhs._itemsize,
@@ -699,9 +720,9 @@ APyFixedArray APyFixedArray::matmul(const APyFixedArray& rhs) const
 
     // Unsupported `__matmul__` dimensionality, raise exception
     throw std::length_error(fmt::format(
-        "APyFixedArray.__matmul__: input shape mismatch, lhs: ({}), rhs: ({})",
-        string_from_vec(_shape),
-        string_from_vec(rhs._shape)
+        "APyFixedArray.__matmul__: input shape mismatch, lhs: {}, rhs: {}",
+        tuple_string_from_vec(_shape),
+        tuple_string_from_vec(rhs._shape)
     ));
 }
 
@@ -739,6 +760,30 @@ APyFixedArray APyFixedArray::transpose() const
  * *                               Other methods                                    * *
  * ********************************************************************************** */
 
+APyFixedArray APyFixedArray::broadcast_to(const std::vector<std::size_t> shape) const
+{
+    if (!is_broadcastable(_shape, shape)) {
+        throw nb::value_error(
+            fmt::format(
+                "Operands could not be broadcast together with shapes: {}, {}",
+                tuple_string_from_vec(_shape),
+                tuple_string_from_vec(shape)
+            )
+                .c_str()
+        );
+    }
+
+    APyFixedArray result(shape, bits(), int_bits());
+    broadcast_data_copy(_data.begin(), result._data.begin(), _shape, shape, _itemsize);
+    return result;
+}
+
+APyFixedArray
+APyFixedArray::broadcast_to_python(const std::variant<nb::tuple, nb::int_> shape) const
+{
+    return broadcast_to(cpp_shape_from_python_shape_like(shape));
+}
+
 std::string APyFixedArray::repr() const
 {
     std::stringstream ss {};
@@ -764,9 +809,9 @@ std::string APyFixedArray::repr() const
 
         ss.seekp(-2, ss.cur);
     }
-    ss << "], shape=(";
-    ss << string_from_vec(_shape);
-    ss << "), "
+    ss << "], shape=";
+    ss << tuple_string_from_vec(_shape);
+    ss << ", "
        << "bits=" << std::dec << bits() << ", "
        << "int_bits=" << std::dec << int_bits() << ")";
     return ss.str();
diff --git a/src/apyfixedarray.h b/src/apyfixedarray.h
index e8acb1550..71aff7846 100644
--- a/src/apyfixedarray.h
+++ b/src/apyfixedarray.h
@@ -134,6 +134,13 @@ class APyFixedArray : public APyBuffer<mp_limb_t> {
      * *                          Public member functions                           * *
      * ****************************************************************************** */
 
+    //! Broadcast to a new shape
+    APyFixedArray broadcast_to(const std::vector<std::size_t> shape) const;
+
+    //! Python-exposed `broadcast_to`
+    APyFixedArray broadcast_to_python(const std::variant<nb::tuple, nb::int_> shape
+    ) const;
+
     //! Python `__repr__()` function
     std::string repr() const;
 
diff --git a/src/apyfixedarray_wrapper.cc b/src/apyfixedarray_wrapper.cc
index 38719ad25..d2be6f265 100644
--- a/src/apyfixedarray_wrapper.cc
+++ b/src/apyfixedarray_wrapper.cc
@@ -1,7 +1,6 @@
 #include "apyfixed.h"
 #include "apyfixedarray.h"
 #include "apyfixedarray_iterator.h"
-#include "apytypes_common.h"
 
 #include <nanobind/nanobind.h>
 #include <nanobind/operators.h>
@@ -343,6 +342,24 @@ void bind_fixed_array(nb::module_& m)
             :class:`APyFixedArray`
             )pbdoc"
         )
+        .def(
+            "broadcast_to",
+            &APyFixedArray::broadcast_to_python,
+            nb::arg("shape"),
+            R"pbdoc(
+            Broadcast array to new shape.
+
+            Parameters
+            ----------
+            shape : tuple or int
+                The shape to broadcast to. A single integer ``i`` is interpreted as ``(i,)``.
+
+            Returns
+            -------
+            :class:`APyFixedArray`
+
+            )pbdoc"
+        )
 
         /*
          * Static methods
@@ -369,11 +386,11 @@ void bind_fixed_array(nb::module_& m)
                 Floating point values to initialize from. The tensor shape will be taken
                 from the sequence shape.
             int_bits : int, optional
-                Number of integer bits in the created fixed-point tensor
+                Number of integer bits in the created fixed-point tensor.
             frac_bits : int, optional
-                Number of fractional bits in the created fixed-point tensor
+                Number of fractional bits in the created fixed-point tensor.
             bits : int, optional
-                Total number of bits in the created fixed-point tensor
+                Total number of bits in the created fixed-point tensor.
 
             Returns
             -------
@@ -421,11 +438,11 @@ void bind_fixed_array(nb::module_& m)
             ndarray : ndarray
                 Values to initialize from. The tensor shape will be taken from the ndarray shape.
             int_bits : int, optional
-                Number of integer bits in the created fixed-point tensor
+                Number of integer bits in the created fixed-point tensor.
             frac_bits : int, optional
-                Number of fractional bits in the created fixed-point tensor
+                Number of fractional bits in the created fixed-point tensor.
             bits : int, optional
-                Total number of bits in the created fixed-point tensor
+                Total number of bits in the created fixed-point tensor.
 
             Returns
             -------
diff --git a/src/apyfloatarray.cc b/src/apyfloatarray.cc
index e6eb46fe7..651be7b21 100644
--- a/src/apyfloatarray.cc
+++ b/src/apyfloatarray.cc
@@ -7,6 +7,7 @@ namespace nb = nanobind;
 #include "apyfloat.h"
 #include "apyfloat_util.h"
 #include "apyfloatarray.h"
+#include "broadcast.h"
 #include "ieee754.h"
 #include "python_util.h"
 #include <algorithm>
@@ -95,13 +96,16 @@ APyFloatArray::APyFloatArray(
 
 APyFloatArray APyFloatArray::operator+(const APyFloatArray& rhs) const
 {
-    // Make sure `_shape` of `*this` and `rhs` are the same
     if (shape != rhs.shape) {
-        throw std::length_error(fmt::format(
-            "APyFloatArray.__add__: shape mismatch, lhs.shape={}, rhs.shape={}",
-            string_from_vec(shape),
-            string_from_vec(rhs.shape)
-        ));
+        auto broadcast_shape = smallest_broadcastable_shape(shape, rhs.shape);
+        if (broadcast_shape.size() == 0) {
+            throw std::length_error(fmt::format(
+                "APyFloatArray.__add__: shape mismatch, lhs.shape={}, rhs.shape={}",
+                tuple_string_from_vec(shape),
+                tuple_string_from_vec(rhs.shape)
+            ));
+        }
+        return broadcast_to(broadcast_shape) + rhs.broadcast_to(broadcast_shape);
     }
 
     const auto quantization = get_float_quantization_mode();
@@ -466,14 +470,18 @@ APyFloatArray APyFloatArray::operator+(const APyFloat& rhs) const
 
 APyFloatArray APyFloatArray::operator-(const APyFloatArray& rhs) const
 {
-    // Make sure `_shape` of `*this` and `rhs` are the same
     if (shape != rhs.shape) {
-        throw std::length_error(fmt::format(
-            "APyFloatArray.__sub__: shape mismatch, lhs.shape={}, rhs.shape={}",
-            string_from_vec(shape),
-            string_from_vec(rhs.shape)
-        ));
+        auto broadcast_shape = smallest_broadcastable_shape(shape, rhs.shape);
+        if (broadcast_shape.size() == 0) {
+            throw std::length_error(fmt::format(
+                "APyFloatArray.__sub__: shape mismatch, lhs.shape={}, rhs.shape={}",
+                tuple_string_from_vec(shape),
+                tuple_string_from_vec(rhs.shape)
+            ));
+        }
+        return broadcast_to(broadcast_shape) - rhs.broadcast_to(broadcast_shape);
     }
+
     return operator+(-rhs);
 }
 
@@ -502,13 +510,16 @@ APyFloatArray APyFloatArray::abs() const
 
 APyFloatArray APyFloatArray::operator*(const APyFloatArray& rhs) const
 {
-    // Make sure `_shape` of `*this` and `rhs` are the same
     if (shape != rhs.shape) {
-        throw std::length_error(fmt::format(
-            "APyFloatArray.__mul__: shape mismatch, lhs.shape={}, rhs.shape={}",
-            string_from_vec(shape),
-            string_from_vec(rhs.shape)
-        ));
+        auto broadcast_shape = smallest_broadcastable_shape(shape, rhs.shape);
+        if (broadcast_shape.size() == 0) {
+            throw std::length_error(fmt::format(
+                "APyFloatArray.__mul__: shape mismatch, lhs.shape={}, rhs.shape={}",
+                tuple_string_from_vec(shape),
+                tuple_string_from_vec(rhs.shape)
+            ));
+        }
+        return broadcast_to(broadcast_shape) * rhs.broadcast_to(broadcast_shape);
     }
 
     // Calculate new format
@@ -874,14 +885,16 @@ APyFloatArray APyFloatArray::operator*(const APyFloat& rhs) const
 
 APyFloatArray APyFloatArray::operator/(const APyFloatArray& rhs) const
 {
-    // Make sure `_shape` of `*this` and `rhs` are the same
     if (shape != rhs.shape) {
-        throw std::length_error(fmt::format(
-            "APyFloatArray.__truediv__: shape mismatch, lhs.shape={}, "
-            "rhs.shape={}",
-            string_from_vec(shape),
-            string_from_vec(rhs.shape)
-        ));
+        auto broadcast_shape = smallest_broadcastable_shape(shape, rhs.shape);
+        if (broadcast_shape.size() == 0) {
+            throw std::length_error(fmt::format(
+                "APyFloatArray.__truediv__: shape mismatch, lhs.shape={}, rhs.shape={}",
+                tuple_string_from_vec(shape),
+                tuple_string_from_vec(rhs.shape)
+            ));
+        }
+        return broadcast_to(broadcast_shape) / rhs.broadcast_to(broadcast_shape);
     }
 
     // Calculate new format
@@ -967,9 +980,9 @@ std::variant<APyFloatArray, APyFloat> APyFloatArray::matmul(const APyFloatArray&
 
     // Unsupported `__matmul__` dimensionality, raise exception
     throw std::length_error(fmt::format(
-        "APyFloatArray.__matmul__: input shape mismatch, lhs: ({}), rhs: ({})",
-        string_from_vec(shape),
-        string_from_vec(rhs.shape)
+        "APyFloatArray.__matmul__: input shape mismatch, lhs: {}, rhs: {}",
+        tuple_string_from_vec(shape),
+        tuple_string_from_vec(rhs.shape)
     ));
 }
 
@@ -994,9 +1007,9 @@ std::string APyFloatArray::repr() const
     } else {
         ss << "[], [], [], ";
     }
-    ss << "shape=(";
-    ss << string_from_vec(shape);
-    ss << "), "
+    ss << "shape=";
+    ss << tuple_string_from_vec(shape);
+    ss << ", "
        << "exp_bits=" << static_cast<unsigned>(exp_bits) << ", "
        << "man_bits=" << static_cast<unsigned>(man_bits) << ", "
        << "bias=" << bias << ")";
@@ -1188,6 +1201,29 @@ void APyFloatArray::_set_values_from_ndarray(const nb::ndarray<nb::c_contig>& nd
     );
 }
 
+APyFloatArray APyFloatArray::broadcast_to(const std::vector<std::size_t> shape) const
+{
+    if (!is_broadcastable(this->shape, shape)) {
+        throw nb::value_error(
+            fmt::format(
+                "Operands could not be broadcast together with shapes: {}, {}",
+                tuple_string_from_vec(this->shape),
+                tuple_string_from_vec(shape)
+            )
+                .c_str()
+        );
+    }
+    APyFloatArray result(shape, exp_bits, man_bits, bias);
+    broadcast_data_copy(data.begin(), result.data.begin(), this->shape, shape);
+    return result;
+}
+
+APyFloatArray
+APyFloatArray::broadcast_to_python(const std::variant<nb::tuple, nb::int_> shape) const
+{
+    return broadcast_to(cpp_shape_from_python_shape_like(shape));
+}
+
 APyFloatArray APyFloatArray::transpose() const
 {
     if (get_ndim() > 2) {
diff --git a/src/apyfloatarray.h b/src/apyfloatarray.h
index 3cf52ddf5..64c6bce6f 100644
--- a/src/apyfloatarray.h
+++ b/src/apyfloatarray.h
@@ -126,6 +126,17 @@ class APyFloatArray {
 
     void _set_values_from_ndarray(const nanobind::ndarray<nanobind::c_contig>& ndarray);
 
+    /* ****************************************************************************** *
+     * *                          Public member functions                           * *
+     * ****************************************************************************** */
+
+    //! Broadcast to a new shape
+    APyFloatArray broadcast_to(const std::vector<std::size_t> shape) const;
+
+    //! Python-exposed `broadcast_to`
+    APyFloatArray broadcast_to_python(const std::variant<nb::tuple, nb::int_> shape
+    ) const;
+
     //! Transposition function. For a 1-D array, returns an exact copy of `*this`. For
     //! a 2-D array, returns the matrix transposition of `*this`.
     APyFloatArray transpose() const;
diff --git a/src/apyfloatarray_wrapper.cc b/src/apyfloatarray_wrapper.cc
index a76b18d51..f5b0e8662 100644
--- a/src/apyfloatarray_wrapper.cc
+++ b/src/apyfloatarray_wrapper.cc
@@ -452,6 +452,24 @@ void bind_float_array(nb::module_& m)
             -------
             :class:`APyFloatArray`
             )pbdoc")
+        .def(
+            "broadcast_to",
+            &APyFloatArray::broadcast_to_python,
+            nb::arg("shape"),
+            R"pbdoc(
+            Broadcast array to new shape.
+
+            Parameters
+            ----------
+            shape : tuple or int
+                The shape to broadcast to. A single integer ``i`` is interpreted as ``(i,)``.
+
+            Returns
+            -------
+            :class:`APyFloatArray`
+
+            )pbdoc"
+        )
 
         // Iteration and friends
         .def("__getitem__", &APyFloatArray::get_item, nb::arg("idx"))
diff --git a/src/apytypes_util.h b/src/apytypes_util.h
index b7da5007c..5685db134 100644
--- a/src/apytypes_util.h
+++ b/src/apytypes_util.h
@@ -8,6 +8,8 @@
 #include <nanobind/nanobind.h>
 #include <nanobind/stl/string.h>
 
+#include "fmt/format.h"
+
 #include <algorithm>        // std::find
 #include <cstddef>          // std::size_t
 #include <cstdint>          // int64_t
@@ -16,11 +18,13 @@
 #include <iomanip>          // std::setfill, std::setw
 #include <ios>              // std::hex
 #include <iterator>         // std::distance
+#include <numeric>          // std::accumulate, std::multiplies
 #include <optional>         // std::optional, std::nullopt
 #include <regex>            // std::regex, std::regex_replace
 #include <sstream>          // std::stringstream
 #include <string>           // std::string
 #include <tuple>            // std::tuple
+#include <variant>          // std::variant
 #include <vector>           // std::vector
 
 /*
@@ -1004,17 +1008,22 @@ uint64_t_from_limb_vector(const std::vector<mp_limb_t>& limb_vec, std::size_t n)
     }
 }
 
-template <typename T> std::string string_from_vec(const std::vector<T>& vec)
+//! Construct a Python tuple-literal-string from a vector of `T`. The type `T` must be
+//! convertible to string through both a `std::stringstream` and `fmt::format({})`.
+template <typename T> std::string tuple_string_from_vec(const std::vector<T>& vec)
 {
     if (vec.size() == 0) {
-        return "";
-    }
-
-    std::stringstream ss;
-    for (auto& d : vec) {
-        ss << d << ", ";
+        return "()";
+    } else if (vec.size() == 1) {
+        return fmt::format("({},)", vec[0]);
+    } else {
+        std::stringstream ss;
+        ss << "(";
+        for (auto& d : vec) {
+            ss << d << ", ";
+        }
+        return ss.str().substr(0, ss.str().length() - 2).append(")");
     }
-    return ss.str().substr(0, ss.str().length() - 2);
 }
 
 [[maybe_unused, nodiscard]] static APY_INLINE mp_limb_t
@@ -1030,6 +1039,43 @@ twos_complement_overflow(mp_limb_t value, int bits)
     return value;
 }
 
+//! Fold a shape under multiplication
+[[maybe_unused]] static APY_INLINE std::size_t
+fold_shape(const std::vector<std::size_t>& shape)
+{
+    return std::accumulate(shape.cbegin(), shape.cend(), 1, std::multiplies {});
+}
+
+//! Compute the strides from a shape
+[[maybe_unused]] static APY_INLINE std::vector<std::size_t>
+strides_from_shape(const std::vector<std::size_t>& shape, std::size_t acc_base = 1)
+{
+    std::vector<std::size_t> strides(shape.size(), 0);
+    for (std::size_t i = 0; i < shape.size(); i++) {
+        strides[shape.size() - 1 - i] = std::accumulate(
+            shape.crbegin(), shape.crbegin() + i, acc_base, std::multiplies {}
+        );
+    }
+    return strides;
+}
+
+//! Create a C++ shape vector (`std::vector<std::size_t>`) from a Python shape object
+//! (`std::variant<nanobind::tuple, nanobind::int_>`).
+static APY_INLINE std::vector<std::size_t> cpp_shape_from_python_shape_like(
+    const std::variant<nanobind::tuple, nanobind::int_>& shape
+)
+{
+    std::vector<std::size_t> cpp_shape {};
+    if (std::holds_alternative<nanobind::tuple>(shape)) {
+        for (const auto& element : std::get<nanobind::tuple>(shape)) {
+            cpp_shape.push_back(nanobind::cast<std::size_t>(element));
+        }
+    } else {
+        cpp_shape.push_back(static_cast<std::size_t>(std::get<nanobind::int_>(shape)));
+    }
+    return cpp_shape;
+}
+
 //! Macro for creating a void-specialization state-less functor `FUNCTOR_NAME` from a
 //! function `FUNC_NAME`. The void-specialization functor allows template argument
 //! deduction to be performed once its function is called. Neat!
diff --git a/src/broadcast.h b/src/broadcast.h
new file mode 100644
index 000000000..a3780feaa
--- /dev/null
+++ b/src/broadcast.h
@@ -0,0 +1,162 @@
+/*
+ * Common general broadcasting functionality for APyTypes
+ *
+ * NumPy fundamental broadcasting:
+ * https://numpy.org/doc/stable/user/basics.broadcasting.html
+ */
+
+#include "apytypes_util.h"
+
+#include <algorithm> // std::any_of, std::copy_n
+#include <tuple>     // std::make_tuple
+#include <vector>    // std::vector
+
+//! Test if `src_shape` can be broadcast to `dst_shape`
+static APY_INLINE bool is_broadcastable(
+    const std::vector<std::size_t>& src_shape, const std::vector<std::size_t>& dst_shape
+)
+{
+    // Can not broadcast if either shape is zero-dimensional
+    if (src_shape.size() == 0 || dst_shape.size() == 0) {
+        return false;
+    }
+
+    // Can not broadcast if either shape has a zero dimension
+    auto is_zero = [](auto n) { return n == 0; };
+    for (const auto& shape : { src_shape, dst_shape }) {
+        if (std::any_of(std::begin(shape), std::end(shape), is_zero)) {
+            return false;
+        }
+    }
+
+    // Can not broadcast if destination shape has fewer dimensions than source shape
+    if (src_shape.size() > dst_shape.size()) {
+        return false;
+    }
+
+    // Iterate shapes from the trailing (right-most) dimensions
+    auto [src_it, dst_it] = std::make_tuple(src_shape.crbegin(), dst_shape.crbegin());
+    while (src_it != src_shape.crend() && dst_it != dst_shape.crend()) {
+        if (*src_it != 1 && *src_it != *dst_it) {
+            return false;
+        }
+        ++src_it;
+        ++dst_it;
+    }
+
+    return true;
+}
+
+//! Get the smallest broadcastable shape from `shape1` and `shape2`. Returns an empty
+//! vector `_.size() == 0` if the shapes can not be broadcast together.
+static APY_INLINE std::vector<std::size_t> smallest_broadcastable_shape(
+    const std::vector<std::size_t>& shape1, const std::vector<std::size_t>& shape2
+)
+{
+    // Can not broadcast if either shape is zero-dimensional
+    if (shape1.size() == 0 || shape2.size() == 0) {
+        return {};
+    }
+
+    // Can not broadcast if any shape dimension is zero
+    auto is_zero = [](auto n) { return n == 0; };
+    for (const auto& shape : { shape1, shape2 }) {
+        if (std::any_of(std::begin(shape), std::end(shape), is_zero)) {
+            return {};
+        }
+    }
+
+    // Iterate shapes from the trailing (right-most) dimensions
+    std::vector<std::size_t> result {};
+    auto [shape1_it, shape2_it] = std::make_tuple(shape1.crbegin(), shape2.crbegin());
+    while (shape1_it != shape1.crend() && shape2_it != shape2.crend()) {
+        if (*shape1_it != 1 && *shape2_it != 1 && *shape1_it != *shape2_it) {
+            // Dimension disallows broadcasting
+            return {};
+        }
+        result.push_back(std::max(*shape1_it, *shape2_it));
+        ++shape1_it;
+        ++shape2_it;
+    }
+
+    // Add remaining shapes from the more high-dimensional shape
+    for (; shape1_it < shape1.crend(); ++shape1_it) {
+        result.push_back(*shape1_it);
+    }
+    for (; shape2_it < shape2.crend(); ++shape2_it) {
+        result.push_back(*shape2_it);
+    }
+
+    // Reverse the dimensions and return
+    std::reverse(std::begin(result), std::end(result));
+    assert(is_broadcastable(shape1, result));
+    assert(is_broadcastable(shape2, result));
+    return result;
+}
+
+//! Compute the destination index of source index `i` based on a `broadcasting_rule`.
+//! The destination index is weighted using `strides`. This function assumes
+//! that `strides.size() == broadcast_rule.size()`
+template <typename FUNC>
+static APY_INLINE std::size_t src_to_dst_idx(
+    std::size_t i, FUNC&& broadcast_rule, const std::vector<std::size_t>& strides
+)
+{
+    std::size_t offset = 0;
+    for (std::size_t j = strides.size(); j--;) {
+        if (i == 0) {
+            break;
+        }
+        offset += strides[j] * (i % broadcast_rule(j));
+        i /= broadcast_rule(j);
+    }
+    return offset;
+}
+
+//! Perform a broadcast by copying data. This function assumes that `src_shape` can be
+//! broadcast to `dst_shape`, i.e., `is_broadcastable(src_shape, dst_shape) == true`. It
+//! further assumes that the data in `src` is stored in C-style order and that `dst` has
+//! enough space to store the broadcast result.
+template <typename RANDOM_ACCESS_CONST_ITERATOR, typename RANDOM_ACCESS_ITERATOR>
+static APY_INLINE void broadcast_data_copy(
+    RANDOM_ACCESS_CONST_ITERATOR src,
+    RANDOM_ACCESS_ITERATOR dst,
+    const std::vector<std::size_t>& src_shape,
+    const std::vector<std::size_t>& dst_shape,
+    std::size_t itemsize = 1
+)
+{
+    std::size_t src_elements = fold_shape(src_shape);
+    std::size_t dst_elements = fold_shape(dst_shape);
+    std::size_t broadcast_elements = dst_elements / src_elements;
+
+    if (src_elements == dst_elements) {
+        std::copy_n(src, src_elements * itemsize, dst);
+        return; // early exit
+    }
+
+    // Destination strides and intermediate woorking coordinate vector
+    std::vector<std::size_t> dst_stride = strides_from_shape(dst_shape);
+
+    // Lambda for left-padded source shape
+    auto src_shape_adjst = [&](std::size_t i) -> std::size_t {
+        std::size_t shape_idx_diff = dst_shape.size() - src_shape.size();
+        return i < shape_idx_diff ? 1 : src_shape[i - shape_idx_diff];
+    };
+
+    // Compute broadcast offsets (stride-weighted permutations of broadcasting vector)
+    auto broadcast_vec = [&](auto j) { return dst_shape[j] - src_shape_adjst(j) + 1; };
+    std::vector<std::size_t> broadcast_offsets(broadcast_elements);
+    for (std::size_t i = 0; i < broadcast_elements; i++) {
+        broadcast_offsets[i] = src_to_dst_idx(i, broadcast_vec, dst_stride);
+    }
+
+    // Loop over elements in the source vector and broadcast to the destination
+    for (std::size_t i = 0; i < src_elements; i++) {
+        auto dst_idx = src_to_dst_idx(i, src_shape_adjst, dst_stride);
+        auto broadcast = [&](auto j) {
+            std::copy_n(src + i * itemsize, itemsize, dst + (dst_idx + j) * itemsize);
+        };
+        std::for_each(broadcast_offsets.begin(), broadcast_offsets.end(), broadcast);
+    }
+}


Test patch:
diff --git a/lib/test/apyfixedarray/test_arithmetic.py b/lib/test/apyfixedarray/test_arithmetic.py
index 8ea302772..1d1c95027 100644
--- a/lib/test/apyfixedarray/test_arithmetic.py
+++ b/lib/test/apyfixedarray/test_arithmetic.py
@@ -360,6 +360,11 @@ def test_array_div():
         )
     )
 
+    # Does not die
+    a = APyFixedArray.from_float([5], bits=100, int_bits=50)
+    b = APyFixedArray.from_float([0], bits=250, int_bits=100)
+    assert (a / b).is_identical(APyFixedArray([0], bits=351, int_bits=201))
+
 
 def test_array_div_scalar():
     a = APyFixedArray.from_float([-5, -6, 7, 8, 9, 5], bits=10, int_bits=5)
diff --git a/lib/test/apyfixedarray/test_matmul.py b/lib/test/apyfixedarray/test_matmul.py
index b26cc43f7..2c053a47b 100644
--- a/lib/test/apyfixedarray/test_matmul.py
+++ b/lib/test/apyfixedarray/test_matmul.py
@@ -17,7 +17,7 @@ def test_dimension_mismatch_raises():
     with pytest.raises(
         ValueError,
         match=re.escape(
-            "APyFixedArray.__matmul__: input shape mismatch, lhs: (3), rhs: (2)"
+            "APyFixedArray.__matmul__: input shape mismatch, lhs: (3,), rhs: (2,)"
         ),
     ):
         _ = a @ b
@@ -25,7 +25,7 @@ def test_dimension_mismatch_raises():
     with pytest.raises(
         ValueError,
         match=re.escape(
-            "APyFixedArray.__matmul__: input shape mismatch, lhs: (2), rhs: (3)"
+            "APyFixedArray.__matmul__: input shape mismatch, lhs: (2,), rhs: (3,)"
         ),
     ):
         _ = b @ a
diff --git a/lib/test/apyfixedarray/test_representation.py b/lib/test/apyfixedarray/test_representation.py
index bcb23cd16..75f7314f3 100644
--- a/lib/test/apyfixedarray/test_representation.py
+++ b/lib/test/apyfixedarray/test_representation.py
@@ -6,11 +6,11 @@
 def test_array_representation():
     assert (
         APyFixedArray([], bits=1, int_bits=-5).__repr__()
-        == "APyFixedArray([], shape=(0), bits=1, int_bits=-5)"
+        == "APyFixedArray([], shape=(0,), bits=1, int_bits=-5)"
     )
     assert (
         APyFixedArray([1, 2, 3], bits=2, int_bits=2).__repr__()
-        == "APyFixedArray([1, 2, 3], shape=(3), bits=2, int_bits=2)"
+        == "APyFixedArray([1, 2, 3], shape=(3,), bits=2, int_bits=2)"
     )
     assert (
         APyFixedArray([range(5, 8), range(3)], bits=10, int_bits=20).__repr__()
diff --git a/lib/test/apyfloatarray/test_arithmetic.py b/lib/test/apyfloatarray/test_arithmetic.py
index f27ab0cba..0c567e7af 100644
--- a/lib/test/apyfloatarray/test_arithmetic.py
+++ b/lib/test/apyfloatarray/test_arithmetic.py
@@ -9,7 +9,7 @@
 
 @pytest.mark.float_array
 def test_array_raises():
-    a = APyFloatArray([1], [5], [4], 10, 10)
+    a = APyFloatArray([1, 2, 3], [5, 2, 3], [4, 2, 3], 10, 10)
     b = APyFloatArray([1, 2], [5, 3], [4, 4], 10, 10)
     with pytest.raises(ValueError, match="APyFloatArray.__add__: shape mismatch"):
         _ = a + b
diff --git a/lib/test/apyfloatarray/test_representation.py b/lib/test/apyfloatarray/test_representation.py
index f2e227ec0..41fb8eac1 100644
--- a/lib/test/apyfloatarray/test_representation.py
+++ b/lib/test/apyfloatarray/test_representation.py
@@ -7,19 +7,19 @@ def test_repr():
     arr = APyFloatArray([], [], [], 6, 7)
     assert (
         repr(arr)
-        == "APyFloatArray([], [], [], shape=(0), exp_bits=6, man_bits=7, bias=31)"
+        == "APyFloatArray([], [], [], shape=(0,), exp_bits=6, man_bits=7, bias=31)"
     )
 
     arr = APyFloatArray([1], [5], [3], 5, 2)
     assert (
         repr(arr)
-        == "APyFloatArray([1], [5], [3], shape=(1), exp_bits=5, man_bits=2, bias=15)"
+        == "APyFloatArray([1], [5], [3], shape=(1,), exp_bits=5, man_bits=2, bias=15)"
     )
 
     arr = APyFloatArray([1, 0], [5, 10], [3, 1], 5, 2)
     assert (
         repr(arr)
-        == "APyFloatArray([1, 0], [5, 10], [3, 1], shape=(2), exp_bits=5, man_bits=2, bias=15)"
+        == "APyFloatArray([1, 0], [5, 10], [3, 1], shape=(2,), exp_bits=5, man_bits=2, bias=15)"
     )
 
     arr = APyFloatArray([[1, 0]], [[5, 10]], [[3, 1]], 5, 2)
diff --git a/lib/test/test_broadcasting.py b/lib/test/test_broadcasting.py
new file mode 100644
index 000000000..dd89db05f
--- /dev/null
+++ b/lib/test/test_broadcasting.py
@@ -0,0 +1,92 @@
+from apytypes import APyFixedArray, APyFloatArray
+
+import pytest
+from typing import List, Tuple
+from itertools import product
+
+
+def all_shapes(ndim: int, dim_elements: int):
+    """
+    Retrieve a list of all possible shapes with `ndim` dimensions, spanning from one to
+    `dim_elements` (inclusive)
+    """
+    return product(range(1, dim_elements + 1), repeat=ndim)
+
+
+# APyTypes array convenience generators
+apytypes_from_float_generators = [
+    lambda f: APyFixedArray.from_float(f, int_bits=200, frac_bits=0),  # multi-limb int
+    lambda f: APyFloatArray.from_float(f, man_bits=52, exp_bits=11),  # double-precision
+]
+
+
+@pytest.mark.parametrize("apyarray_from_float", apytypes_from_float_generators)
+def test_array_broadcast_to(apyarray_from_float):
+    MAX_NDIM = 3
+    MAX_DIM_ELEMENTS = 3
+    np = pytest.importorskip("numpy")
+
+    # Iterate over all possible source and destination shapes and test
+    for src_ndim, dst_ndim in product(range(1, MAX_NDIM + 1), range(1, MAX_NDIM + 1)):
+        src_shapes = all_shapes(src_ndim, MAX_DIM_ELEMENTS)
+        dst_shapes = all_shapes(dst_ndim, MAX_DIM_ELEMENTS)
+        for src_shape, dst_shape in product(src_shapes, dst_shapes):
+            numpy_src = np.arange(np.prod(src_shape), dtype=float).reshape(src_shape)
+            apy_src = apyarray_from_float(numpy_src)
+            try:
+                numpy_res = np.broadcast_to(numpy_src, dst_shape)
+                apy_res = apy_src.broadcast_to(dst_shape).to_numpy()
+                assert np.all(numpy_res == apy_res)
+            except ValueError:
+                # `ValueError` raised by `np.broadcast_to()` or `apy_src.broadcast_to()`
+                # for illegal broadcast. Make sure *both* arrays raises `ValueError`
+                with pytest.raises(ValueError):
+                    np.broadcast_to(numpy_src, dst_shape)
+                with pytest.raises(ValueError, match="Operands could not be broadcast"):
+                    apy_src.broadcast_to(dst_shape)
+
+
+@pytest.mark.parametrize("apyarray_from_float", apytypes_from_float_generators)
+def test_array_broadcast_to_integer(apyarray_from_float):
+    np = pytest.importorskip("numpy")
+    assert np.all(
+        apyarray_from_float([1.0]).broadcast_to(9).to_numpy() == np.ones((9,))
+    )
+
+
+@pytest.mark.parametrize("apyarray_from_float", apytypes_from_float_generators)
+def test_array_broadcast_to_raises(apyarray_from_float):
+    # Can not broadcast if any of the destination dimensions are zero
+    with pytest.raises(ValueError, match="Operands could not be broadcast together"):
+        apyarray_from_float([1.0]).broadcast_to((3, 2, 1, 0, 1))
+
+    # Can not broadcast if the destination number-of-dimensions is zero
+    with pytest.raises(ValueError, match="Operands could not be broadcast together"):
+        apyarray_from_float([1.0]).broadcast_to(tuple())
+
+
+@pytest.mark.parametrize(
+    "a_shape, b_shape",
+    [
+        ((5, 1, 3, 4), (5, 5, 1, 4)),
+        ((1, 3, 4), (5, 5, 1, 4)),
+        ((5, 1, 3, 4), (5, 1, 4)),
+    ],
+)
+@pytest.mark.parametrize(
+    "bin_func",
+    [
+        lambda a, b: a + b,
+        lambda a, b: a - b,
+        lambda a, b: a / b,
+        lambda a, b: a * b,
+    ],
+)
+@pytest.mark.parametrize("from_float", apytypes_from_float_generators)
+def test_array_broadcast_arithmetic(bin_func, from_float, a_shape, b_shape):
+    np = pytest.importorskip("numpy")
+    a = np.arange(1, np.prod(a_shape) + 1).reshape(a_shape)
+    b = np.ones(b_shape)
+    ref = bin_func(a, b)
+    acc = bin_func(from_float(a), from_float(b))
+    assert np.all(ref == acc)