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

Question

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
    // N is an std::ratio, not an int
    template<class N , class K = std::ratio<4>, std::intmax_t RecursionDepth = 5>
    struct ratio_sqrt_impl
    {
        static_assert(/* N is std::ratio */);

        // 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)

2. 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?

EDIT:

I'm not looking for any solutions with constexpr, as the compilers I have to support don't uniformly have it.

The problem with increasing the recursion depth is that the num/denom of the std::ratio overflow after only a couple of recursions. The accuracy of the represented square root is actually OK, but I need to find a generic solution that limits the recursion depth to the point where the ratios don't overflow (and thus, don't compile). E.g. ratio_sqrt<std::ratio<2>> can go to depth 5 before overflowing, but ratio_sqrt<std::ratio<1000>> is limited to 4.

Answer 1

Update

Fixed overflow with super-large inputs in iSqrt

test_dump's input pack no longer needs a terminating bool.

This code calculates the rational square root of a std::ratio It works with Visual Studio 2013 and g++ at IdeOne.

This uses a mediant search to converge on the target value. The search passes through the same convergents as a continued fraction, but with a few more iterations. This is similar to walking the Stern-Brocot tree, where each node is a closer rational approximation to some real target value.

This keeps converging until either the numerator or denominator of the last mediant passes a limit, close to the point of triggering an integer overflow.

At the bottom, a few test cases are calculated.

#include <iostream>
#include <iomanip>
#include <ratio>
#include <cmath>  // for sqrt(double) -- only used for testing

// old versions of std::numeric_limits can't be used at compile-time
// integer_limits is a replacement, for two's complement integer types
template <typename T>
struct integer_limits {
    static const bool is_integer = T(1)/T(3) == T(0);
    static_assert(is_integer, "integer_limits got non-integer type");
    static const bool is_signed = T(1)-T(3) < T(1);
    static const int  bit_count = sizeof(T) * 8; // 8-bit char safely assumed
    static const T    min_value = is_signed ? T(1) << (bit_count-1) : T(0);
    static const T    max_value = ~min_value;
};

// iLog2: static intmax_t log2, used by iSqrt to calculate an initial guess
template <intmax_t n>
struct iLog2 {
    template <intmax_t b, intmax_t c=0>
    struct L2 {
        static const intmax_t value = L2<b/2, c+1>::value;
    };
    template <intmax_t c>
    struct L2<1,c> {
        static const intmax_t value = c;
    };
    static_assert(n > 0, "Log2 got n <= 0");
    static const intmax_t value = L2<n>::value;
};

// iSqrt: static intmax_t square root
// starts off using iLog2 to get an approximate result range
// uses standard binary search to narrow the range and converge
// on the square root's floor value
template <intmax_t n>
struct iSqrt {
    static_assert(n >= 0, "iSqrt got a negative n");

    // main recursive binary search loop
    template <intmax_t low, intmax_t high, bool end=(low >= high)>
    struct srch {
        static const intmax_t mid       = low + (high - low)/2;
        static const intmax_t ndivm     = n / mid;
        static const int      mid_cmp   = int(mid > ndivm) - int(mid < ndivm);
        static const intmax_t next_high = mid_cmp < 0 ? high - (mid == low)  : mid;
        static const intmax_t next_low  = mid_cmp > 0 ? low                  : mid;
        static const intmax_t value     = srch<next_low, next_high>::value;
    };
    // terminating specialization for the binary search (end is true)
    template <intmax_t low, intmax_t high>
    struct srch<low, high, true> {
        static_assert(low*low <= n && (low+1)*(low+1) > n, "bad result");
        static const intmax_t value = low;
    };
    // initial range and kickoff
    static const intmax_t low_approx = intmax_t(1) << iLog2<n>::value / 2;
    static const intmax_t high_approx = low_approx << 1;
    static const intmax_t value = srch<low_approx, high_approx>::value;
};
// Specializations handling sqrt(0) and sqrt(1)
template <> struct iSqrt<0> { static const intmax_t value = 0; };
template <> struct iSqrt<1> { static const intmax_t value = 1; };

