legendre_roots.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 <limits>
#include <tuple>
 
#include <Eigen/Core>
 
#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/pi.h"    // IWYU pragma: keep
#include "num_collect/constants/zero.h"  // IWYU pragma: keep
#include "num_collect/functions/legendre.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/roots/newton_raphson.h"
#include "num_collect/util/assert.h"
 
namespace num_collect::functions {
 
namespace impl {
 
template <base::concepts::real_scalar T>
class legendre_for_newton {
public:
    using variable_type = T;
 
    using jacobian_type = T;
 
    explicit legendre_for_newton(index_type degree) : degree_(degree) {
        NUM_COLLECT_PRECONDITION(
            degree >= 1, "Degree of Legendre function must be at least one.");
    }
 
    void evaluate_on(const variable_type& variable) {
        std::tie(value_, jacobian_) = legendre_with_diff(variable, degree_);
    }
 
    [[nodiscard]] auto value() const noexcept -> const variable_type& {
        return value_;
    }
 
    [[nodiscard]] auto jacobian() const noexcept -> const jacobian_type& {
        return jacobian_;
    }
 
private:
    index_type degree_;
 
    variable_type value_{};
 
    jacobian_type jacobian_{};
};
 
}  // namespace impl
 
template <base::concepts::real_scalar T>
class legendre_roots {
public:
    using variable_type = T;
 
    explicit legendre_roots(index_type degree = 0) : degree_(degree) {
        if (degree == 0) {
            roots_.resize(0);
            return;
        }
        compute(degree);
    }
 
    void compute(index_type degree) {
        NUM_COLLECT_PRECONDITION(
            degree >= 1, "Degree of Legendre function must be at least one.");
        degree_ = degree;
 
        roots_.resize(degree_);
        const index_type roots_to_solve = degree_ / 2;
 
        using function_type = impl::legendre_for_newton<variable_type>;
        using solver_type = roots::newton_raphson<function_type>;
        auto solver = solver_type(function_type(degree_));
        constexpr auto tol_last_change =
            std::numeric_limits<variable_type>::epsilon() *
            static_cast<variable_type>(1e+2);
        constexpr auto tol_value_norm =
            std::numeric_limits<variable_type>::epsilon() *
            static_cast<variable_type>(1e+2);
        solver.tol_last_change(tol_last_change);
        solver.tol_value_norm(tol_value_norm);
 
        for (index_type i = 0; i < roots_to_solve; ++i) {
            using std::cos;
            constexpr auto offset_in_num = static_cast<variable_type>(0.75);
            constexpr auto offset_in_den = static_cast<variable_type>(0.5);
            const variable_type init_var = cos(constants::pi<variable_type> *
                (static_cast<variable_type>(i) + offset_in_num) /
                (static_cast<variable_type>(degree_) + offset_in_den));
 
            solver.init(init_var);
            solver.solve();
            roots_[i] = solver.variable();
        }
        const index_type center = (degree_ - 1) / 2;
        if (degree_ % 2 == 1) {
            roots_[center] = constants::zero<variable_type>;
        }
        for (index_type i = center + 1; i < degree_; ++i) {
            roots_[i] = -roots_[degree_ - 1 - i];
        }
    }
 
    [[nodiscard]] auto degree() const noexcept -> index_type { return degree_; }
 
    [[nodiscard]] auto size() const noexcept -> index_type {
        return roots_.size();
    }
 
    [[nodiscard]] auto root(index_type i) const -> variable_type {
        NUM_COLLECT_DEBUG_ASSERT(i >= 0);
        NUM_COLLECT_DEBUG_ASSERT(i < roots_.size());
        return roots_[i];
    }
 
    [[nodiscard]] auto operator[](index_type i) const -> variable_type {
        return root(i);
    }
 
private:
    index_type degree_;
 
    Eigen::Matrix<variable_type, Eigen::Dynamic, 1> roots_{};
};
 
}  // namespace num_collect::functions