I'm trying to find a rational approximation of the square root of a std::ratio
at compile time. It would be extremely useful to derive ellipsoid parameters from defined parameters for coordinate conversions, which are themselves defined as std::ratio
.
There is a question about finding powers/roots of std::ratio, but as a condition of that question, it was OK to fail if the ratio had no integral root, which is the opposite of what I want. I would instead like to find the closest reasonable approximation I can.
I've come up with the following meta-program that calculates the roots based on the Newton-Raphson Method, which is known to produce (relatively) accurate results with only a few iterations:
namespace detail
{
/// implementation of ratio_sqrt
template<class N, class K = std::ratio<4>, std::intmax_t RecursionDepth = 5>
struct ratio_sqrt_impl
{
// Recursive Newton-Raphson
// EQUATION: K_{n+1} = (K_{n} - N / K_{n}) / 2
// WHERE:
// K_{n+1} : square root approximation
// K_{n} : previous square root approximation
// N : ratio whose square root we are finding
using type = typename ratio_sqrt_impl<N, std::ratio_subtract<K,
std::ratio_divide<std::ratio_subtract<std::ratio_multiply<K, K>, N>,
std::ratio_multiply<std::ratio<2>, K>>>, RecursionDepth - 1>::type;
};
template<class N, class K>
struct ratio_sqrt_impl<N, K, 1>
{
using type = K;
};
}
template<class Ratio>
using ratio_sqrt = typename detail::ratio_sqrt_impl<Ratio>::type;
With the example usage:
// Error calculations
using rt2 = ratio_sqrt<std::ratio<2>>;
std::cout << (sqrt(2) - ((double)rt2::num / rt2::den))/sqrt(2) << std::endl;
scalar_t result = pow<2>(scalar_t((double)rt2::num / rt2::den));
std::cout << (2 - result.toDouble()) / 2 << std::endl;
using rt4 = ratio_sqrt<std::ratio<4>>;
std::cout << (sqrt(4) - ((double)rt4::num / rt4::den)) / sqrt(4) << std::endl;
using rt10 = ratio_sqrt<std::ratio<10>>;
std::cout << (sqrt(10) - ((double)rt10::num / rt10::den)) / sqrt(10) << std::endl;
Producing the results:
1.46538e-05 // sqrt(2)
4.64611e-08 // sqrt(4)
2.38737e-15 // sqrt(10)
which could certainly be decent for some applications.
The Problems
- The biggest problem here is the fixed Recursion depth. These ratios get BIG, very quickly, and so for roots > 100, this overflows like crazy. However, too small of a recursion depth, and you lose all the accuracy.
Is there a good way that the recursion could be adapted to the overflow depth limit, and then have the type set to be an iteration or two before that? (I say a few iterations because it might be nice to keep headroom in the integer sizes to do further calculations later)
- The initial condition of 4 seemed to be pretty magical in terms of producing the lowest error for roots < 100, but is there a more methodical way to set that?
Aucun commentaire:
Enregistrer un commentaire