2022-11-04 03:40:36 +01:00
|
|
|
/*
|
|
|
|
|
RATIONAL.hpp - support for precise rational arithmetics
|
|
|
|
|
|
|
|
|
|
Copyright (C) Lumiera.org
|
|
|
|
|
2022, Hermann Vosseler <Ichthyostega@web.de>
|
|
|
|
|
|
|
|
|
|
This program is free software; you can redistribute it and/or
|
|
|
|
|
modify it under the terms of the GNU General Public License as
|
|
|
|
|
published by the Free Software Foundation; either version 2 of
|
|
|
|
|
the License, or (at your option) any later version.
|
|
|
|
|
|
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
|
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
|
GNU General Public License for more details.
|
|
|
|
|
|
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
|
|
|
along with this program; if not, write to the Free Software
|
|
|
|
|
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
|
|
|
|
|
|
|
|
|
*/
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/** @file rational.hpp
|
|
|
|
|
** Rational number support, based on `boost::rational`.
|
|
|
|
|
** As an extension to integral arithmetics, rational numbers can be defined
|
|
|
|
|
** as a pair (numerator, denominator); since most calculations imply multiplication
|
|
|
|
|
** by common factors, each calculation will be followed by normalisation to greatest
|
|
|
|
|
** common denominator, to keep numbers within value range. Obviously, this incurs
|
|
|
|
|
** a significant performance penalty — while on the other hand allowing for lossless
|
|
|
|
|
** computations on fractional scales, which can be notoriously difficult to handle
|
|
|
|
|
** with floating point numbers. The primary motivation for using this number format
|
|
|
|
|
** is for handling fractional time values properly, e.g 1/30 sec or 1/44100 sec.
|
|
|
|
|
**
|
|
|
|
|
** The underlying implementation from boost::rational can be parametrised with various
|
|
|
|
|
** integral data types; since our time handling is based on 64bit integers, we mainly
|
|
|
|
|
** use the specialisation `boost::rational<int64_t>`.
|
|
|
|
|
**
|
|
|
|
|
** @note all compatible integral types can be automatically converted to rational
|
|
|
|
|
** numbers, which is a lossless conversion. The opposite is not true: to get
|
|
|
|
|
** a "ordinary" number — be it integral or floating point — an explicit
|
|
|
|
|
** conversion using `rational_cast<NUM> (fraction)` is necessary, which
|
|
|
|
|
** performs the division of `numerator/denominator` in the target value domain.
|
|
|
|
|
**
|
2022-11-06 23:56:56 +01:00
|
|
|
** # Perils of fractional arithmetics
|
|
|
|
|
**
|
|
|
|
|
** While the always precise results of integral numbers might seem compelling, the
|
|
|
|
|
** danger of _numeric overflow_ is significantly increased by fractional computations.
|
|
|
|
|
** Most notably, this danger is *not limited to large numbers*. Adding two fractional
|
|
|
|
|
** number requires multiplications with both denominators, which can overflow easily.
|
|
|
|
|
** Thus, for every given fractional number, there is a class of »dangerous counterparts«
|
|
|
|
|
** which can not be added without derailing the computation, leading to arbitrary wrong
|
|
|
|
|
** results without detectable failure. And these problematic counterparts are distributed
|
|
|
|
|
** _over the whole valid numeric range._ To give an extreme example, any numbers of the
|
|
|
|
|
** form `n / INT_MAX` can not be added or subtracted with any other rational number > 1,
|
|
|
|
|
** while being themselves perfectly valid and representable.
|
|
|
|
|
** \par rule of thumb
|
|
|
|
|
** Use fractional arithmetics only where it is possible to control the denominators involved.
|
|
|
|
|
** Never use them for computations drawing from arbitrary (external) input.
|
|
|
|
|
**
|
2022-11-14 05:20:37 +01:00
|
|
|
** @see Rational_test
|
2022-11-04 03:40:36 +01:00
|
|
|
** @see zoom-window.hpp
|
|
|
|
|
** @see timevalue.hpp
|
|
|
|
|
*/
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
#ifndef LIB_RATIONAL_H
|
|
|
|
|
#define LIB_RATIONAL_H
|
|
|
|
|
|
|
|
|
|
|
2022-11-13 16:52:12 +01:00
|
|
|
#include <cmath>
|
|
|
|
|
#include <limits>
|
2022-11-04 03:40:36 +01:00
|
|
|
#include <stdint.h>
|
|
|
|
|
#include <boost/rational.hpp>
|
|
|
|
|
|
Timeline: safely calculate sum/difference of large fractional times
...in a similar vein as done for the product calculation.
In this case, we need to check the dimensions carefully and pick
the best calculation path, but as long as the overall result can
be represented, it should be possible to carry out the calculation
with fractional values, albeit introducing a small error.
As a follow-up, I have now also refactored the re-quantisation
functions, to be usable for general requantisation to another grid,
and I used these to replace the *naive* implementation of the
conversion FSecs -> µ-Grid, which caused a lot of integer-wrap-around
However, while the test now works basically without glitch or wrap,
the window position is still numerically of by 1e-6, which becomes
quite noticeably here due to the large overall span used for the test.
2022-12-01 23:23:50 +01:00
|
|
|
#include "lib/util-quant.hpp"
|
2022-11-04 03:40:36 +01:00
|
|
|
|
|
|
|
|
|
|
|
|
|
namespace util {
|
|
|
|
|
|
|
|
|
|
using Rat = boost::rational<int64_t>;
|
|
|
|
|
using boost::rational_cast;
|
2022-11-13 16:52:12 +01:00
|
|
|
using std::abs;
|
2022-11-04 03:40:36 +01:00
|
|
|
|
2022-11-14 05:20:37 +01:00
|
|
|
// these are neither standard C++ nor POSIX, yet pretty much common place...
|
|
|
|
|
using uchar = unsigned char;
|
|
|
|
|
using uint = unsigned int;
|
|
|
|
|
|
2022-11-13 16:52:12 +01:00
|
|
|
|
|
|
|
|
/**
|
|
|
|
|
* Integral binary logarithm (disregarding fractional part)
|
|
|
|
|
* @return index of the largest bit set in `num`; -1 for `num==0`
|
|
|
|
|
* @todo C++20 will provide `std::bit_width(i)` — run a microbenchmark!
|
|
|
|
|
* @remark The implementation uses an unrolled loop to break down the given number
|
|
|
|
|
* in a logarithmic search, subtracting away the larger powers of 2 first.
|
|
|
|
|
* Explained 10/2021 by user «[ToddLehman]» in this [stackoverflow].
|
2022-11-14 05:20:37 +01:00
|
|
|
* @note Microbenchmarks indicate that this function and `std::ilogb(double)` perform
|
|
|
|
|
* in the same order of magnitude (which is surprising). This function gets
|
|
|
|
|
* slightly faster for smaller data types. The naive bitshift-count implementation
|
|
|
|
|
* is always significantly slower (8 times for int64_t, 1.6 times for int8_t)
|
2022-11-15 02:13:57 +01:00
|
|
|
* @see Rational_test::verify_intLog2()
|
2022-11-13 16:52:12 +01:00
|
|
|
* @see ZoomWindow_test
|
|
|
|
|
*
|
|
|
|
|
* [ToddLehman]: https://stackoverflow.com/users/267551/todd-lehman
|
|
|
|
|
* [stackoverflow]: https://stackoverflow.com/a/24748637 "How to do an integer log2()"
|
|
|
|
|
*/
|
2022-11-14 05:20:37 +01:00
|
|
|
template<typename I>
|
2022-11-13 16:52:12 +01:00
|
|
|
inline int
|
2022-11-14 05:20:37 +01:00
|
|
|
ilog2 (I num)
|
2022-11-13 16:52:12 +01:00
|
|
|
{
|
2022-11-14 05:20:37 +01:00
|
|
|
if (num <= 0)
|
|
|
|
|
return -1;
|
|
|
|
|
const I MAX_POW = sizeof(I)*CHAR_BIT - 1;
|
2022-11-13 16:52:12 +01:00
|
|
|
int logB{0};
|
2022-11-14 05:20:37 +01:00
|
|
|
auto remove_power = [&](I pow)
|
2022-11-13 16:52:12 +01:00
|
|
|
{
|
2022-11-14 05:20:37 +01:00
|
|
|
if (pow > MAX_POW) return;
|
|
|
|
|
if (num >= I{1} << pow)
|
2022-11-13 16:52:12 +01:00
|
|
|
{
|
|
|
|
|
logB += pow;
|
|
|
|
|
num >>= pow;
|
|
|
|
|
}
|
|
|
|
|
};
|
|
|
|
|
remove_power(32);
|
|
|
|
|
remove_power(16);
|
|
|
|
|
remove_power (8);
|
|
|
|
|
remove_power (4);
|
|
|
|
|
remove_power (2);
|
|
|
|
|
remove_power (1);
|
|
|
|
|
|
|
|
|
|
return logB;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
inline bool
|
|
|
|
|
can_represent_Product (int64_t a, int64_t b)
|
|
|
|
|
{
|
2022-11-14 05:20:37 +01:00
|
|
|
return ilog2(abs(a))+1
|
|
|
|
|
+ ilog2(abs(b))+1
|
2022-11-13 16:52:12 +01:00
|
|
|
< 63;
|
|
|
|
|
}
|
|
|
|
|
|
2022-11-29 02:00:41 +01:00
|
|
|
inline bool
|
|
|
|
|
can_represent_Product (Rat a, Rat b)
|
|
|
|
|
{
|
|
|
|
|
return can_represent_Product(a.numerator(), b.numerator())
|
|
|
|
|
and can_represent_Product(a.denominator(), b.denominator());
|
|
|
|
|
}
|
|
|
|
|
|
2022-11-13 16:52:12 +01:00
|
|
|
inline bool
|
|
|
|
|
can_represent_Sum (Rat a, Rat b)
|
|
|
|
|
{
|
|
|
|
|
return can_represent_Product(a.numerator(), b.denominator())
|
|
|
|
|
and can_represent_Product(b.numerator(), a.denominator());
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
Timeline: safely calculate sum/difference of large fractional times
...in a similar vein as done for the product calculation.
In this case, we need to check the dimensions carefully and pick
the best calculation path, but as long as the overall result can
be represented, it should be possible to carry out the calculation
with fractional values, albeit introducing a small error.
As a follow-up, I have now also refactored the re-quantisation
functions, to be usable for general requantisation to another grid,
and I used these to replace the *naive* implementation of the
conversion FSecs -> µ-Grid, which caused a lot of integer-wrap-around
However, while the test now works basically without glitch or wrap,
the window position is still numerically of by 1e-6, which becomes
quite noticeably here due to the large overall span used for the test.
2022-12-01 23:23:50 +01:00
|
|
|
/**
|
|
|
|
|
* Re-Quantise a number into a new grid, truncating to the next lower grid point.
|
|
|
|
|
* @remark Grid-aligned values can be interpreted as rational numbers (integer fractions),
|
|
|
|
|
* where the quantiser corresponds to the denominator and the numerator counts
|
|
|
|
|
* the number of grid steps. To work both around precision problems and the
|
|
|
|
|
* danger of integer wrap-around, the integer division is performed on the
|
|
|
|
|
* old value and then the re-quantisation done on the remainder, using
|
|
|
|
|
* 128bit floating point for maximum precision. This operation can
|
|
|
|
|
* also be used to re-form a fraction to be cast in terms of the
|
|
|
|
|
* new quantiser; this introduces a tiny error, but typically
|
|
|
|
|
* allows for safe or simplified calculations.
|
|
|
|
|
* @param num the count in old grid-steps (#den) or the numerator
|
|
|
|
|
* @param den the old quantiser or the denominator of a fraction
|
|
|
|
|
* @param u the new quantiser or the new denominator to use
|
|
|
|
|
* @return the adjusted numerator, so that the fraction with u
|
|
|
|
|
* will be almost the same than dividing #num by #den
|
|
|
|
|
*/
|
|
|
|
|
inline int64_t
|
|
|
|
|
reQuant (int64_t num, int64_t den, int64_t u)
|
|
|
|
|
{
|
|
|
|
|
u = 0!=u? u:1;
|
|
|
|
|
auto [d,r] = util::iDiv (num, den);
|
|
|
|
|
using f128 = long double;
|
|
|
|
|
// round to smallest integer fraction, to shake off "number dust"
|
|
|
|
|
f128 const ROUND_ULP = 1 + 1/(f128(std::numeric_limits<int64_t>::max()) * 2);
|
|
|
|
|
|
|
|
|
|
// construct approximation quantised to 1/u
|
|
|
|
|
f128 frac = f128(r) / den;
|
|
|
|
|
int64_t res = d*u + int64_t(frac*u * ROUND_ULP);
|
Library: rectify clipping of time::Duration (see #1263)
This is a deep refactoring to allow to represent the distance
between all valid time points as a time::Offset or time::Duration.
By design this is possible, since Time::MAX was defined as 1/30 of
the maximum value technically representable as int64_t. However,
introducing a different limiter for offsets and durations turns
out difficult, due to the inconsistencies in the exiting hierarchy
of temporal entities. Which in turn seems to stem from the unfortunate
decision to make time entities immutable, see #1261
Since the limiter is hard wired into the `time::TimeValue` constructor,
we are forced to create a "backdoor" of sorts, to pass up values
with different limiting from child classes. This would not be so
much of a problem if calculations weren't forced to go through `TimeVar`,
which does not distinguish between time points and time durations.
This solution rearranges all checks to be performed now by time::Offset,
while time::Duration will only take the absolute value at construction,
based on the fact that there is no valid construction path to yield
a duration which does not go through an offset first.
Later, when we're ready to sort out the implementation base of time values
(see #1258), this design issue should be revisited
- either we'll allow derived classes explicitly to invoke the limiter functions
- or we may be able to have an automatic conversion path from clearly
marked base implementation types, in which case we wouldn't use the
buildRaw_() and _raw() "backdoor" functions any more...
2022-12-05 00:58:32 +01:00
|
|
|
ENSURE (abs (f128(res)/u - rational_cast<f128>(Rat{num,den})) <= 1.0/abs(u)
|
2022-12-10 04:26:22 +01:00
|
|
|
,"Requantisation error exceeded num=%li / den=%li -> res=%li / quant=%li"
|
Library: rectify clipping of time::Duration (see #1263)
This is a deep refactoring to allow to represent the distance
between all valid time points as a time::Offset or time::Duration.
By design this is possible, since Time::MAX was defined as 1/30 of
the maximum value technically representable as int64_t. However,
introducing a different limiter for offsets and durations turns
out difficult, due to the inconsistencies in the exiting hierarchy
of temporal entities. Which in turn seems to stem from the unfortunate
decision to make time entities immutable, see #1261
Since the limiter is hard wired into the `time::TimeValue` constructor,
we are forced to create a "backdoor" of sorts, to pass up values
with different limiting from child classes. This would not be so
much of a problem if calculations weren't forced to go through `TimeVar`,
which does not distinguish between time points and time durations.
This solution rearranges all checks to be performed now by time::Offset,
while time::Duration will only take the absolute value at construction,
based on the fact that there is no valid construction path to yield
a duration which does not go through an offset first.
Later, when we're ready to sort out the implementation base of time values
(see #1258), this design issue should be revisited
- either we'll allow derived classes explicitly to invoke the limiter functions
- or we may be able to have an automatic conversion path from clearly
marked base implementation types, in which case we wouldn't use the
buildRaw_() and _raw() "backdoor" functions any more...
2022-12-05 00:58:32 +01:00
|
|
|
, num, den, res, u);
|
Timeline: safely calculate sum/difference of large fractional times
...in a similar vein as done for the product calculation.
In this case, we need to check the dimensions carefully and pick
the best calculation path, but as long as the overall result can
be represented, it should be possible to carry out the calculation
with fractional values, albeit introducing a small error.
As a follow-up, I have now also refactored the re-quantisation
functions, to be usable for general requantisation to another grid,
and I used these to replace the *naive* implementation of the
conversion FSecs -> µ-Grid, which caused a lot of integer-wrap-around
However, while the test now works basically without glitch or wrap,
the window position is still numerically of by 1e-6, which becomes
quite noticeably here due to the large overall span used for the test.
2022-12-01 23:23:50 +01:00
|
|
|
return res;
|
|
|
|
|
}
|
|
|
|
|
|
2022-11-13 16:52:12 +01:00
|
|
|
/**
|
|
|
|
|
* re-Quantise a rational number to a (typically smaller) denominator.
|
|
|
|
|
* @param u the new denominator to use
|
|
|
|
|
* @warning this is a lossy operation and possibly introduces an error
|
|
|
|
|
* of up to 1/u
|
|
|
|
|
* @remark Rational numbers with large numerators can be »poisonous«,
|
|
|
|
|
* causing numeric overflow when used, even just additively.
|
|
|
|
|
* This function can thus be used to _"sanitise"_ a number,
|
|
|
|
|
* and thus accept a small error while preventing overflow.
|
2022-11-15 02:13:57 +01:00
|
|
|
* @note using extended-precision floating point to get close to the
|
|
|
|
|
* quantiser even for large denominator. Some platforms
|
|
|
|
|
* fall back to double, leading to extended error bounds
|
2022-11-13 16:52:12 +01:00
|
|
|
*/
|
|
|
|
|
inline Rat
|
|
|
|
|
reQuant (Rat src, int64_t u)
|
|
|
|
|
{
|
Timeline: safely calculate sum/difference of large fractional times
...in a similar vein as done for the product calculation.
In this case, we need to check the dimensions carefully and pick
the best calculation path, but as long as the overall result can
be represented, it should be possible to carry out the calculation
with fractional values, albeit introducing a small error.
As a follow-up, I have now also refactored the re-quantisation
functions, to be usable for general requantisation to another grid,
and I used these to replace the *naive* implementation of the
conversion FSecs -> µ-Grid, which caused a lot of integer-wrap-around
However, while the test now works basically without glitch or wrap,
the window position is still numerically of by 1e-6, which becomes
quite noticeably here due to the large overall span used for the test.
2022-12-01 23:23:50 +01:00
|
|
|
return Rat{reQuant (src.numerator(), src.denominator(), u), u};
|
2022-11-13 16:52:12 +01:00
|
|
|
}
|
2022-11-04 03:40:36 +01:00
|
|
|
} // namespace util
|
|
|
|
|
|
|
|
|
|
|
2022-11-13 16:52:12 +01:00
|
|
|
|
2022-11-06 23:56:56 +01:00
|
|
|
/**
|
2022-11-04 03:40:36 +01:00
|
|
|
* user defined literal for constant rational numbers.
|
|
|
|
|
* \code
|
|
|
|
|
* Rat twoThirds = 2_r/3;
|
2022-11-06 23:56:56 +01:00
|
|
|
* \endcode
|
2022-11-04 03:40:36 +01:00
|
|
|
*/
|
|
|
|
|
inline util::Rat
|
|
|
|
|
operator""_r (unsigned long long num)
|
|
|
|
|
{
|
|
|
|
|
return util::Rat{num};
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
#endif /*LIB_RATIONAL_H*/
|