mirror of
https://github.com/saymrwulf/uhd.git
synced 2026-05-16 21:10:10 +00:00
ee520a9df5
Note: template_lvbitx.{cpp,hpp} need to be excluded from the
list of files that clang-format gets applied against.
host/lib/dep is also excluded from this change.
101 lines
3.2 KiB
C++
101 lines
3.2 KiB
C++
//
|
|
// Copyright 2018 Ettus Research, a National Instruments Company
|
|
//
|
|
// SPDX-License-Identifier: GPL-3.0-or-later
|
|
//
|
|
|
|
// More math, but not meant for public API
|
|
|
|
#ifndef INCLUDED_UHDLIB_UTILS_MATH_HPP
|
|
#define INCLUDED_UHDLIB_UTILS_MATH_HPP
|
|
|
|
#include <uhdlib/utils/narrow.hpp>
|
|
#include <cmath>
|
|
#include <vector>
|
|
|
|
namespace uhd { namespace math {
|
|
|
|
/*! log2(num), rounded up to the nearest integer.
|
|
*/
|
|
template <class T>
|
|
T ceil_log2(T num)
|
|
{
|
|
return std::ceil(std::log(num) / std::log(T(2)));
|
|
}
|
|
|
|
/**
|
|
* Function which attempts to find integer values a and b such that
|
|
* a / b approximates the decimal represented by f within max_error and
|
|
* b is not greater than maximum_denominator.
|
|
*
|
|
* If the approximation cannot achieve the desired error without exceeding
|
|
* the maximum denominator, b is set to the maximum value and a is set to
|
|
* the closest value.
|
|
*
|
|
* @param f is a positive decimal to be converted, must be between 0 and 1
|
|
* @param maximum_denominator maximum value allowed for b
|
|
* @param max_error how close to f the expression a / b should be
|
|
*/
|
|
template <typename IntegerType>
|
|
std::pair<IntegerType, IntegerType> rational_approximation(
|
|
const double f, const IntegerType maximum_denominator, const double max_error)
|
|
{
|
|
static constexpr IntegerType MIN_DENOM = 1;
|
|
static constexpr size_t MAX_APPROXIMATIONS = 64;
|
|
|
|
UHD_ASSERT_THROW(maximum_denominator <= std::numeric_limits<IntegerType>::max());
|
|
UHD_ASSERT_THROW(f < 1 and f >= 0);
|
|
|
|
// This function uses a successive approximations formula to attempt to
|
|
// find a "best" rational to use to represent a decimal. This algorithm
|
|
// finds a continued fraction of up to 64 terms, or such that the last
|
|
// term is less than max_error
|
|
|
|
if (f < max_error) {
|
|
return {0, MIN_DENOM};
|
|
}
|
|
|
|
double c = f;
|
|
std::vector<double> saved_denoms = {c};
|
|
|
|
// Create the continued fraction by taking the reciprocal of the
|
|
// fractional part, expressing the denominator as a mixed number,
|
|
// then repeating the algorithm on the fractional part of that mixed
|
|
// number until a maximum number of terms or the fractional part is
|
|
// nearly zero.
|
|
for (size_t i = 0; i < MAX_APPROXIMATIONS; ++i) {
|
|
double x = std::floor(1.0 / c);
|
|
c = (1.0 / c) - x;
|
|
|
|
saved_denoms.push_back(x);
|
|
|
|
if (std::abs(c) < max_error)
|
|
break;
|
|
}
|
|
|
|
double num = 1.0;
|
|
double denom = saved_denoms.back();
|
|
|
|
// Calculate a single rational which will be equivalent to the
|
|
// continued fraction created earlier. Because the continued fraction
|
|
// is composed of only integers, the final rational will be as well.
|
|
for (auto it = saved_denoms.rbegin() + 1; it != saved_denoms.rend() - 1; ++it) {
|
|
double new_denom = denom * (*it) + num;
|
|
if (new_denom > maximum_denominator) {
|
|
// We can't do any better than using the maximum denominator
|
|
num = std::round(f * maximum_denominator);
|
|
denom = maximum_denominator;
|
|
break;
|
|
}
|
|
|
|
num = denom;
|
|
denom = new_denom;
|
|
}
|
|
|
|
return {uhd::narrow<IntegerType>(num), uhd::narrow<IntegerType>(denom)};
|
|
}
|
|
|
|
|
|
}} /* namespace uhd::math */
|
|
|
|
#endif /* INCLUDED_UHDLIB_UTILS_MATH_HPP */
|