Arithmetic.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <stratax/core/dtypes/Concepts.hpp>
#include <stratax/exceptions/Exceptions.hpp>
#include <stratax/ops/Broadcasting.hpp>
#include <stratax/ops/detail/Compound.hpp>
#include <stratax/core/validation/NumericValidation.hpp>
 
#include <functional>
 
template<Array L, Array R, typename Op>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
auto binary_op(
    const L& lhs,
    const R& rhs,
    Op op,
    bool check_zero_divisor = false)
{
    auto checked_op = [&](const auto& left, const auto& right)
    {
        if (check_zero_divisor)
        {
            stratax::core::validation::require_valid_division(left, right);
        }
 
        return op(left, right);
    };
 
    using result_value_type =
        stratax::core::promote_t<
            typename L::value_type,
            typename R::value_type>;
 
    return stratax::core::broadcasted_op<result_value_type>(
        lhs,
        rhs,
        checked_op);
}
 
template<Array A, Numeric Scalar, typename Op>
auto binary_scalar_op(const A& lhs, const Scalar& rhs, Op op, bool check_zero_divisor = false)
{
    if (check_zero_divisor && rhs == Scalar{})
    {
        throw Exceptions::ZeroDivisionError("Division by zero.");
    }
 
    auto checked_op = [&](const auto& left, const auto& right)
    {
        if (check_zero_divisor)
        {
            stratax::core::validation::require_valid_division(left, right);
        }
        return op(left, right);
    };
    return broadcasted_op(lhs, rhs, checked_op);
}
 
template<Numeric Scalar, Array A, typename Op>
auto binary_scalar_op(const Scalar& lhs, const A& rhs, Op op, bool check_zero_divisor = false)
{
    auto checked_op = [&](const auto& left, const auto& right)
    {
        if (check_zero_divisor)
        {
            stratax::core::validation::require_valid_division(left, right);
        }
 
        return op(left, right);
    };
 
    return broadcasted_op(lhs, rhs, checked_op);
}
 
template<Array L, Array R, typename Op>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
L& compound_op(
    L& lhs,
    const R& rhs,
    Op op,
    bool check_zero_divisor = false)
{
    const auto result_shape =
        stratax::core::broadcast_detail::array_result_shape(lhs, rhs);
 
    // Compound assignment cannot change the lhs shape.
    if (result_shape != lhs.shape())
    {
        throw Exceptions::BroadcastError(
            "In-place broadcasting cannot change the left operand's shape.");
    }
 
    // Validate all used divisors before changing any elements.
    if (check_zero_divisor)
    {
        for (std::size_t i = 0; i < lhs.size(); ++i)
        {
            const auto rhs_index = stratax::core::broadcast_detail::flat_operand_index(
                i, lhs.shape(), rhs.shape());
            stratax::core::validation::require_valid_division(lhs[i], rhs[rhs_index]);
        }
    }
 
    constexpr bool stage = !stratax::core::compound_detail::owns_storage<L> ||
        !stratax::core::compound_detail::owns_storage<R>;
    return stratax::core::compound_detail::write_results<stage>(lhs, [&](std::size_t i)
    {
        const auto rhs_index = stratax::core::broadcast_detail::flat_operand_index(
            i, lhs.shape(), rhs.shape());
        return op(lhs[i], rhs[rhs_index]);
    });
}
 
template<Array A, Numeric S, typename Op>
A& compound_scalar_op(
    A& lhs,
    const S& rhs,
    Op op,
    bool check_zero_divisor = false)
{
    const S scalar = rhs;
    if (check_zero_divisor)
    {
        if (scalar == S{})
        {
            throw Exceptions::ZeroDivisionError("Division by zero.");
        }
        for (std::size_t i = 0; i < lhs.size(); ++i)
        {
            stratax::core::validation::require_valid_division(lhs[i], scalar);
        }
    }
 
    return stratax::core::compound_detail::write_results<
        !stratax::core::compound_detail::owns_storage<A>>(lhs, [&](std::size_t i)
    {
        return op(lhs[i], scalar);
    });
}
 
