LUMIERA.clone/src/lib/util-quant.hpp
Ichthyostega 289f92da7e 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

159 lines
4.5 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 <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
*/
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
}
} // namespace util
#endif /*UTIL_QUANT_H*/