This is a continuation from a previous question; I updated the code with the suggestions and added an additional iteration.
- Templates are not an option because of the usage: The objects are stored as a class members where the sizes are not known at compile time
- Iterators are not used, not even for
scan_kernel
because I could not find a way to introduce them; I'm still open to suggestions - Because of the new constructor I've had some trouble with the rule of 5, and I'm not sure about the duplication either, but there is a usecase where a new object needs to be constructed from another with a different padding so this is the best I could support it..
- The biggest problem with this is I do not know how to test
scan_kernel
; Otherwise how could this be improved?
Thank you very much for the great suggestions so far!
header:
#include <cstdint>
#include <vector>
#include <optional>
#include <functional>
class NDArrayIndex{
public:
NDArrayIndex(
const std::vector<std::uint32_t>& dimensions, const std::vector<std::int32_t>& padding = {},
const std::vector<std::uint32_t>& position = {}
);
NDArrayIndex(const NDArrayIndex& other, const std::vector<std::int32_t>& padding = {});
/** @brief Sets the position to all zeroes
*
* @return Reference to the object
*/
NDArrayIndex& reset(){
return set(std::vector<std::uint32_t>(size(), 0));
}
/** @brief updates the position of the object based on the given argument
*
* @param[in] position The position to move the object to; Must be within bounds!
*
* @return Reference to the object
*/
NDArrayIndex& set(const std::vector<std::uint32_t>& position);
/** @brief updates the position of the object based on the given arguments
*
* @param[in] dimension The dimension to set the position of; Must be within bounds!
* @param[in] position The position to move the object to; Must be within bounds!
*
* @return Reference to the object
*/
NDArrayIndex& set(std::uint32_t dimension, std::uint32_t position);
/** @brief updates the position of the object to go to the next position in the buffer range
*
* @return The index of the highest dimension the step modified
*/
std::uint32_t step();
/** @brief updates the position of the object to go to direction given in the arguments if the new position is inside bounds,
* throws an exception of the new position is not inside bounds
*
* @param[in] dimension The dimension to move the object on
* @param[in] delta The number of steps to move the objects position in the given dimension
*
* @return Reference to the object
*/
NDArrayIndex& step(std::uint32_t dimension, std::int32_t delta = 1);
const std::vector<std::uint32_t>& position() const{
return m_position;
}
/** @brief updates the position of the object based on the given argument
*
* @param[in] position The position to base the calculations on
*
* @return The value of the mapped index, if there is any
*/
std::optional<std::uint32_t> calculate_mapped_position(const std::vector<std::uint32_t>& position) const;
/** @brief Provides the index of the current position in the underlying buffer, if there is any
*
* @return the index pointing to the actual position, if it is mappable to the internal buffer
*/
std::optional<std::uint32_t> mapped_position() const{
return m_mappedIndex;
}
/** @brief Tells if the given position is inside the indexing interval provided by the dimensions and padding of the object
*
* @param[in] position The position to base the calculations on
* @param[in] dimension An optional dimension parameter to help shift the parameter
* @param[in] delta The number of steps to move the objects position temporarily in the given dimension for the relevant check
*
* @return true if the given position is inside bounds
*/
bool inside_bounds(const std::vector<std::uint32_t>& position, std::uint32_t dimension = 0u, std::int32_t delta = 0) const;
/** @brief Tells if the stored position is inside the indexing interval provided by the dimensions and padding of the object
*
* @param[in] dimension An optional dimension parameter to help shift the parameter
* @param[in] delta The number of steps to move the objects position temporarily in the given dimension for the relevant check
*
* @return true if the given position is inside bounds
*/
bool inside_bounds(std::uint32_t dimension = 0u, std::int32_t delta = 0) const{
return inside_bounds(m_position, dimension, delta);
}
/** @brief Tells if the given position is inside the indexing interval provided by the dimensions and padding of the object
*
* @param[in] index The index object supporting the position to base the check upon
* @param[in] dimension An optional dimension parameter to help shift the parameter
* @param[in] delta The number of steps to move the objects position temporarily in the given dimension for the relevant check
*
* @return true if the given position is inside bounds
*/
bool inside_bounds(const NDArrayIndex& index, std::uint32_t dimension = 0u, std::int32_t delta = 0) const{
return inside_bounds(index.position(), dimension, delta);
}
/** @brief Tells if the given position is mappable to the buffer range determined by the dimensions and padding
*
* @param[in] position The position to base the calculations on
* @param[in] dimension An optional dimension parameter to help shift the parameter
* @param[in] delta The number of steps to move the objects position temporarily in the given dimension for the relevant check
*
* @return true if the given position is inside bounds
*/
bool inside_content(const std::vector<std::uint32_t>& position, std::uint32_t dimension = 0u, std::int32_t delta = 0) const;
/** @brief Tells if the stored position is mappable to the buffer range determined by the dimensions and padding
*
* @param[in] dimension An optional dimension parameter to help shift the parameter
* @param[in] delta The number of steps to move the objects position temporarily in the given dimension for the relevant check
*
* @return true if the given position is inside bounds
*/
bool inside_content(std::uint32_t dimension = 0u, std::int32_t delta = 0) const{
return inside_content(m_position, dimension, delta);
}
/** @brief Tells if the given position is mappable to the buffer range determined by the dimensions and padding
*
* @param[in] index The index object supporting the position to base the check upon
* @param[in] dimension An optional dimension parameter to help shift the parameter
* @param[in] delta The number of steps to move the objects position temporarily in the given dimension for the relevant check
*
* @return true if the given position is inside bounds
*/
bool inside_content(const NDArrayIndex& index, std::uint32_t dimension = 0u, std::int32_t delta = 0) const{
return inside_content(index.position(), dimension, delta);
}
/** @struct IntervalPart
* @brief Describes part of an interval excluding the direction it lies on
* @var position_start the absolute starting position of the interval relevan part
* @var steps_inside_target the size of the interval's relevant part
*/
struct IntervalPart{
std::uint32_t position_start;
std::uint32_t steps_inside_target;
};
/** @brief Tells which parts of the provided range relative to the stored direction are mappable to
* the buffer range determined by the dimensions and padding
*
* @param[in] dimension The direction of the range relative to the currently stored position
* @param[in] delta The size of the range starting from the stored position
*
* @return A vector of the parts of the interval inside the bounds of the objects buffer range:
* {position, size}:
* |->absolute position inside the given dimension,
* |--------->number of steps still inside the defined ranges in the direction of the given dimension
*/
std::vector<IntervalPart> mappable_parts_of(
const std::vector<std::uint32_t>& position, std::uint32_t dimension, std::int32_t delta
) const;
/** @brief Tells which parts of the stored range relative to the stored direction are mappable to
* the buffer range determined by the dimensions and padding
*
* @param[in] dimension The direction of the range relative to the currently stored position
* @param[in] delta The size of the range starting from the stored position
*
* @return A vector of the parts of the interval inside the bounds of the objects buffer range:
* {position, size}:
* |->absolute position inside the given dimension,
* |--------->number of steps still inside the defined ranges in the direction of the given dimension
*/
std::vector<IntervalPart> mappable_parts_of(std::uint32_t dimension, std::int32_t delta) const{
return mappable_parts_of(m_position, dimension, delta);
}
/** @brief Tells the size of the internal buffer range, which maps every item in the NDArray to a one dimensional array
*
* @return Reference to the object
*/
std::uint32_t buffer_size() const{
return m_bufferSize;
}
/** @brief Returns the number of dimensions
*
* @return The number of dimensions
*/
std::uint32_t size() const{
return m_dimensions.size();
}
/** @brief Returns the number of elements inside bounds under the given dimension
*
* @param[in] dimension the dimension to query the size for
*
* @return the number of elements ( including padding ) the given dimension contains
*/
std::uint32_t operator[](std::int32_t dimension) const{
return m_padding[dimension] + m_dimensions[dimension] + m_padding[dimension];
}
/** @brief Tells if the Object contains any padding at all
*
* @return true, if any padding dimension is non-zero
*/
bool has_padding(){
return (static_cast<std::uint32_t>(std::count(m_padding.begin(), m_padding.end(), 0)) < m_padding.size());
}
/** @brief Runs a given function through a kernel, starting from the stored position.
* The provided function is being called once every time kernel iteration hits the
* beginning of dimension[0] in the provided kernel. The arguments are called with
* the mapped index values inside this object, along with the count of elements available
* from the start until the end of dimension[0]. The position of the object is updated along
* with the position of the kernel during iteration, and is restored after iteration is finished.
*
* @param kernel the kernel dimensions to use for iteration
* @param[in] fun the function to call for each kernel iteration.
* Arguments: void(mapped_position, interval size)
*/
void scan_kernel(NDArrayIndex& kernel, std::function<void(std::uint32_t, std::uint32_t)> fun);
private:
const std::vector<std::uint32_t> m_dimensions;
const std::vector<std::int32_t> m_padding;
const std::vector<std::uint32_t> m_strides;
const std::uint32_t m_bufferSize;
std::vector<std::uint32_t> m_position;
std::optional<std::uint32_t> m_mappedIndex;
};
source:
#include <cassert>
#include <numeric>
#include <algorithm>
#include <cmath>
std::vector<std::int32_t> init_padding(
const std::vector<std::uint32_t>& dimensions, const std::vector<std::int32_t>& padding
){
if(1 == padding.size())
return std::vector<std::int32_t>(dimensions.size(), padding[0]);
if(1 < padding.size()){
assert(dimensions.size() == padding.size());
return padding;
}
return std::vector<std::int32_t>(dimensions.size(), 0);
}
std::vector<std::uint32_t> init_strides(
const std::vector<std::uint32_t>& dimensions, const std::vector<std::int32_t>& padding
){
assert(dimensions.size() == padding.size());
std::vector<std::uint32_t> strides;
std::uint32_t prev_stride = 1u;
std::int32_t prev_padding = padding[0];
std::uint32_t dim = 0;
for(const std::uint32_t& dimension : dimensions){
strides.push_back(prev_stride);
prev_stride *= dimension + 2 * std::min(0, prev_padding);
prev_padding = padding[dim++];
}
return strides;
}
std::vector<std::uint32_t> init_position(
const std::vector<std::uint32_t>& dimensions, const std::vector<std::uint32_t>& position
){
if(0 < position.size()){
assert(dimensions.size() == position.size());
return {position};
}
return std::vector<std::uint32_t>(dimensions.size(), 0);
}
} /* namespace */
namespace rafko_utilities {
NDArrayIndex::NDArrayIndex(
const std::vector<std::uint32_t>& dimensions, const std::vector<std::int32_t>& padding,
const std::vector<std::uint32_t>& position
)
: m_dimensions(dimensions)
, m_padding(init_padding(m_dimensions, padding))
, m_strides(init_strides(dimensions, m_padding))
, m_bufferSize(std::accumulate(m_dimensions.begin(), m_dimensions.end(), 1.0,
[](const std::uint32_t& partial, const std::uint32_t& element){ return partial * element; }
))
, m_position(init_position(m_dimensions, position))
, m_mappedIndex(calculate_mapped_position(m_position))
{
assert(0 == std::count(m_dimensions.begin(), m_dimensions.end(), 0));
assert(inside_bounds(m_position));
}
NDArrayIndex::NDArrayIndex(const NDArrayIndex& other, const std::vector<std::int32_t>& padding)
: m_dimensions(other.m_dimensions)
, m_padding(init_padding(m_dimensions, padding))
, m_strides(init_strides(m_dimensions, m_padding))
, m_bufferSize(other.m_bufferSize)
, m_position(other.m_position)
, m_mappedIndex(calculate_mapped_position(m_position))
{
assert(0 == std::count(m_dimensions.begin(), m_dimensions.end(), 0));
assert(inside_bounds(m_position));
}
NDArrayIndex& NDArrayIndex::set(const std::vector<std::uint32_t>& position){
assert(position.size() == m_position.size());
assert(inside_bounds(position));
m_position = position;
m_mappedIndex = calculate_mapped_position(m_position);
assert( (!m_mappedIndex.has_value())||(m_mappedIndex.value() < m_bufferSize) );
return *this;
}
NDArrayIndex& NDArrayIndex::set(std::uint32_t dimension, std::uint32_t position){
assert(dimension < size());
m_position[dimension] = position;
m_mappedIndex = calculate_mapped_position(m_position);
assert( (!m_mappedIndex.has_value())||(m_mappedIndex.value() < m_bufferSize) );
assert(inside_bounds(m_position));
return *this;
}
std::uint32_t NDArrayIndex::step(){
std::uint32_t dim = 0;
bool changed = false;
while(dim < m_dimensions.size()){
if(inside_bounds(dim, 1)){
step(dim, 1);
break;
}else{
changed = true;
m_position[dim] = 0u;
}
++dim;
}
if(dim >= m_dimensions.size()){
m_mappedIndex = 0u; /* Overflow happened, start from the beginning */
return m_dimensions.size() - 1;
}else{
if(changed)m_mappedIndex = calculate_mapped_position(m_position);
assert(m_mappedIndex < m_bufferSize);
return dim;
}
}
NDArrayIndex& NDArrayIndex::step(std::uint32_t dimension, std::int32_t delta){
const std::int32_t new_position = static_cast<std::int32_t>(m_position[dimension]) + delta;
if(
(new_position < 0)
||(
new_position >= static_cast<std::int32_t>(m_dimensions[dimension] + (2 * std::max(0, m_padding[dimension])))
)
)throw std::runtime_error(
"Current position d[" + std::to_string(dimension) + "] + " + std::to_string(delta) + " out of bounds!"
);
m_position[dimension] = new_position;
bool new_position_is_inside_content = inside_content(m_position);
if(m_mappedIndex.has_value() && new_position_is_inside_content){ /* m_mappedIndex has a value if the previous position was valid */
m_mappedIndex.value() += m_strides[dimension] * delta;
assert(m_mappedIndex < m_bufferSize);
}else if(new_position_is_inside_content){ /* if the new position is inside bounds, then the mapped index can be caluclated */
m_mappedIndex = calculate_mapped_position(m_position);
}else m_mappedIndex = {}; /* No mapped index for positions inside the padding */
return *this;
}
std::optional<std::uint32_t> NDArrayIndex::calculate_mapped_position(const std::vector<std::uint32_t>& position) const{
assert(position.size() == m_strides.size());
if(!inside_content(position))
return {};
std::uint32_t result_index = 0u;
for(std::uint32_t dim = 0; dim < position.size(); ++dim){
result_index += (position[dim] - std::max(m_padding[dim], -m_padding[dim])) * m_strides[dim];
}
return result_index;
}
bool NDArrayIndex::inside_bounds(const std::vector<std::uint32_t>& position, std::uint32_t dimension, std::int32_t delta) const{
std::uint32_t dim = 0;
return std::all_of(position.begin(), position.end(),
[this, &dim, dimension, delta](const std::uint32_t& pos){
std::int32_t position = static_cast<std::int32_t>(pos);
if(dim == dimension) position += delta;
++dim;
return(
(0 <= position)
&&( position < (2 * std::max(int64_t{0}, static_cast<int64_t>(m_padding[dim - 1])) + static_cast<int64_t>(m_dimensions[dim - 1])) )
);
}
);
}
bool NDArrayIndex::inside_content(const std::vector<std::uint32_t>& position, std::uint32_t dimension, std::int32_t delta) const{
std::uint32_t dim = 0;
return std::all_of(position.begin(), position.end(),
[this, &dim, dimension, delta](const std::uint32_t& pos){
std::int32_t actual_position = static_cast<std::int32_t>(pos);
if(dim == dimension) actual_position += delta;
++dim;
return(
(std::max(m_padding[dim - 1], -m_padding[dim - 1]) <= actual_position)
&&(actual_position < static_cast<std::int32_t>(m_dimensions[dim - 1] + m_padding[dim - 1]))
);
}
);
}
std::vector<NDArrayIndex::IntervalPart> NDArrayIndex::mappable_parts_of(
const std::vector<std::uint32_t>& position, std::uint32_t dimension, std::int32_t delta
) const{
std::vector<NDArrayIndex::IntervalPart> result;
bool part_in_progress = false;
for(std::int32_t delta_index = 0; delta_index < delta; delta_index += std::copysign(1, delta)){
const bool current_position_in_inside_content = inside_content(position, dimension, delta_index);
if(current_position_in_inside_content && part_in_progress){
assert(0 < result.size());
++result.back().steps_inside_target; /* Increase the size of the current part of the interval */
}else if(current_position_in_inside_content){ /* If the interval iteration became inside bounds */
result.push_back({(position[dimension] + delta_index), 1}); /* Add the new part as a result */
part_in_progress = true;
}else part_in_progress = false;
}
return result;
}
void NDArrayIndex::scan_kernel(NDArrayIndex& kernel, std::function<void(std::uint32_t, std::uint32_t)> fun){
assert(!kernel.has_padding());
assert(size() == kernel.size());
std::vector<std::uint32_t> original_position = m_position;
kernel.reset();
do{ /* Acquire interval inside bounds of the mappable buffer and call the function with it */
std::vector<NDArrayIndex::IntervalPart> parts_inside_content = mappable_parts_of(m_position, 0, kernel[0]);
if(mapped_position().has_value() && (0 < parts_inside_content.size())){
auto& [start_pos, interval_size] = parts_inside_content[0];
fun(mapped_position().value() - m_position[0] + start_pos, interval_size);
}
try{ /* Step in dimension[1], because dimension[0] should stay at position 0 */
kernel.step(1,1); /* throws exception if the position goes out of bounds */
step(1,1);
std::cout << "^";
}catch(...){ /* In case dimension[1] is out of bounds, step to the next position by the interface with wrap around */
std::cout << ">";
kernel.set(0, kernel[0] - 1); /* set the first dimensions position inside the kernel to the last */
std::uint32_t modified_dimension = kernel.step();
/*!Note: Because of overflow, the position in dimension[0] will reset to zero */
if((kernel.mapped_position().has_value())&&(kernel.mapped_position().value() != 0)){
for(std::uint32_t dim = 0; dim < modified_dimension; ++dim){
m_position[dim] = original_position[dim];
}
++m_position[modified_dimension];
m_mappedIndex = calculate_mapped_position(m_position);
}
}
}while( /* when the mapped position for the kernel points to the start of the first dimension, and the end of the others */
(kernel.mapped_position().has_value()) /* the kernel is iterated through */
&&(0 != kernel.mapped_position().value())
);
set(original_position);
}
test:
TEST_CASE("Testing NDArray Indexing with a 2D array without padding", "[NDArray]"){
std::uint32_t width = rand()%100;
std::uint32_t height = rand()%100;
rafko_utilities::NDArrayIndex idx({width, height});
REQUIRE(!idx.has_padding());
for(std::uint32_t variant = 0; variant < 5; ++variant){
std::uint32_t x = rand()%width;
std::uint32_t y = rand()%height;
idx.set({x,y});
REQUIRE(idx.inside_bounds());
REQUIRE(idx.mapped_position().has_value());
REQUIRE(idx.mapped_position().value() == (x + (y * width)));
std::uint32_t elements_after_x_row = width - x;
REQUIRE(1 == idx.mappable_parts_of(0,width).size());
REQUIRE(x == idx.mappable_parts_of(0,width)[0].position_start);
REQUIRE(elements_after_x_row == idx.mappable_parts_of(0,width)[0].steps_inside_target);
/*!Note: using width in the above interfaces because it is guaranteed
* that an interval of that size spans over the relevant dimension
* */
if(y < (height - 1u))REQUIRE( idx.step(1,1).mapped_position() == (x + ((y + 1) * width)) );
else CHECK_THROWS(idx.step(1,1));
}
REQUIRE(idx.buffer_size() == (width * height));
idx.set({0,0});
for(std::uint32_t i = 0; i < idx.buffer_size(); ++i){
REQUIRE(idx.inside_bounds());
REQUIRE(idx.inside_content());
REQUIRE(idx.mapped_position().has_value() == true);
REQUIRE(idx.mapped_position().value() == i);
idx.step();
}
}
TEST_CASE("Testing NDArray Indexing with a 2D array with positive padding", "[NDArray][padding]"){
std::uint32_t width = 1 + rand()%20;
std::uint32_t height = 1 + rand()%20;
std::int32_t padding_x = rand()%5;
std::int32_t padding_y = rand()%5;
rafko_utilities::NDArrayIndex idx({width, height}, {padding_x, padding_y});
REQUIRE(idx.has_padding());
for(std::uint32_t variant = 0; variant < 5; ++variant){
std::uint32_t x = padding_x + rand()%(width);
std::uint32_t y = padding_y + rand()%(height);
idx.set({x,y});
REQUIRE(idx.inside_bounds());
REQUIRE(idx.mapped_position().has_value());
REQUIRE( idx.mapped_position().value() == (x - padding_x + ((y - padding_y) * width)) );
std::uint32_t elements_after_x_row = padding_x + width - x;
REQUIRE(1 == idx.mappable_parts_of(0,width).size());
REQUIRE(x == idx.mappable_parts_of(0,width)[0].position_start);
REQUIRE(elements_after_x_row == idx.mappable_parts_of(0,width)[0].steps_inside_target);
if((static_cast<std::int32_t>(y) >= padding_y) && (y < (height + padding_y - 1)))
REQUIRE( idx.step(1,1).mapped_position() == (x - padding_x + ((y - padding_y + 1) * width)) );
else CHECK_NOTHROW(idx.step(1,1));
}
REQUIRE(idx.buffer_size() == (width * height));
std::uint32_t x = 0u;
std::uint32_t y = 0u;
std::uint32_t reference_mapped_position = 0u;
idx.set({0,0});
for(std::uint32_t i = 0; i < idx.buffer_size(); ++i){
if(
(padding_x <= static_cast<std::int32_t>(x) && x < (padding_x + width))
&&(padding_y <= static_cast<std::int32_t>(y) && y < (padding_y + height))
){
REQUIRE(idx.inside_bounds());
REQUIRE(idx.inside_content());
REQUIRE(idx.mapped_position().has_value() == true);
REQUIRE(idx.mapped_position().value() == reference_mapped_position);
++reference_mapped_position;
}else{
REQUIRE(idx.inside_bounds());
REQUIRE(idx.mapped_position().has_value() == false);
}
if(x < padding_x + width + padding_x - 1){
REQUIRE(idx.step() == 0u);
++x;
}else{
REQUIRE(idx.step() == 1u);
x = 0;
++y;
}
}
}
TEST_CASE("Testing NDArray Indexing with a 2D array with negative padding", "[NDArray][padding]"){
std::uint32_t width = 11 + rand()%20;
std::uint32_t height = 11 + rand()%20;
std::int32_t padding_x = -rand()%5;
std::int32_t padding_y = -rand()%5;
rafko_utilities::NDArrayIndex idx({width, height}, {padding_x, padding_y});
REQUIRE(idx.has_padding());
for(std::uint32_t variant = 0; variant < 5; ++variant){
std::uint32_t x = -padding_x + rand()%(width + 2 * padding_x);
std::uint32_t y = -padding_y + rand()%(height + 2 * padding_y);
idx.set({x,y});
REQUIRE(idx.inside_bounds());
REQUIRE(idx.mapped_position().has_value());
REQUIRE( idx.mapped_position().value() == (x + padding_x + ((y + padding_y) * (width + 2 * padding_x))) );
std::uint32_t elements_after_x_row = padding_x + width - x;
REQUIRE(1 == idx.mappable_parts_of(0,width).size());
REQUIRE(x == idx.mappable_parts_of(0,width)[0].position_start);
REQUIRE(elements_after_x_row == idx.mappable_parts_of(0,width)[0].steps_inside_target);
if((static_cast<std::int32_t>(y) > -padding_y) && (y < (height + padding_y - 1)))
REQUIRE( idx.step(1,1).mapped_position() == (x + padding_x + ((y + padding_y + 1) * (width + 2 * padding_x))) );
else CHECK_NOTHROW(idx.step(1,1));
}
REQUIRE(idx.buffer_size() == (width * height));
std::uint32_t x = 0u;
std::uint32_t y = 0u;
std::uint32_t reference_mapped_position = 0u;
idx.set({0,0});
for(std::uint32_t i = 0; i < idx.buffer_size(); ++i){
if(
(-padding_x <= static_cast<std::int32_t>(x) && x < (padding_x + width))
&&(-padding_y <= static_cast<std::int32_t>(y) && y < (padding_y + height))
){
REQUIRE(idx.inside_bounds());
REQUIRE(idx.inside_content());
REQUIRE(idx.mapped_position().has_value() == true);
REQUIRE(idx.mapped_position().value() == reference_mapped_position);
++reference_mapped_position;
}else{
REQUIRE(idx.inside_bounds());
REQUIRE(idx.mapped_position().has_value() == false);
}
if(x < (width - 1)){
REQUIRE(idx.step() == 0u);
++x;
}else{
REQUIRE(idx.step() == 1u);
x = 0;
++y;
}
}
}
std::mdspan
. You might want to see if that is something that would be useful for you, and implement something with a similar API. Also have a look at its reference implementation. \$\endgroup\$