Ichthyostega
597d8191c7
...turns out challenging, since our intention here is borderline to the intended design of the Lumiera ETD. It ''should work'' though, when combined with a Variant-visitor...
207 lines
6.4 KiB
C++
207 lines
6.4 KiB
C++
/*
|
|
UTIL-QUANT.hpp - helper functions to deal with quantisation and comparison
|
|
|
|
Copyright (C) Lumiera.org
|
|
2011, 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 util-quant.hpp
|
|
** Utilities for quantisation (grid alignment) and comparisons.
|
|
*/
|
|
|
|
|
|
#ifndef LIB_UTIL_QUANT_H
|
|
#define LIB_UTIL_QUANT_H
|
|
|
|
#include <cstdlib>
|
|
#include <climits>
|
|
#include <cfloat>
|
|
#include <cmath>
|
|
|
|
|
|
|
|
namespace util {
|
|
|
|
|
|
/** helper to treat int or long division uniformly */
|
|
template<typename I>
|
|
struct IDiv
|
|
{
|
|
I quot;
|
|
I rem;
|
|
|
|
IDiv (I num, I den)
|
|
: quot(num/den)
|
|
, rem(num - quot*den)
|
|
{ }
|
|
};
|
|
|
|
template<>
|
|
struct IDiv<int>
|
|
: div_t
|
|
{
|
|
IDiv<int> (int num, int den)
|
|
: div_t(div (num,den))
|
|
{ }
|
|
};
|
|
|
|
template<>
|
|
struct IDiv<long>
|
|
: ldiv_t
|
|
{
|
|
IDiv<long> (long num, long den)
|
|
: ldiv_t(ldiv (num,den))
|
|
{ }
|
|
};
|
|
|
|
template<>
|
|
struct IDiv<long long>
|
|
: lldiv_t
|
|
{
|
|
IDiv<long long> (long long num, long long den)
|
|
: lldiv_t(lldiv (num,den))
|
|
{ }
|
|
};
|
|
|
|
template<typename I>
|
|
inline IDiv<I>
|
|
iDiv (I num, I den) ///< support type inference and auto typing...
|
|
{
|
|
return IDiv<I>{num,den};
|
|
}
|
|
|
|
|
|
/** floor function for integer arithmetics.
|
|
* Unlike the built-in integer division, this function
|
|
* always rounds towards the _next smaller integer,_
|
|
* even for negative numbers.
|
|
* @warning floor on doubles performs way better
|
|
* @see UtilFloordiv_test
|
|
*/
|
|
template<typename I>
|
|
inline I
|
|
floordiv (I num, I den)
|
|
{
|
|
if (0 < (num^den))
|
|
return num/den;
|
|
else
|
|
{ // truncate similar to floor()
|
|
IDiv<I> res(num,den);
|
|
return (res.rem)? res.quot-1 // negative results truncated towards next smaller int
|
|
: res.quot; //..unless the division result not truncated at all
|
|
}
|
|
}
|
|
|
|
/** scale wrapping operation.
|
|
* Quantises the numerator value into the scale given by the denominator.
|
|
* Unlike the built-in integer division, this function always rounds towards
|
|
* the _next smaller integer_ and also relates the remainder (=modulo) to
|
|
* this next lower scale grid point.
|
|
* @return quotient and remainder packed into a struct
|
|
* @see UtilFloorwarp_test
|
|
*/
|
|
template<typename I>
|
|
inline IDiv<I>
|
|
floorwrap (I num, I den)
|
|
{
|
|
IDiv<I> res(num,den);
|
|
if (0 > (num^den) && res.rem)
|
|
{ // negative results
|
|
// wrapped similar to floor()
|
|
--res.quot;
|
|
res.rem = den - (-res.rem);
|
|
}
|
|
return res;
|
|
}
|
|
|
|
|
|
|
|
/**
|
|
* epsilon comparison of doubles.
|
|
* @remarks Floating point calculations are only accurate up to a certain degree,
|
|
* and we need to adjust for the magnitude of the involved numbers, since
|
|
* floating point numbers are scaled by the exponent. Moreover, we need
|
|
* to be careful with very small numbers (close to zero), where calculating
|
|
* the difference could yield coarse grained 'subnormal' values.
|
|
* @param ulp number of grid steps to allow for difference (default = 2).
|
|
* Here, a 'grid step' is the smallest difference to 1.0 which can be
|
|
* represented in floating point ('units in the last place')
|
|
* @warning don't use this for comparison against zero, rather use an absolute epsilon then.
|
|
* @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
|
|
* @see http://en.cppreference.com/w/cpp/types/numeric_limits/epsilon
|
|
* @see https://en.wikipedia.org/wiki/Unit_in_the_last_place
|
|
* @todo 3/2024 seems we have solved this problem several times meanwhile /////////////////////////////////TICKET #1360 sort out floating-point rounding and precision
|
|
*/
|
|
inline bool
|
|
almostEqual (double d1, double d2, unsigned int ulp =2)
|
|
{
|
|
using std::fabs;
|
|
return fabs (d1-d2) < DBL_EPSILON * fabs (d1+d2) * ulp
|
|
|| fabs (d1-d2) < DBL_MIN; // special treatment for subnormal results
|
|
}
|
|
|
|
|
|
|
|
/**
|
|
* 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].
|
|
* @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)
|
|
* @see Rational_test::verify_intLog2()
|
|
* @see ZoomWindow_test
|
|
*
|
|
* [ToddLehman]: https://stackoverflow.com/users/267551/todd-lehman
|
|
* [stackoverflow]: https://stackoverflow.com/a/24748637 "How to do an integer log2()"
|
|
*/
|
|
template<typename I>
|
|
inline constexpr int
|
|
ilog2 (I num)
|
|
{
|
|
if (num <= 0)
|
|
return -1;
|
|
const I MAX_POW = sizeof(I)*CHAR_BIT - 1;
|
|
int logB{0};
|
|
auto remove_power = [&](I pow) constexpr
|
|
{
|
|
if (pow > MAX_POW) return;
|
|
if (num >= I{1} << pow)
|
|
{
|
|
logB += pow;
|
|
num >>= pow;
|
|
}
|
|
};
|
|
remove_power(32);
|
|
remove_power(16);
|
|
remove_power (8);
|
|
remove_power (4);
|
|
remove_power (2);
|
|
remove_power (1);
|
|
|
|
return logB;
|
|
}
|
|
|
|
|
|
} // namespace util
|
|
#endif /*UTIL_QUANT_H*/
|