template<Array L, Array R>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
auto operator+(const L& lhs, const R& rhs)
{
    return binary_op(lhs, rhs, std::plus<>{});
}
 
template<Array L, Array R>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
auto operator-(const L& lhs, const R& rhs)
{
    return binary_op(lhs, rhs, std::minus<>{});
}
 
template<Array L, Array R>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
auto operator*(const L& lhs, const R& rhs)
{
    return binary_op(lhs, rhs, std::multiplies<>{});
}
 
template<Array L, Array R>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
auto operator/(const L& lhs, const R& rhs)
{
    return binary_op(lhs, rhs, std::divides<>{}, true);
}
 
template<Array A, Numeric Scalar>
requires Numeric<typename A::value_type>
auto operator+(const A& lhs, const Scalar& rhs)
{
    return binary_scalar_op(lhs, rhs, std::plus<>{});
}
 
template<Array A, Numeric Scalar>
requires Numeric<typename A::value_type>
auto operator-(const A& lhs, const Scalar& rhs)
{
    return binary_scalar_op(lhs, rhs, std::minus<>{});
}
 
template<Array A, Numeric Scalar>
requires Numeric<typename A::value_type>
auto operator*(const A& lhs, const Scalar& rhs)
{
    return binary_scalar_op(lhs, rhs, std::multiplies<>{});
}
 
template<Array A, Numeric Scalar>
requires Numeric<typename A::value_type>
auto operator/(const A& lhs, const Scalar& rhs)
{
    return binary_scalar_op(lhs, rhs, std::divides<>{}, true);
}
 
template<Numeric Scalar, Array A>
requires Numeric<typename A::value_type>
auto operator+(const Scalar& lhs, const A& rhs)
{
    return rhs + lhs;
}
 
template<Numeric Scalar, Array A>
auto operator-(const Scalar& lhs, const A& rhs)
{
    return binary_scalar_op(lhs, rhs, std::minus<>{});
}
 
template<Numeric Scalar, Array A>
requires Numeric<typename A::value_type>
auto operator*(const Scalar& lhs, const A& rhs)
{
    return rhs * lhs;
}
 
template<Numeric Scalar, Array A>
requires Numeric<typename A::value_type>
auto operator/(const Scalar& lhs, const A& rhs)
{
    return binary_scalar_op(lhs, rhs, std::divides<>{}, true);
}
 
template<Array L, Array R>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
L& operator+=(L& lhs, const R& rhs)
{
    return compound_op(lhs, rhs, std::plus<>{});
}
 
template<Array L, Array R>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
L& operator-=(L& lhs, const R& rhs)
{
    return compound_op(lhs, rhs, std::minus<>{});
}
 
template<Array L, Array R>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
L& operator*=(L& lhs, const R& rhs)
{
    return compound_op(lhs, rhs, std::multiplies<>{});
}
 
template<Array L, Array R>
requires (
    Numeric<typename L::value_type> &&
    Numeric<typename R::value_type>
)
L& operator/=(L& lhs, const R& rhs)
{
    return compound_op(lhs, rhs, std::divides<>{}, true);
}
 
template<Array A, Numeric S>
requires Numeric<typename A::value_type>
A& operator+=(A& lhs, const S& rhs)
{
    return compound_scalar_op(lhs, rhs, std::plus<>{});
}
 
template<Array A, Numeric S>
requires Numeric<typename A::value_type>
A& operator-=(A& lhs, const S& rhs)
{
    return compound_scalar_op(lhs, rhs, std::minus<>{});
}
 
template<Array A, Numeric S>
requires Numeric<typename A::value_type>
A& operator*=(A& lhs, const S& rhs)
{
    return compound_scalar_op(lhs, rhs, std::multiplies<>{});
}
 
template<Array A, Numeric S>
requires Numeric<typename A::value_type>
A& operator/=(A& lhs, const S& rhs)
{
    return compound_scalar_op(lhs, rhs, std::divides<>{}, true);
}
 
template<Array A>
requires Numeric<typename A::value_type>
A operator-(const A& arr)
{
    return arr * typename A::value_type{-1};
}
 
template<Array A>
requires Numeric<typename A::value_type>
A operator+(const A& arr)
{
    return arr;
}