编译时std :: ratio的平方根的有理逼近

尼古拉斯·霍尔特豪斯(Nicolas Holthaus)

我试图std::ratio在编译时找到a的平方根的有理近似值从定义的参数中导出椭球参数以进行坐标转换将非常有用,坐标本身定义为std::ratio

关于找到std :: ratio的幂/根有一个问题,但是作为该问题的条件,如果比率没有整数根,则可以失败,这与我想要的相反。相反,我想找到我能找到的最接近的合理近似值。

我已经提出了以下元程序,该程序基于牛顿-拉夫森法Newton-Raphson Method)来计算根,已知该方法只需几次迭代即可产生(相对)准确的结果:

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;

使用示例用法:

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

产生结果:

1.46538e-05 // sqrt(2)

4.64611e-08 // sqrt(4)

2.38737e-15 // sqrt(10)

对于某些应用程序来说肯定是不错的选择。

问题所在

  1. 这里最大的问题是固定的递归深度。这些比率很快就变得很大,因此对于大于100的根,这会疯狂地溢出。但是,递归深度太小会损失所有精度。

有没有一种好的方法可以使递归适应溢出深度限制,然后在此之前将类型设置为一个或两个迭代?(我说了几次迭代,因为将净空保持在整数大小可能会更好,以便以后进行进一步的计算)

  1. 就产生根数小于100的最低误差而言,初始条件4似乎是非常神奇的,但是有没有更系统的方法来设置?

编辑:

我没有使用寻找任何解决方案constexpr,因为我必须支持的编译器并没有统一的解决方案

增加递归深度的问题是std::ratio仅几次递归后溢出的数字/单位所表示的平方根的精度实际上是可以的,但是我需要找到一个通用的解决方案,将递归深度限制到比率不溢出(因此不编译)的程度。例如ratio_sqrt<std::ratio<2>>,在溢出之前可以到达深度5,但ratio_sqrt<std::ratio<1000>>限于4。

Kan Li

问题是您使用牛顿法计算平方根,如果要获得数值逼近,这是正确的,但是,如果要找到最佳有理逼近,则必须使用连续分数这是我的程序的结果:

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

我的程序找到了(几乎)最小的分子和分母满足误差范围的近似值。如果int在我的代码中将更改为更长的整数类型,则可以使用更小的整数类型eps

我的代码(它在g ++-4.8.4,恕我直言下编译,如果允许使用constexpr,则代码会简单得多):

#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>();
}

BinarySearch为MSVC 2013修改

由于Visual Studio 2013编译器的模板推导实现中存在错误,因此BinarySearch在该平台上构建时必须进行修改:

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

本文收集自互联网,转载请注明来源。

如有侵权,请联系 [email protected] 删除。

编辑于
0

我来说两句

0 条评论
登录 后参与评论

相关文章

如何计算Coq中自然数或有理数的平方根?

快速整数平方根逼近

为什么在编译时已知的常数的平方根不需要链接数学库?

平方根正数时math.sqrt域错误

std :: norm(std :: complex)使用平方根而不是快速实现

禁用Matlab符号有理逼近

使用Nasm x86 Assembly计算平方根时输出不正确

牛顿-拉夫森法计算平方根时的While循环条件

当要求更精确的答案时,平方根计算器会中断

Java:使用Unicode上划线显示平方根时的字符串长度?

理性的平方根

角的平方根

平方根计算

用牛顿法将R循环逼近正实数的平方根

使用现代C ++时,带有私有构造函数的类的std :: vector无法编译

MSVC在最大及附近产生std :: uint64_t的不正确平方根。错误或预期?

平方系列的平方根

编译器通过错误:尝试使用std :: pair向量时没有匹配函数

为什么此std :: string C ++代码没有给出编译时错误?

在C ++中,如果必须在编译时确定大小,则std :: array有什么意义?

为什么在转换为std :: function时没有编译错误

确保std :: map在编译时具有具体大小的优雅方法

在另一个数组中找到n的平方根时,布尔结果不正确

在 Windows 10 命令提示符下运行 PHP 脚本时如何打印出平方根符号?

我在将平方根函数应用于由 float (C#) 组成的方程时遇到问题

Java中的平方根和平方根

是否有BigDecimal上平方根的库?

没有得到正确的平方根

PHP平方根计算器(带有HTML)

TOP 榜单

热门标签

归档