quad.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 <tuple>
#include <type_traits>
 
#include "num_collect/multi_double/impl/basic_operations.h"
 
namespace num_collect::multi_double {
 
class quad {
public:
    constexpr quad() noexcept = default;
 
    constexpr quad(double high, double low) noexcept : high_(high), low_(low) {}
 
    template <typename Scalar,
        std::enable_if_t<(std::is_integral_v<Scalar> ||
                             std::is_floating_point_v<Scalar>) &&
            (!std::is_same_v<Scalar,
                quad>)>* = nullptr>
    constexpr quad(Scalar value) noexcept  // NOLINT
        : high_(static_cast<double>(value)) {}
 
    [[nodiscard]] auto high() const noexcept -> double { return high_; }
 
    [[nodiscard]] auto low() const noexcept -> double { return low_; }
 
    auto operator-() const noexcept -> quad { return quad(-high_, -low_); }
 
    auto operator+=(const quad& right) noexcept -> quad& {
        auto [x_h, x_l] = impl::two_sum(high_, right.high_);
        x_l += low_;
        x_l += right.low_;
        std::tie(high_, low_) = impl::quick_two_sum(x_h, x_l);
        return *this;
    }
 
    auto operator-=(const quad& right) noexcept -> quad& {
        return operator+=(-right);
    }
 
    auto operator*=(const quad& right) noexcept -> quad& {
        auto [x_h, x_l] = impl::two_prod(high_, right.high_);
        x_l += high_ * right.low_;
        x_l += low_ * right.high_;
        std::tie(high_, low_) = impl::quick_two_sum(x_h, x_l);
        return *this;
    }
 
    auto operator/=(const quad& right) noexcept -> quad& {
        const double inv_right_h = 1.0 / right.high_;
        const double rate_right = right.low_ * inv_right_h;
        const double x_h = high_ * inv_right_h;
        const auto [r_1, r_2] = impl::two_prod(x_h, right.high_);
        double x_l = ((high_ - r_1) - r_2) * inv_right_h;
        x_l += x_h * ((low_ / high_) - rate_right);
        std::tie(high_, low_) = impl::quick_two_sum(x_h, x_l);
        return *this;
    }
 
private:
    double high_{0.0};
 
    double low_{0.0};
};
 
inline auto operator+(const quad& left, const quad& right) -> quad {
    return quad(left) += right;
}
 
inline auto operator-(const quad& left, const quad& right) -> quad {
    return quad(left) -= right;
}
 
inline auto operator*(const quad& left, const quad& right) -> quad {
    return quad(left) *= right;
}
 
inline auto operator/(const quad& left, const quad& right) -> quad {
    return quad(left) /= right;
}
 
}  // namespace num_collect::multi_double