const nb::typed<nb::iterable, nb::any>& bit_pattern_sequence,

std::optional<int> int_bits,

std::optional<int> frac_bits,

std::optional<int> bits

)

: APyFixedArray(

python_iterable_extract_shape(bit_pattern_sequence, "APyFixedArray.__init__"),

int_bits,

frac_bits,

bits

)

{

// Specialized initialization for NDArray

if (nb::isinstance<nb::ndarray<>>(bit_pattern_sequence)) {

auto ndarray = nb::cast<nb::ndarray<nb::c_contig>>(bit_pattern_sequence);

_set_bits_from_ndarray(ndarray);

return; // initialization completed

}

// 1D vector of Python int object (`nb::int_` objects)

auto python_ints = python_iterable_walk<nb::int_>(

bit_pattern_sequence, "APyFixedArray.__init__"

);

for (std::size_t i = 0; i < _data.size() / _itemsize; i++) {

nb::int_ python_int = nb::cast<nb::int_>(python_ints[i]);

auto limb_vec = python_long_to_limb_vec(python_int, _itemsize);

_overflow_twos_complement(

std::begin(limb_vec), std::end(limb_vec), _bits, _int_bits

);

std::copy_n(std::begin(limb_vec), _itemsize, std::begin(_data) + i * _itemsize);

}

const std::vector<std::size_t>& shape, int bits, int int_bits

)

: APyArray(shape, bits_to_limbs(bits))

, _bits { bits }

, _int_bits { int_bits }

{

const std::vector<std::size_t>& shape, int bits, int int_bits, vector_type&& v

)

: APyArray(shape, bits_to_limbs(bits), std::move(v))

, _bits { bits }

, _int_bits { int_bits }

{

const std::vector<std::size_t>& shape,

std::optional<int> int_bits,

std::optional<int> frac_bits,

std::optional<int> bits

)

: APyArray(shape, bits_to_limbs(bits_from_optional(bits, int_bits, frac_bits)))

, _bits { bits.has_value() ? *bits : *int_bits + *frac_bits }

, _int_bits { int_bits.has_value() ? *int_bits : *bits - *frac_bits }

