de_infinite_integrator.h

Go to the documentation of this file.

/*
 * Copyright 2021 MusicScience37 (Kenta Kabashima)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
 
#include <cmath>
#include <type_traits>
 
#include "num_collect/base/concepts/invocable_as.h"
#include "num_collect/base/concepts/real_scalar.h"
#include "num_collect/base/exception.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/precondition.h"
#include "num_collect/constants/half.h"  // IWYU pragma: keep
#include "num_collect/constants/one.h"   // IWYU pragma: keep
#include "num_collect/constants/pi.h"    // IWYU pragma: keep
#include "num_collect/constants/zero.h"  // IWYU pragma: keep
#include "num_collect/logging/logging_macros.h"
#include "num_collect/logging/logging_mixin.h"
 
namespace num_collect::integration {

constexpr auto de_infinite_integrator_tag =
    logging::log_tag_view("num_collect::integration::de_infinite_integrator");

template <typename Signature>
class de_infinite_integrator;

template <typename Result, base::concepts::real_scalar Variable>
class de_infinite_integrator<Result(Variable)> : public logging::logging_mixin {
public:
    using variable_type = std::decay_t<Variable>;

    using result_type = std::decay_t<Result>;

    de_infinite_integrator()
        : logging::logging_mixin(de_infinite_integrator_tag) {}

    template <base::concepts::invocable_as<result_type(variable_type)> Function>
    [[nodiscard]] auto integrate(const Function& function) const
        -> result_type {
        const variable_type interval =
            max_point_ / static_cast<variable_type>(points_);
 
        constexpr variable_type diff_coeff_center = half_pi;
        result_type sum =
            function(constants::zero<variable_type>) * diff_coeff_center;
 
        for (index_type i = 1; i < points_; ++i) {
            const variable_type changed_var =
                interval * static_cast<variable_type>(i);
            const variable_type half_pi_sinh = half_pi * std::sinh(changed_var);
            const variable_type var = std::sinh(half_pi_sinh);
            const variable_type diff_coeff =
                half_pi * std::cosh(half_pi_sinh) * std::cosh(changed_var);
            sum += function(var) * diff_coeff;
            sum += function(-var) * diff_coeff;
        }
 
        return sum * interval;
    }

    template <base::concepts::invocable_as<result_type(variable_type)> Function>
    [[nodiscard]] auto operator()(const Function& function) const
        -> result_type {
        return integrate(function);
    }

    auto max_point(variable_type val) -> de_infinite_integrator& {
        NUM_COLLECT_PRECONDITION(val > static_cast<variable_type>(0),
            this->logger(), "Maximum point must be a positive value.");
        max_point_ = val;
        return *this;
    }

    auto points(index_type val) -> de_infinite_integrator& {
        NUM_COLLECT_PRECONDITION(val > 0, this->logger(),
            "Number of points must be a positive integer.");
        points_ = val;
        return *this;
    }
 
private:
    static constexpr variable_type half_pi =
        constants::half<variable_type> * constants::pi<variable_type>;

    static constexpr auto default_max_point = static_cast<variable_type>(4);

    variable_type max_point_{default_max_point};

    static constexpr index_type default_points = 20;

    index_type points_{default_points};
};
 
}  // namespace num_collect::integration