c – 编译时std :: ratio的平方根的有理逼近
我试图在编译时找到std :: ratio的平方根的有理逼近.从定义的坐标转换参数中导出
ellipsoid parameters非常有用,坐标转换本身定义为std :: ratio.
有一个question about finding powers/roots of 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; 产生结果:
对某些应用程序来说肯定是不错的. 问题 >这里最大的问题是固定的递归深度.这些比率非常快,因此对于根>> 100,这疯狂溢出.但是,递归深度太小,而且你失去了所有的准确性. 是否有一种好的方法可以将递归调整到溢出深度限制,然后将类型设置为之前的一两次迭代? (我说了一些迭代,因为将余量保持在整数大小以便以后进行进一步的计算可能会很好) > 4的初始条件在产生根的最低误差方面似乎是非常神奇的. 100,但有没有更有条理的方法来设置它? 编辑: 我不是在寻找constexpr的任何解决方案,因为我必须支持的编译器不能统一使用它. 增加递归深度的问题是std :: ratio的num / denom仅在几次递归后溢出.所表示的平方根的准确性实际上是正常的,但我需要找到一个通用的解决方案,将递归深度限制在比率不会溢出的点(因此,不要编译).例如. ratio_sqrt<的std ::比< 2 – ;>可以在溢出之前进入深度5,但是ratio_sqrt< std :: ratio< 1000>>限于4. 解决方法
问题是你使用牛顿的方法计算平方根,如果你想得到一个数值近似是正确的,但是,如果你想找到最合理的近似,你必须使用
continued fraction.这是我的程序的结果:
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,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,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,eps=0.0001 sqrt(10000/1) ~ 100/1,eps=0.0001 sqrt(0/20) ~ 0/1,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,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,intmax_t N> struct Sqrt_<R,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 Enabled = void> struct Sqrt { static_assert(ratio_greater_equal<R,"R can't be negative"); }; template <typename R,typename Eps> struct Sqrt<R,typename enable_if<ratio_greater_equal<R,Zero>::value && IsPerfectSquare<R>::value>::type> { using type = typename IsPerfectSquare<R>::Sqrt; }; template <typename R,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,10000000000>(); test<N,10000000000000>(); test<N,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>(); } 为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,typename Condition1 = void,typename Condition2 = void> struct DoubleSidedSearch_ : DoubleSidedSearch_<Lower,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::false_type>::type> {}; template <intmax_t Lower,intmax_t Upper> struct DoubleSidedSearch_<Lower,std::false_type,std::false_type> : DoubleSidedSearch_<Lower,Lower + (Upper - Lower) / 2> {}; template <intmax_t Lower,typename Condition2> struct DoubleSidedSearch_<Lower,Condition2> : std::integral_constant<intmax_t,Lower>{}; template <intmax_t Lower,typename Condition1> struct DoubleSidedSearch_<Lower,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::false_type>::type>{}; template <std::intmax_t Lower> struct SingleSidedSearch_<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; }; (编辑:李大同) 【声明】本站内容均来自网络,其相关言论仅代表作者个人观点,不代表本站立场。若无意侵犯到您的权利,请及时与联系站长删除相关内容! |