{

{

// Increase word length of result by one

const int res_int_bits = std::max(rhs.int_bits(), int_bits()) + 1;

const int res_frac_bits = std::max(rhs.frac_bits(), frac_bits());

const int res_bits = res_int_bits + res_frac_bits;

// Resulting vector

APyFixedArray result(_shape, res_bits, res_int_bits);

// Special case #1: Operands and result fit in single limb

if (unsigned(res_bits) <= APY_LIMB_SIZE_BITS) {

if (frac_bits() == rhs.frac_bits()) {

// Operands have equally many fractional bits.

simd_op {}(

_data.begin(),

rhs._data.begin(),

result._data.begin(),

result._data.size()

);

} else {

auto rhs_shift_amount = unsigned(res_frac_bits - rhs.frac_bits());

auto lhs_shift_amount = unsigned(res_frac_bits - frac_bits());

simd_shift_op {}(

_data.begin(),

rhs._data.begin(),

result._data.begin(),

lhs_shift_amount,

rhs_shift_amount,

result._data.size()

);

}

return result; // early exit

}

// Special case #2: Operands and result have equally many limbs

if (result._itemsize == _itemsize && result._itemsize == rhs._itemsize) {

const apy_limb_t* src1_ptr;

const apy_limb_t* src2_ptr;

if (frac_bits() == rhs.frac_bits()) {

// Right-hand side and left-hand side have equally many fractional bits

src1_ptr = &_data[0];

src2_ptr = &rhs._data[0];

} else if (frac_bits() <= rhs.frac_bits()) {

// Right-hand side has more fractional bits. Upsize `*this`

_cast_no_quantize_no_overflow(

std::begin(_data), // src

std::begin(result._data), // dst

_itemsize, // src_limbs

result._itemsize, // dst_limbs

_nitems, // n_items

result.frac_bits() - frac_bits() // left_shift_amount

);

src1_ptr = &result._data[0];

src2_ptr = &rhs._data[0];

} else {

// Left-hand side has more fractional bits. Upsize `rhs`

_cast_no_quantize_no_overflow(

std::begin(rhs._data), // src

std::begin(result._data), // dst

rhs._itemsize, // src_limbs

result._itemsize, // dst_limbs

rhs._nitems, // n_items

result.frac_bits() - rhs.frac_bits() // left_shift_amount

);

src1_ptr = &_data[0];

src2_ptr = &result._data[0];

}

for (std::size_t i = 0; i < result._data.size(); i += result._itemsize) {

ripple_carry_op {}(

&result._data[i], // dst

&src1_ptr[i], // src1

&src2_ptr[i], // src2

result._itemsize // limb vector length

);

}

return result; // early exit

}

// Most general case: Works in any situation, but is slowest

APyFixedArray imm(_shape, res_bits, res_int_bits);

_cast_no_quantize_no_overflow(

std::begin(_data), // src

std::begin(result._data), // dst

_itemsize, // src_limbs

result._itemsize, // dst_limbs

_nitems, // n_items

result.frac_bits() - frac_bits() // left_shift_amount

);

_cast_no_quantize_no_overflow(

std::begin(rhs._data), // src

std::begin(imm._data), // dst

rhs._itemsize, // src_limbs

imm._itemsize, // dst_limbs

rhs._nitems, // n_items

imm.frac_bits() - rhs.frac_bits() // left_shift_amount

);

// Perform ripple-carry operation for each element

for (std::size_t i = 0; i < result._data.size(); i += result._itemsize) {

ripple_carry_op {}(

&result._data[i], // dst

&result._data[i], // src1

&imm._data[i], // src2

result._itemsize // limb vector length

);

}

return result;

{

// Increase word length of result by one

const int res_int_bits = std::max(rhs.int_bits(), int_bits()) + 1;

const int res_frac_bits = std::max(rhs.frac_bits(), frac_bits());

const int res_bits = res_int_bits + res_frac_bits;

// Resulting vector

APyFixedArray result(_shape, res_bits, res_int_bits);

// Special case #1: Operands and result fit in single limb

if (unsigned(res_bits) <= APY_LIMB_SIZE_BITS) {

if (frac_bits() == rhs.frac_bits()) {

// Operands have equally many fractional bits.

simd_op_const {}(

_data.begin(), rhs._data[0], result._data.begin(), result._data.size()

);

} else {

auto rhs_shift_amount = unsigned(res_frac_bits - rhs.frac_bits());

auto lhs_shift_amount = unsigned(res_frac_bits - frac_bits());

simd_shift_op_const {}(

std::begin(_data), // src1

rhs._data[0] << rhs_shift_amount, // constant

std::begin(result._data), // dst

lhs_shift_amount, // src1 shift amount

result._data.size()

);

}

return result; // early exit

}

// Most general case: Works in any situation, but is slowest

APyFixed imm(res_bits, res_int_bits);

auto rhs_shift_amount = unsigned(res_frac_bits - rhs.frac_bits());

auto lhs_shift_amount = unsigned(res_frac_bits - frac_bits());

_cast_no_quantize_no_overflow(

std::begin(_data), // src

std::begin(result._data), // dst

_itemsize, // src_limbs

result._itemsize, // dst_limbs

_nitems, // n_items

lhs_shift_amount // left_shift_amount

);

_cast_no_quantize_no_overflow(

std::begin(rhs._data),

std::end(rhs._data),

std::begin(imm._data),

std::end(imm._data),

rhs_shift_amount

);

for (std::size_t i = 0; i < result._data.size(); i += result._itemsize) {

// Perform ripple-carry operation

ripple_carry_op {}(

&result._data[i], // dst

&result._data[i], // src1

&imm._data[0], // src2

result._itemsize // limb vector length

);

}

return result;

{

// Increase word length of result by one

const int res_int_bits = std::max(lhs.int_bits(), int_bits()) + 1;

const int res_frac_bits = std::max(lhs.frac_bits(), frac_bits());

const int res_bits = res_int_bits + res_frac_bits;

// Adjust binary point

APyFixedArray result(_shape, res_bits, res_int_bits);

_cast_no_quantize_no_overflow(

std::begin(_data), // src

std::begin(result._data), // dst

_itemsize, // src_limbs

result._itemsize, // dst_limbs

_nitems, // n_items

result.frac_bits() - frac_bits() // left_shift_amount

);

auto lhs_shift_amount = unsigned(res_frac_bits - lhs.frac_bits());

if (unsigned(res_bits) <= APY_LIMB_SIZE_BITS) {

apy_limb_t operand = lhs._data[0] << lhs_shift_amount;

simd::vector_rsub_const(

result._data.begin(), operand, result._data.begin(), result._data.size()

);

} else {

APyFixed imm(res_bits, res_int_bits);

_cast_no_quantize_no_overflow(

std::begin(lhs._data), // src_begin

std::end(lhs._data), // src_end

std::begin(imm._data), // dst_begin

std::end(imm._data), // dst_end

lhs_shift_amount

);

// Perform subtraction

for (std::size_t i = 0; i < result._data.size(); i += result._itemsize) {

apy_inplace_reversed_subtraction_same_length(

&result._data[i], // dst/src2

&imm._data[0], // src1

result._itemsize // limb vector length

);

}

}

// Return result

return result;

{

const int res_int_bits = lhs.int_bits() + frac_bits() + 1;

const int res_frac_bits = lhs.frac_bits() + int_bits();

const int res_bits = res_int_bits + res_frac_bits;

APyFixedArray result(_shape, res_bits, res_int_bits);

// Special case #1: The resulting number of bits fit in a single limb

if (unsigned(res_bits) <= APY_LIMB_SIZE_BITS) {

simd::vector_rdiv_const_signed(

std::begin(_data), // src2 (denominator)

lhs._data[0] << bits(), // src1 (constant numerator)

std::begin(result._data), // dst

_data.size() // vector elements

);

return result;

}

// General case: This always works but is slower than the special cases.

// Absolute value denominator

ScratchVector<apy_limb_t> abs_den(_itemsize);

// Absolute value left-shifted numerator. As `apy_unsigned_division` alters the

// numerator on call, we allocate twice its size. The first half [ `0`,

// `result._itemsize` ) is passed to `apy_unsigned_division` and [

// `result._item_size`, `2*result._itemsize` ) is used to cache the left-shifted

// absolute numerator.

ScratchVector<apy_limb_t> abs_num(2 * result._itemsize);

bool num_sign = limb_vector_abs(

std::begin(lhs._data),

std::end(lhs._data),

std::begin(abs_num) + result._itemsize

);

limb_vector_lsl(std::begin(abs_num) + result._itemsize, std::end(abs_num), bits());

for (std::size_t i = 0; i < _nitems; i++) {

bool den_sign = limb_vector_abs(

std::begin(_data) + (i + 0) * _itemsize,

std::begin(_data) + (i + 1) * _itemsize,

std::begin(abs_den)

);

std::copy_n(

std::begin(abs_num) + result._itemsize,

result._itemsize,

std::begin(abs_num)

);

// `apy_unsigned_division` requires the number of *significant* limbs in

// denominator

std::size_t den_significant_limbs

= significant_limbs(std::begin(abs_den), std::end(abs_den));

apy_unsigned_division(

&result._data[i * result._itemsize], // Quotient

&abs_num[0], // Numerator

result._itemsize, // Numerator limbs

&abs_den[0], // Denominator

den_significant_limbs // Denominator significant limbs

);

// Negate result if negative

if (num_sign ^ den_sign) {

limb_vector_negate_inplace(

std::begin(result._data) + (i + 0) * result._itemsize,

std::begin(result._data) + (i + 1) * result._itemsize

);

}

}

return result;

{

if (ndim() == 1 && rhs.ndim() == 1) {

if (_shape[0] == rhs._shape[0]) {

// Dimensionality for a standard scalar inner product checks out. Perform

// the checked inner product.

return _checked_inner_product(rhs, get_accumulator_mode_fixed());

}

}

if (ndim() == 2 && (rhs.ndim() == 2 || rhs.ndim() == 1)) {

if (_shape[1] == rhs._shape[0]) {

// Dimensionality for a standard 2D matrix multiplication checks out.

// Perform the checked 2D matrix

return _checked_2d_matmul(rhs, get_accumulator_mode_fixed());

}

}

// Unsupported `__matmul__` dimensionality, raise exception

throw std::length_error(

fmt::format(

"APyFixedArray.__matmul__: input shape mismatch, lhs: {}, rhs: {}",

tuple_string_from_vec(_shape),

tuple_string_from_vec(rhs._shape)

)

);

{

// Extract axes to sum over

std::vector<std::size_t> axes = cpp_axes_from_python(py_axis, _ndim);

// Retrieve how many elements will be summed together

std::size_t n_elems = array_fold_get_elements(axes);

// Compute the result word length

int pad_bits = n_elems ? bit_width(n_elems - 1) : 0;

int bits = _bits + pad_bits;

int int_bits = _int_bits + pad_bits;

std::size_t res_limbs = bits_to_limbs(bits);

// Accumulation function

auto fold = fold_accumulate<vector_type>(_itemsize, res_limbs);

return array_fold(axes, fold, std::nullopt /* no init */, bits, int_bits);

{

std::optional<std::size_t> axis = py_axis.has_value()

? std::optional<std::size_t>(nb::cast<std::size_t>(*py_axis))

: std::nullopt;

if (axis.has_value() && axis >= _ndim) {

std::string msg = fmt::format(

"APyFixedArray.cumsum: axis {} out of range (ndim = {})", *axis, _ndim

);

throw nb::index_error(msg.c_str());

}

// Retrieve how many elements will be summed together

std::size_t n_elems = axis.has_value() ? _shape[*axis] : _nitems;

// Compute the result word length

int pad_bits = n_elems ? bit_width(n_elems - 1) : 0;

int bits = _bits + pad_bits;

int int_bits = _int_bits + pad_bits;

std::size_t res_limbs = bits_to_limbs(bits);

// Accumulation function

auto fold = fold_accumulate<vector_type>(_itemsize, res_limbs);

auto post_proc = [](auto, auto) { /* no post processing */ };

return array_fold_cumulative(axis, fold, post_proc, std::nullopt, bits, int_bits);

{

// Extract axes to sum over

std::vector<std::size_t> axes = cpp_axes_from_python(py_axis, _ndim);

// Retrieve how many elements will be summed together

std::size_t n_elems = array_fold_get_elements(axes);

// Compute the result word length

if (n_elems == 0) {

// Empty array, return scalar one (NumPy semantics)

return APyFixed::one(_bits, _int_bits);

} else {

// Non-empty array

int int_bits = n_elems * _int_bits;

int frac_bits = n_elems * (_bits - _int_bits);

int bits = int_bits + frac_bits;

std::size_t res_limbs = bits_to_limbs(bits);

// Multiplicative fold function function

ScratchVector<apy_limb_t, 32> scratch(2 * res_limbs + 2 * _itemsize);

auto fold_func = fold_multiply<vector_type>(_itemsize, res_limbs, scratch);

APyFixed init_one(_bits, _int_bits, { apy_limb_t(1) });

return array_fold(axes, fold_func, init_one, bits, int_bits);

}

{

std::optional<std::size_t> axis = py_axis.has_value()

? std::optional<std::size_t>(nb::cast<std::size_t>(*py_axis))

: std::nullopt;

if (axis.has_value() && axis >= _ndim) {

std::string msg = fmt::format(

"APyFixedArray.cumprod: axis {} out of range (ndim = {})", *axis, _ndim

);

throw nb::index_error(msg.c_str());

}

// Retrieve how many elements will be folded together

std::size_t n_elems = axis.has_value() ? _shape[*axis] : _nitems;

// Compute the result word length

int int_bits = std::max(int(_int_bits * n_elems), _int_bits);

int frac_bits = std::max(int((_bits - _int_bits) * n_elems), _bits - _int_bits);

int bits = int_bits + frac_bits;

std::size_t res_limbs = bits_to_limbs(bits);

// Multiplicative fold function

ScratchVector<apy_limb_t, 32> scratch(2 * res_limbs + 2 * _itemsize);

auto fold_func = fold_multiply<vector_type>(_itemsize, res_limbs, scratch);

// Post processing: adjust the binary point of each partial product

auto post_proc = [&](auto dst_it, std::size_t n) {

int frac_bits = _bits - _int_bits;

int shift_amnt = frac_bits < 0 ? -frac_bits * n : frac_bits * (n_elems - 1 - n);

limb_vector_lsl(dst_it, dst_it + res_limbs, shift_amnt);

};

APyFixed init_one(_bits, _int_bits, { apy_limb_t(1) });

return array_fold_cumulative(axis, fold_func, post_proc, init_one, bits, int_bits);

{

if (ndim() != 1 || other.ndim() != 1) {

auto msg = fmt::format(

"can only convolve 1D arrays (lhs.ndim = {}, rhs.ndim = {})",

ndim(),

other.ndim()