【问题标题】:Seeding the Newton iteration for cube root efficiently有效地为立方根播种牛顿迭代
【发布时间】:2011-09-18 18:30:31
【问题描述】:

如何有效地找到数字的立方根? 我认为可以使用 Newton-Raphson 方法,但我不知道如何以编程方式猜测初始解决方案以最小化迭代次数。

【问题讨论】:

  • “数字”是什么意思?

标签: algorithm math


【解决方案1】:

这是一个看似复杂的问题。 Here 对一些可能的方法进行了很好的调查。

【讨论】:

  • 我可能听起来很粗鲁,但我认为如果链接失效,您应该在答案中添加 一些 内容。 (如果需要,复制和粘贴)
【解决方案2】:

鉴于超过 Accepted Answer 的“链接失效”,我将给出一个更加独立的答案,重点关注快速获得适合超线性迭代的初始猜测的主题。

metomerist (Wayback link) 的“调查”为各种起始值/迭代组合(包括牛顿法和哈雷法)提供了一些时序比较。它引用了W. Kahan,“计算一个真正的立方体根”和K. Turkowski,“计算立方体根”的作品。

metamarist 用这个 sn-p 更新了 W. Kahan 的 DEC-VAX 时代比特摆弄技术,它“假定 32 位整数”并依赖 IEEE 754 格式的双精度“生成具有 5 位精度的初始估计":

inline double cbrt_5d(double d) 
{ 
   const unsigned int B1 = 715094163; 
   double t = 0.0; 
   unsigned int* pt = (unsigned int*) &t; 
   unsigned int* px = (unsigned int*) &d; 
   pt[1]=px[1]/3+B1; 
   return t; 
} 

K. Turkowski 的代码通过在 float fr 上的传统二次方缩放提供了稍微更高的精度(“大约 6 位”),然后在区间 [0.125,1.0) 上对其立方根进行二次逼近:

/* Compute seed with a quadratic qpproximation */
fr = (-0.46946116F * fr + 1.072302F) * fr + 0.3812513F;/* 0.5<=fr<1 */

以及随后恢复二的指数(调整为三分之一)。指数/尾数提取和恢复利用math library callsfrexpldexp

与其他立方根“种子”近似值的比较

要了解这些立方根近似值,我们需要将它们与其他可能的形式进行比较。首先判断的标准:我们考虑区间[1/8,1]的近似值,我们使用best(最小化最大)相对误差。

也就是说,如果f(x)x^{1/3} 的近似值,我们会发现它的相对误差:

        error_rel = max | f(x)/x^(1/3) - 1 | on [1/8,1]

最简单的近似当然是在区间上使用单个常数,在这种情况下,最好的相对误差是通过选择f_0(x) = sqrt(2)/2(端点值的几何平均值)来实现的。这给出了 1.27 位的相对精度,这是牛顿迭代的快速但肮脏的起点。

更好的近似是最好的一次多项式:

 f_1(x) = 0.6042181313*x + 0.4531635984

这给出了 4.12 位的相对精度,这是一个很大的改进,但比 Kahan 和 Turkowski 各自的方法所承诺的 5-6 位的相对精度还差。但它只是在球场上,只使用一次乘法(和一次加法)。

最后,如果我们允许自己进行除法而不是乘法呢?事实证明,通过一次除法和两次“加法”,我们可以得到最好的线性分数函数:

 f_M(x) = 1.4774329094 - 0.8414323527/(x+0.7387320679) 

提供 7.265 位的相对准确度。

乍一看,这似乎是一种有吸引力的方法,但一个古老的经验法则是将 FP 除法的成本视为三个 FP 乘法(并且几乎忽略了加法和减法)。然而,对于当前的 FPU 设计,这是不现实的。虽然乘法到加法/减法的相对成本已经下降,但在大多数情况下是两倍甚至相等,但除法的成本并没有下降,而是经常上升到乘法成本的 7-10 倍。因此,我们必须吝啬我们的部门运作。

【讨论】:

    【解决方案3】:
    static double cubeRoot(double num) {
        double x = num;
    
        if(num >= 0) {
            for(int i = 0; i < 10 ; i++) {
                x = ((2 * x * x * x) + num ) / (3 * x * x); 
            }
        } 
        return x;
    }
    

    【讨论】:

      【解决方案4】:

      似乎优化问题已经得到解决,但我想对此处发布的 cubeRoot() 函数进行改进,以便其他人在此页面上绊倒寻找快速立方根算法。

      现有算法运行良好,但在 0-100 范围之外它会给出错误的结果。

      这是一个修订版,适用于 -/+1 万亿 (1E15) 之间的数字。如果您需要处理更大的数字,只需使用更多的迭代。

      static double cubeRoot( double num ){
          boolean neg = ( num < 0 );
          double x = Math.abs( num );
          for( int i = 0, iterations = 60; i < iterations; i++ ){
              x = ( ( 2 * x * x * x ) + num ) / ( 3 * x * x ); 
          }
          if( neg ){ return 0 - x; }
          return x;
      }
      

      关于优化,我猜原始发帖人是在询问如何在给定任意输入大小的情况下预测准确结果的最小迭代次数。但对于大多数一般情况,优化带来的收益似乎不值得增加复杂性。即使使用上述功能,100 次迭代在普通消费类硬件上花费的时间也不到 0.2 毫秒。如果速度是最重要的,我会考虑使用预先计算的查找表。但这来自桌面开发人员,而不是嵌入式系统工程师。

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 1970-01-01
        • 2022-01-16
        • 1970-01-01
        • 2015-04-26
        • 2019-06-27
        • 1970-01-01
        • 1970-01-01
        • 2014-03-07
        相关资源
        最近更新 更多