// std::ratio mediant calculation, used in the mediant search
template <typename a, typename b>
struct mediant {
    typedef std::ratio<a::num + b::num, a::den + b::den> type;
};

// conditional type selection utility
template <typename A, typename B, bool sel_a> struct select_t { typedef A type; };
template <typename A, typename B> struct select_t<A, B, false>{ typedef B type; };

// Sqrt calculates the square root of a std::ratio
template <typename r>
struct Sqrt {
    // log_limit limits how big (in bits) the intmax_t can get before ending the mediant search
    // the fraction here comes from trial and error.  Some implementations of
    // std::ratio might be more likely to overflow than others, so you can experiement with
    // this to try to get more precision.
    // This search tests convergents against the input by squaring them, so there is
    // a hard limit at 50% (half the bits).  My std::ratio implementation also
    // seems to need a lot of cushioning when it does comparisons.
    static const int log_limit = (integer_limits<intmax_t>::bit_count-1) * 31/100;
    static const intmax_t root_limit = intmax_t(1) << log_limit;

    // The main mediant search loop.  This is similar to walking the Stern-Brocot tree
    // toward a real target.  Each iteration brings the ratio closer to the target by
    // increasing both the numerator and denominator.  The convergents of a mediant
    // search intersect the convergents of a continued fraction search, but with more
    // iterations and intermediate convergents.  However, it is probably simpler to implement
    // for compile-time than the Euclidean algorithm used for continued fractions.
    // This is not a binary search, although it looks similar.  A binary search would
    // not work well for this.
    template <typename low, typename high, bool done=false>
    struct srch {
        static_assert(std::ratio_greater_equal<high, low>::value, "Sqrt::srch got low > high");
        typedef typename mediant<low, high>::type med;
        typedef typename std::ratio_multiply<med, med> med_sq;
        static const bool med_over = std::ratio_greater<med_sq, r>::value;
        typedef typename select_t<low, med,  med_over>::type next_low;
        typedef typename select_t<med, high, med_over>::type next_high;
        static const bool stop = med::num > root_limit || med::den > root_limit;
        typedef typename srch<next_low, next_high, stop>::type type;
    };
    // The terminating specialization for the mediant search, triggered by the
    // last mediant's numerator or denominator exceeding root_limit
    // The low and high convergents are squared and compared with the input
    // in order to select the best result.
    template <typename low, typename high>
    struct srch<low, high, true> {
        typedef typename std::ratio_multiply<low,  low>  low_sq;
        typedef typename std::ratio_multiply<high, high> high_sq;
        typedef std::ratio_subtract<r, low_sq>  err_sq_low;
        typedef std::ratio_subtract<high_sq, r> err_sq_high;
        static const bool low_closer = std::ratio_less<err_sq_low, err_sq_high>::value;
        typedef typename select_t<low, high, low_closer>::type type;
    };
    // init checks the first approximation if it the root of a perfect square or if the numerator is zero
    // (in case on a non-canonical zero-valued ratio (unlikely, though))
    // when the input (r) is not a perfect square, this sets up the high and low ratios for the mediant
    // search.  These are set up in a way which ensures coprime numerators and denominators.
    // Then init kicks off the search loop.
    template <intmax_t num, intmax_t den, bool done = (num*num==r::num && den*den==r::den) || num==0>
    struct init {
        typedef std::ratio<num,     den + 1> low;
        typedef std::ratio<num + 1, den>     high;
        typedef typename srch<low, high>::type type;
    };
    // the init perfect square shortcut specialization
    template <intmax_t num, intmax_t den> struct init<num, den, true> {typedef std::ratio<num, den> type; };

