You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
85 lines
2.1 KiB
85 lines
2.1 KiB
#ifndef ALGOIM_UTILITY_HPP
|
|
#define ALGOIM_UTILITY_HPP
|
|
|
|
// Minor utility methods used throughout Algoim
|
|
|
|
namespace algoim::util
|
|
{
|
|
static_assert(std::is_same_v<real, double>, "Warning: pi constant may require redefining when real != double");
|
|
static constexpr double pi = 3.141592653589793238462643383279502884197169399375105820974944592307816;
|
|
|
|
// square of u
|
|
template<typename T>
|
|
constexpr auto sqr(T u)
|
|
{
|
|
return u*u;
|
|
}
|
|
|
|
// cube of u
|
|
template<typename T>
|
|
constexpr auto cube(T u)
|
|
{
|
|
return u*u*u;
|
|
}
|
|
|
|
// sign of u:
|
|
// -1, if u < 0,
|
|
// +1, if u > 0,
|
|
// 0, otherwise
|
|
template <typename T>
|
|
constexpr int sign(T u) noexcept
|
|
{
|
|
return (T(0) < u) - (u < T(0));
|
|
}
|
|
|
|
// collapse an N-dimensional multi-index into a scalar integer according to its location
|
|
// in an N-dimensional grid of the given extent, such that the lowest dimension iterates
|
|
// the slowest, and the highest dimension corresponds to the inner-most loop, consistent
|
|
// with MultiLoop semantics. For example,
|
|
// furl( {i,j}, {m,n} ) = i*n + j
|
|
// furl( {i,j,k}, {m,n,o} ) = i*n*o + j*o + k
|
|
template<int N>
|
|
int furl(const uvector<int,N>& i, const uvector<int,N>& ext)
|
|
{
|
|
int ind = i(0);
|
|
for (int j = 1; j < N; ++j)
|
|
ind = ext(j) * ind + i(j);
|
|
return ind;
|
|
}
|
|
|
|
// compute a Givens rotation
|
|
template<typename T>
|
|
void givens_get(const T& a, const T& b, T& c, T& s)
|
|
{
|
|
using std::abs;
|
|
using std::sqrt;
|
|
if (b == 0.0)
|
|
{
|
|
c = 1.0;
|
|
s = 0.0;
|
|
}
|
|
else if (abs(b) > abs(a))
|
|
{
|
|
T tmp = a / b;
|
|
s = T(1) / sqrt(1.0 + tmp*tmp);
|
|
c = tmp * s;
|
|
}
|
|
else
|
|
{
|
|
T tmp = b / a;
|
|
c = T(1) / sqrt(1.0 + tmp*tmp);
|
|
s = tmp * c;
|
|
}
|
|
};
|
|
|
|
// apply a Givens rotation
|
|
template<typename T>
|
|
void givens_rotate(T& x, T& y, T c, T s)
|
|
{
|
|
T a = x, b = y;
|
|
x = c * a + s * b;
|
|
y = -s * a + c * b;
|
|
};
|
|
} // namespace algoim::util
|
|
|
|
#endif
|
|
|