product_integral.hpp

Go to the documentation of this file.

#pragma once
 
#include <type_traits>
#include <cmath>
#include <utility>
#include "../algorithms/core/result.hpp"
#include "analysis.hpp"
 
namespace limes::expr {
 
// Forward declarations
template<typename E, std::size_t Dim, typename Lo, typename Hi>
struct Integral;
 
template<typename E, std::size_t Dims>
struct BoxIntegral;
 
template<typename I1, typename I2>
struct ProductIntegral;
 
// Helper to get variable set from any evaluable integral type
namespace detail {
 
template<typename T, typename = void>
struct integral_variable_set {
    static constexpr std::uint64_t value = 0;
};
 
// For types with integrand_type (Integral, BoxIntegral)
template<typename T>
struct integral_variable_set<T, std::void_t<typename T::integrand_type>> {
    static constexpr std::uint64_t value = variable_set_v<typename T::integrand_type>;
};
 
// For ProductIntegral: combine the variable sets of both children
template<typename I1, typename I2>
struct integral_variable_set<ProductIntegral<I1, I2>> {
    static constexpr std::uint64_t value =
        integral_variable_set<I1>::value | integral_variable_set<I2>::value;
};
 
template<typename T>
inline constexpr std::uint64_t integral_variable_set_v = integral_variable_set<T>::value;
 
} // namespace detail
 
// =============================================================================
// ProductIntegral: Composition of independent integrals
// =============================================================================
 
template<typename I1, typename I2>
struct ProductIntegral {
    using value_type = typename I1::value_type;
 
    I1 left;
    I2 right;
 
    // Verify independence at compile time
    // The integrals must depend on disjoint sets of variables
    static_assert(
        (detail::integral_variable_set_v<I1> & detail::integral_variable_set_v<I2>) == 0,
        "ProductIntegral requires integrals over independent variables. "
        "The integrand variable sets must be disjoint."
    );
 
    constexpr ProductIntegral(I1 i1, I2 i2) noexcept
        : left{i1}, right{i2} {}
 
    [[nodiscard]] algorithms::integration_result<value_type>
    eval() const {
        return combine_results(left.eval(), right.eval());
    }
 
    template<typename... Args>
    [[nodiscard]] algorithms::integration_result<value_type>
    eval(Args&&... args) const {
        return combine_results(
            left.eval(std::forward<Args>(args)...),
            right.eval(std::forward<Args>(args)...));
    }
 
private:
    static algorithms::integration_result<value_type>
    combine_results(algorithms::integration_result<value_type> const& r1,
                    algorithms::integration_result<value_type> const& r2) {
        value_type value = r1.value() * r2.value();
        value_type error = std::abs(r2.value()) * r1.error()
                         + std::abs(r1.value()) * r2.error();
        return {value, error,
                r1.iterations() + r2.iterations(),
                r1.evaluations() + r2.evaluations()};
    }
 
public:
 
    [[nodiscard]] std::string to_string() const {
        return "(product-integral " + left.to_string() + " " + right.to_string() + ")";
    }
 
    template<typename I3>
    [[nodiscard]] constexpr auto operator*(I3 const& other) const {
        return ProductIntegral<ProductIntegral, I3>{*this, other};
    }
};
 
// =============================================================================
// Type traits for ProductIntegral
// =============================================================================
 
template<typename T>
struct is_product_integral : std::false_type {};
 
template<typename I1, typename I2>
struct is_product_integral<ProductIntegral<I1, I2>> : std::true_type {};
 
template<typename T>
inline constexpr bool is_product_integral_v = is_product_integral<T>::value;
 
// =============================================================================
// Concepts for integral types
// =============================================================================
 
template<typename T>
concept EvaluableIntegral = requires(T t) {
    { t.eval() } -> std::same_as<algorithms::integration_result<typename T::value_type>>;
};
 
// =============================================================================
// operator* for integrals (creates ProductIntegral)
// =============================================================================
 
template<typename I1, typename I2>
    requires (is_integral_v<I1> && is_integral_v<I2>)
[[nodiscard]] constexpr auto operator*(I1 const& i1, I2 const& i2) {
    return ProductIntegral<I1, I2>{i1, i2};
}
 
template<typename I1, typename I2, typename I3>
    requires is_integral_v<I3>
[[nodiscard]] constexpr auto operator*(ProductIntegral<I1, I2> const& pi, I3 const& i3) {
    return ProductIntegral<ProductIntegral<I1, I2>, I3>{pi, i3};
}
 
template<typename I1, typename I2, typename I3>
    requires is_integral_v<I3>
[[nodiscard]] constexpr auto operator*(I3 const& i3, ProductIntegral<I1, I2> const& pi) {
    return ProductIntegral<I3, ProductIntegral<I1, I2>>{i3, pi};
}
 
// =============================================================================
// Separable integral detection and optimization
// =============================================================================
 
template<typename I1, typename I2>
inline constexpr bool are_independent_integrals_v =
    (detail::integral_variable_set_v<I1> & detail::integral_variable_set_v<I2>) == 0;
 
// =============================================================================
// Convenience: product() function
// =============================================================================
 
template<typename I1, typename I2>
    requires (EvaluableIntegral<I1> && EvaluableIntegral<I2>)
[[nodiscard]] constexpr auto product(I1 const& i1, I2 const& i2) {
    return ProductIntegral<I1, I2>{i1, i2};
}
 
template<typename I1, typename I2, typename... Is>
    requires (EvaluableIntegral<I1> && EvaluableIntegral<I2> && (EvaluableIntegral<Is> && ...))
[[nodiscard]] constexpr auto product(I1 const& i1, I2 const& i2, Is const&... is) {
    return product(ProductIntegral<I1, I2>{i1, i2}, is...);
}
 
 
} // namespace limes::expr