    // An approximate numerator and denominator of the result is calculated 
    // by taking the integer square roots of the input ratio's num and den.
    // Then these values are passed to init for processing.
    static const intmax_t approx_num = iSqrt<r::num>::value;
    static const intmax_t approx_den = iSqrt<r::den>::value;
    typedef typename init<approx_num, approx_den>::type type;
};

// Diagnostic junk follows....
template <typename r>
void dump_ratio() {
    std::cout << r::num << '/' << r::den;
}

template <intmax_t N, intmax_t D>
std::ostream& operator << (std::ostream &os, typename std::ratio<N,D> r) {
    os << r.num << '/' << r.den;
    return os;
}

template <typename r>
double get_double() {
    return double(r::num) / r::den;
}

template <typename ...Rs> struct test_dump;

template <typename R, typename ...Rs>
struct test_dump<R, Rs...> {
    static void go() {
        std::cout << std::setprecision(14);
        double flt_in = get_double<R>();
        std::cout << " in:   " << R{} << " : " << flt_in << " (rat => double)\n";
        using rat_sqrt = typename Sqrt<R>::type;
        std::cout << " sqrt: " << rat_sqrt{} << '\n';
        double flt_rat_sqrt = get_double<rat_sqrt>();
        std::cout << flt_rat_sqrt << " (rat => double)\n";
        double sqrt_flt_in = sqrt(flt_in);
        double err = flt_rat_sqrt - sqrt_flt_in;
        std::cout << sqrt_flt_in << " : (sqrt in) ";
        std::cout << std::setprecision(4);
        std::cout << "err: " << err << "\n\n";
        test_dump<Rs...>::go();
    }
};

template <>
struct test_dump<> { static void go() {} };


int main() {
    test_dump<
        std::ratio<1060,   83>,
        std::ratio<1,      12494234>,
        std::ratio<82378,  1>,
        std::ratio<478723, 23423554>,
        std::ratio<2,      1>,
        std::ratio<2,      72>,
        std::ratio<2,      10001>,
        std::ratio<1000,   1>,
        std::ratio<10001,  2>
    >::go();
}

This is what the testing code at the bottom outputs:

 in:   1060/83 : 12.771084337349 (rat => double)
 sqrt: 28986/8111
3.5736653926766 (rat => double)
3.5736653924716 : (sqrt in) err: 2.05e-010

 in:   1/12494234 : 8.0036919430195e-008 (rat => double)
 sqrt: 174/615041
0.00028290796873704 (rat => double)
0.00028290796989515 : (sqrt in) err: -1.158e-012

 in:   82378/1 : 82378 (rat => double)
 sqrt: 585799/2041
287.01567858893 (rat => double)
287.01567901423 : (sqrt in) err: -4.253e-007

 in:   68389/3346222 : 0.020437675683203 (rat => double)
 sqrt: 88016/615667
0.14296039904689 (rat => double)
0.14296039900337 : (sqrt in) err: 4.352e-011

 in:   2/1 : 2 (rat => double)
 sqrt: 665857/470832
1.4142135623747 (rat => double)
1.4142135623731 : (sqrt in) err: 1.595e-012

 in:   1/36 : 0.027777777777778 (rat => double)
 sqrt: 1/6
0.16666666666667 (rat => double)
0.16666666666667 : (sqrt in) err: 0

 in:   2/10001 : 0.0001999800019998 (rat => double)
 sqrt: 9899/700000
0.014141428571429 (rat => double)
0.014141428569978 : (sqrt in) err: 1.45e-012

 in:   1000/1 : 1000 (rat => double)
 sqrt: 702247/22207
31.622776601972 (rat => double)
31.622776601684 : (sqrt in) err: 2.886e-010

 in:   10001/2 : 5000.5 (rat => double)
 sqrt: 700000/9899
70.714213556925 (rat => double)
70.714213564177 : (sqrt in) err: -7.252e-009

Answer 2

