mercredi 30 mars 2016

rational approximation of square root of std::ratio at compile-time

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

  1. 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)

  1. 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