1 //
2 // Copyright 2018 Ettus Research, a National Instruments Company
3 //
4 // SPDX-License-Identifier: GPL-3.0-or-later
5 //
6
7 // More math, but not meant for public API
8
9 #pragma once
10
11 #include <uhdlib/utils/narrow.hpp>
12 #include <cmath>
13 #include <vector>
14
15 namespace uhd { namespace math {
16
17 /*! log2(num), rounded up to the nearest integer.
18 */
19 template <class T>
ceil_log2(T num)20 T ceil_log2(T num)
21 {
22 return std::ceil(std::log(num) / std::log(T(2)));
23 }
24
25 /**
26 * Function which attempts to find integer values a and b such that
27 * a / b approximates the decimal represented by f within max_error and
28 * b is not greater than maximum_denominator.
29 *
30 * If the approximation cannot achieve the desired error without exceeding
31 * the maximum denominator, b is set to the maximum value and a is set to
32 * the closest value.
33 *
34 * @param f is a positive decimal to be converted, must be between 0 and 1
35 * @param maximum_denominator maximum value allowed for b
36 * @param max_error how close to f the expression a / b should be
37 */
38 template <typename IntegerType>
rational_approximation(const double f,const IntegerType maximum_denominator,const double max_error)39 std::pair<IntegerType, IntegerType> rational_approximation(
40 const double f, const IntegerType maximum_denominator, const double max_error)
41 {
42 static constexpr IntegerType MIN_DENOM = 1;
43 static constexpr size_t MAX_APPROXIMATIONS = 64;
44
45 UHD_ASSERT_THROW(maximum_denominator <= std::numeric_limits<IntegerType>::max());
46 UHD_ASSERT_THROW(f < 1 and f >= 0);
47
48 // This function uses a successive approximations formula to attempt to
49 // find a "best" rational to use to represent a decimal. This algorithm
50 // finds a continued fraction of up to 64 terms, or such that the last
51 // term is less than max_error
52
53 if (f < max_error) {
54 return {0, MIN_DENOM};
55 }
56
57 double c = f;
58 std::vector<double> saved_denoms = {c};
59
60 // Create the continued fraction by taking the reciprocal of the
61 // fractional part, expressing the denominator as a mixed number,
62 // then repeating the algorithm on the fractional part of that mixed
63 // number until a maximum number of terms or the fractional part is
64 // nearly zero.
65 for (size_t i = 0; i < MAX_APPROXIMATIONS; ++i) {
66 double x = std::floor(1.0 / c);
67 c = (1.0 / c) - x;
68
69 saved_denoms.push_back(x);
70
71 if (std::abs(c) < max_error)
72 break;
73 }
74
75 double num = 1.0;
76 double denom = saved_denoms.back();
77
78 // Calculate a single rational which will be equivalent to the
79 // continued fraction created earlier. Because the continued fraction
80 // is composed of only integers, the final rational will be as well.
81 for (auto it = saved_denoms.rbegin() + 1; it != saved_denoms.rend() - 1; ++it) {
82 double new_denom = denom * (*it) + num;
83 if (new_denom > maximum_denominator) {
84 // We can't do any better than using the maximum denominator
85 num = std::round(f * maximum_denominator);
86 denom = maximum_denominator;
87 break;
88 }
89
90 num = denom;
91 denom = new_denom;
92 }
93
94 return {uhd::narrow<IntegerType>(num), uhd::narrow<IntegerType>(denom)};
95 }
96
97
98 }} /* namespace uhd::math */
99