The problem is you use Newton's method to compute square root, which is correct if you want to get a numerical approximation, however, if you want to find best rational approximation, you gotta use continued fraction. Here is the result of my program:

hidden $ g++ -std=c++11 sqrt.cpp && ./a.out
sqrt(2/1) ~ 239/169, error=1.23789e-05, eps=0.0001
sqrt(2/1) ~ 114243/80782, error=5.41782e-11, eps=1e-10
sqrt(2/1) ~ 3880899/2744210, error=4.68514e-14, eps=1e-13
sqrt(2/1) ~ 131836323/93222358, error=0, eps=1e-16
sqrt(2/10001) ~ 1/71, error=5.69215e-05, eps=0.0001
sqrt(2/10001) ~ 1977/139802, error=2.18873e-11, eps=1e-10
sqrt(2/10001) ~ 13860/980099, error=7.36043e-15, eps=1e-13
sqrt(2/10001) ~ 1950299/137913860, error=3.64292e-17, eps=1e-16
sqrt(10001/2) ~ 495/7, error=7.21501e-05, eps=0.0001
sqrt(10001/2) ~ 980099/13860, error=3.68061e-11, eps=1e-10
sqrt(10001/2) ~ 415701778/5878617, error=1.42109e-14, eps=1e-13
sqrt(10001/2) ~ 970297515/13721393, error=0, eps=1e-16
sqrt(1060/83) ~ 461/129, error=2.19816e-05, eps=0.0001
sqrt(1060/83) ~ 2139943/598809, error=9.07718e-13, eps=1e-10
sqrt(1060/83) ~ 6448815/1804538, error=1.77636e-14, eps=1e-13
sqrt(1060/83) ~ 545951360/152770699, error=4.44089e-16, eps=1e-16
sqrt(1/12494234) ~ 1/3534, error=5.75083e-08, eps=0.0001
sqrt(1/12494234) ~ 32/113111, error=2.9907e-11, eps=1e-10
sqrt(1/12494234) ~ 419/1481047, error=6.02961e-14, eps=1e-13
sqrt(1/12494234) ~ 129879/459085688, error=4.49944e-18, eps=1e-16
sqrt(82378/1) ~ 18369/64, error=5.40142e-05, eps=0.0001
sqrt(82378/1) ~ 37361979/130174, error=1.16529e-11, eps=1e-10
sqrt(82378/1) ~ 1710431766/5959367, error=5.68434e-14, eps=1e-13
sqrt(82378/1) ~ 15563177213/54224136, error=0, eps=1e-16
sqrt(68389/3346222) ~ 197/1378, error=4.13769e-07, eps=0.0001
sqrt(68389/3346222) ~ 17801/124517, error=2.17069e-11, eps=1e-10
sqrt(68389/3346222) ~ 581697/4068938, error=4.30211e-15, eps=1e-13
sqrt(68389/3346222) ~ 16237871/113583000, error=2.77556e-17, eps=1e-16
sqrt(2/72) ~ 1/6, error=0, eps=0.0001
sqrt(10000/1) ~ 100/1, error=0, eps=0.0001
sqrt(0/20) ~ 0/1, error=0, eps=0.0001

My program finds the approximation with (almost) smallest numerator and denominator that satisfies the error bound. If you change the int to longer integer type in my code, it can go with much smaller eps.

My code (it compiles under g++-4.8.4, IMHO, the code would have been much simpler if constexpr is allowed):

#include <cstdint>
#include <cmath>
#include <iostream>
#include <ratio>
#include <type_traits>
#include <utility>

using namespace std;

using Zero = ratio<0>;
using One = ratio<1>;
template <typename R> using Square = ratio_multiply<R, R>;

// Find the largest integer N such that Predicate<N>::value is true.
template <template <intmax_t N> class Predicate, typename Enabled = void>
struct BinarySearch {
  template <intmax_t N>
  struct SafeDouble_ {
    const intmax_t static value = 2 * N;
    static_assert(value > 0, "Overflows when computing 2 * N");
  };

