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

#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