  template <intmax_t Lower, intmax_t Upper, typename Enabled1 = void>
  struct DoubleSidedSearch_ : DoubleSidedSearch_<Lower, Lower+(Upper-Lower)/2> {};

  template <intmax_t Lower, intmax_t Upper>
  struct DoubleSidedSearch_<Lower, Upper, typename enable_if<Upper-Lower==1>::type> : integral_constant<intmax_t, Lower> {};

  template <intmax_t Lower, intmax_t Upper>
  struct DoubleSidedSearch_<Lower, Upper, typename enable_if<(Upper-Lower>1 && Predicate<Lower+(Upper-Lower)/2>::value)>::type>
      : DoubleSidedSearch_<Lower+(Upper-Lower)/2, Upper> {};

  template <intmax_t Lower, typename Enabled1 = void>
  struct SingleSidedSearch_ : DoubleSidedSearch_<Lower, SafeDouble_<Lower>::value> {};

  template <intmax_t Lower>
  struct SingleSidedSearch_<Lower, typename enable_if<Predicate<SafeDouble_<Lower>::value>::value>::type>
      : SingleSidedSearch_<SafeDouble_<Lower>::value> {};

  const static intmax_t value = SingleSidedSearch_<1>::value;
};

template <template <intmax_t N> class Predicate>
struct BinarySearch<Predicate, typename enable_if<!Predicate<1>::value>::type> : integral_constant<intmax_t, 0> {};

// Find largest integer N such that N<=sqrt(R)
template <typename R>
struct Integer {
  template <intmax_t N> using Predicate_ = ratio_less_equal<ratio<N>, ratio_divide<R, ratio<N>>>;
  const static intmax_t value = BinarySearch<Predicate_>::value;
};

template <typename R>
struct IsPerfectSquare {
  const static intmax_t DenSqrt_ = Integer<ratio<R::den>>::value;
  const static intmax_t NumSqrt_ = Integer<ratio<R::num>>::value;
  const static bool value = DenSqrt_ * DenSqrt_ == R::den && NumSqrt_ * NumSqrt_ == R::num;
  using Sqrt = ratio<NumSqrt_, DenSqrt_>;
};

// Represents sqrt(P)-Q.
template <typename Tp, typename Tq>
struct Remainder {
  using P = Tp;
  using Q = Tq;
};

// Represents 1/R = I + Rem where R is a Remainder.
template <typename R>
struct Reciprocal {
  using P_ = typename R::P;
  using Q_ = typename R::Q;
  using Den_ = ratio_subtract<P_, Square<Q_>>;
  using A_ = ratio_divide<Q_, Den_>;
  using B_ = ratio_divide<P_, Square<Den_>>;
  const static intmax_t I_ = (A_::num + Integer<ratio_multiply<B_, Square<ratio<A_::den>>>>::value) / A_::den;
  using I = ratio<I_>;
  using Rem = Remainder<B_, ratio_subtract<I, A_>>;
};

// Expands sqrt(R) to continued fraction:
// f(x)=C1+1/(C2+1/(C3+1/(...+1/(Cn+x)))) = (U*x+V)/(W*x+1) and sqrt(R)=f(Rem).
// The error |f(Rem)-V| = |(U-W*V)x/(W*x+1)| <= |U-W*V|*Rem <= |U-W*V|/I' where
// I' is the integer part of reciprocal of Rem.
template <typename R, intmax_t N>
struct ContinuedFraction {
  template <typename T>
  using Abs_ = typename conditional<ratio_less<T, Zero>::value, ratio_subtract<Zero, T>, T>::type;

  using Last_ = ContinuedFraction<R, N-1>;
  using Reciprocal_ = Reciprocal<typename Last_::Rem>;
  using Rem = typename Reciprocal_::Rem;
  using I_ = typename Reciprocal_::I;
  using Den_ = ratio_add<typename Last_::W, I_>;
  using U = ratio_divide<typename Last_::V, Den_>;
  using V = ratio_divide<ratio_add<typename Last_::U, ratio_multiply<typename Last_::V, I_>>, Den_>;
  using W = ratio_divide<One, Den_>;
  using Error = Abs_<ratio_divide<ratio_subtract<U, ratio_multiply<V, W>>, typename Reciprocal<Rem>::I>>;
};

template <typename R>
struct ContinuedFraction<R, 1> {
  using U = One;
  using V = ratio<Integer<R>::value>;
  using W = Zero;
  using Rem = Remainder<R, V>;
  using Error = ratio_divide<One, typename Reciprocal<Rem>::I>;
};

template <typename R, typename Eps, intmax_t N=1, typename Enabled = void>
struct Sqrt_ : Sqrt_<R, Eps, N+1> {};

template <typename R, typename Eps, intmax_t N>
struct Sqrt_<R, Eps, N, typename enable_if<ratio_less_equal<typename ContinuedFraction<R, N>::Error, Eps>::value>::type> {
  using type = typename ContinuedFraction<R, N>::V;
};

template <typename R, typename Eps, typename Enabled = void>
struct Sqrt {
  static_assert(ratio_greater_equal<R, Zero>::value, "R can't be negative");
};

template <typename R, typename Eps>
struct Sqrt<R, Eps, typename enable_if<ratio_greater_equal<R, Zero>::value && IsPerfectSquare<R>::value>::type> {
  using type = typename IsPerfectSquare<R>::Sqrt;
};

template <typename R, typename Eps>
struct Sqrt<R, Eps, typename enable_if<(ratio_greater_equal<R, Zero>::value && !IsPerfectSquare<R>::value)>::type> : Sqrt_<R, Eps> {};

// Test finding sqrt(N/D) with error 1/Eps
template <intmax_t N, intmax_t D, intmax_t Eps>
void test() {
  using T = typename Sqrt<ratio<N, D>, ratio<1, Eps>>::type;
  cout << "sqrt(" << N << "/" << D << ") ~ " << T::num << "/" << T::den << ", "
       << "error=" << abs(sqrt(N/(double)D) - T::num/(double)T::den) << ", "
       << "eps=" << 1/(double)Eps << endl;
}

template <intmax_t N, intmax_t D>
void testAll() {
  test<N, D, 10000>();
  test<N, D, 10000000000>();
  test<N, D, 10000000000000>();
  test<N, D, 10000000000000000>();
}

int main() {
  testAll<2, 1>();
  testAll<2, 10001>();
  testAll<10001, 2>();

  testAll<1060, 83>();
  testAll<1, 12494234>();
  testAll<82378, 1>();
  testAll<68389, 3346222>();

  test<2, 72, 10000>();
  test<10000, 1, 10000>();
  test<0, 20, 10000>();
  // static assertion failure.
  // test<-10001, 2, 100000>();
}

---

Modified BinarySearch for MSVC 2013

Due to bugs in the template deduction implementation of the Visual Studio 2013 compiler, the BinarySearch has to be modified when building on that platform:

template <template <std::intmax_t N> class Predicate, typename enabled = void>
struct BinarySearch {
    template <std::intmax_t N>
    struct SafeDouble_ {
        static const std::intmax_t value = 2 * N;
        static_assert(value > 0, "Overflows when computing 2 * N");
    };

    template <intmax_t Lower, intmax_t Upper, typename Condition1 = void, typename Condition2 = void>
    struct DoubleSidedSearch_ : DoubleSidedSearch_<Lower, Upper,
        typename std::conditional<(Upper - Lower == 1), std::true_type, std::false_type>::type,
        typename std::conditional<((Upper - Lower>1 && Predicate<Lower + (Upper - Lower) / 2>::value)), std::true_type, std::false_type>::type> {};

    template <intmax_t Lower, intmax_t Upper>
    struct DoubleSidedSearch_<Lower, Upper, std::false_type, std::false_type> : DoubleSidedSearch_<Lower, Lower + (Upper - Lower) / 2> {};

    template <intmax_t Lower, intmax_t Upper, typename Condition2>
    struct DoubleSidedSearch_<Lower, Upper, std::true_type, Condition2> : std::integral_constant<intmax_t, Lower>{};

    template <intmax_t Lower, intmax_t Upper, typename Condition1>
    struct DoubleSidedSearch_<Lower, Upper, Condition1, std::true_type> : DoubleSidedSearch_<Lower + (Upper - Lower) / 2, Upper>{};

    template <std::intmax_t Lower, class enabled1 = void>
    struct SingleSidedSearch_ : SingleSidedSearch_<Lower, typename std::conditional<Predicate<SafeDouble_<Lower>::value>::value, std::true_type, std::false_type>::type>{};

    template <std::intmax_t Lower>
    struct SingleSidedSearch_<Lower, std::false_type> : DoubleSidedSearch_<Lower, SafeDouble_<Lower>::value> {};

    template <std::intmax_t Lower>
    struct SingleSidedSearch_<Lower, std::true_type> : SingleSidedSearch_<SafeDouble_<Lower>::value>{};

    const static std::intmax_t value = SingleSidedSearch_<1>::value;
};

Answer 3

I think I understand the problem, what you are saying is that due to (multiplication) overflow in the integers, you can not do more than 4 or 5 iterations before the overflow... and that can still be a bad approximation.

I think the only thing you can do about it is to have better initial estimates. And probably have that hard coded. In the example I show the range is large enough, but you can hard code more initial guesses to enlarge the range even more.

#include<ratio>
#include<iostream>

namespace detail
{

    template<class R>
    struct estimatimate{
        typedef 
        typename std::conditional<
            std::ratio_less_equal<R, std::ratio<1,100>>::value, 
            std::ratio<1,10>,
            typename std::conditional<
                std::ratio_less_equal<R, std::ratio<25,100>>::value, 
                std::ratio<5,10>,
                typename std::conditional<
                    std::ratio_less_equal<R, std::ratio<1>>::value, 
                    std::ratio<1>,
                    typename std::conditional<
                        std::ratio_less_equal<R, std::ratio<25>>::value, 
                        std::ratio<5>,
                        typename std::conditional<
                            std::ratio_less_equal<R, std::ratio<100>>::value,
                            std::ratio<10>,
                            typename std::conditional<
                                std::ratio_less_equal<R, std::ratio<2500>>::value,
                                std::ratio<50>,
                                typename std::conditional<
                                    std::ratio_less_equal<R, std::ratio<10000>>::value,
                                    std::ratio<100>,
                                    std::ratio<500>
                                >::type
                            >::type
                        >::type
                    >::type
                >::type
            >::type
        >::type type;
    };

    template<class N , class K = typename estimatimate<N>::type, std::intmax_t RecursionDepth = 4>
    struct ratio_sqrt_impl
    {
    //    static_assert(/* N is std::ratio */);

        // 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;

template<class Ratio>
double value(){
    return ((double)Ratio::num)/Ratio::den;
}

#include<cmath>
using namespace std;
int main(){
//  std::cout << value<std::ratio<1,10>>() << std::endl;
//  std::cout << value<ratio_sqrt<std::ratio<4,1>>>() << std::endl;

    std::cout << value<ratio_sqrt<std::ratio<1,128>>>() << " " << sqrt(value<std::ratio<1,128>>()) << std::endl;
    std::cout << value<ratio_sqrt<std::ratio<5000,1>>>() << " " << sqrt(value<std::ratio<5000,1>>()) << std::endl;

}
//output:
// 0.0883883 0.0883883
// 70.7107 70.7107

Unfortunately given the overflow, if you want to have less iterations you have to make the initial guess even better.