使用const non优化pow()

我在我的代码中有热点,我正在做pow()占用了我执行时间的10-20%。

我对pow(x,y)输入是非常具体的,所以我想知道是否有办法以更高的性能滚动两个pow()近似值(每个指数一个):

  • 我有两个常量指数:2.4和1 / 2.4。
  • 当指数为2.4时,x将在范围内(0.090473935,1.0)。
  • 当指数为1 / 2.4时,x将在范围内(0.0031308,1.0)。
  • 我正在使用SSE / AVX float向量。 如果平台的细节可以被利用,就在!
  • 尽管我对全精度( float )算法也很感兴趣,但最大误差率在0.01%左右是理想的。

    我已经使用了一个快速的pow()近似,但它并没有考虑到这些限制。 有没有可能做得更好?


    在IEEE 754黑客手中,这是另一种解决方案,速度更快,不那么“神奇”。 它在大约十几个时钟周期内实现了0.08%的误差范围(对于p = 2.4,在Intel Merom CPU上)。

    浮点数最初是作为对数的近似值发明的,因此您可以使用整数值作为log2的近似值。 通过将convert-from-integer指令应用于浮点值,可以实现这一点,以获得另一个浮点值。

    要完成pow计算,可以乘以一个常数因子,并将对数转换为convert-to-integer指令。 在SSE上,相关说明是cvtdq2pscvtps2dq

    虽然这并不那么简单。 IEEE 754中的指数字段有符号,偏差值为127,表示指数为零。 在您乘以对数之前必须删除此偏差,并在您取幂前重新添加。 此外,通过减法进行的偏差调整不会在零点上起作用。 幸运的是,两种调整都可以通过预先乘以一个常数来实现。

    x^p
    = exp2( p * log2( x ) )
    = exp2( p * ( log2( x ) + 127 - 127 ) - 127 + 127 )
    = cvtps2dq( p * ( log2( x ) + 127 - 127 - 127 / p ) )
    = cvtps2dq( p * ( log2( x ) + 127 - log2( exp2( 127 - 127 / p ) ) )
    = cvtps2dq( p * ( log2( x * exp2( 127 / p - 127 ) ) + 127 ) )
    = cvtps2dq( p * ( cvtdq2ps( x * exp2( 127 / p - 127 ) ) ) )
    

    exp2( 127 / p - 127 )是常数。 这个函数相当专业:它不适用于小分数指数,因为常数因数随着指数的倒数呈指数增长并且会溢出。 它不适用于负指数。 大指数导致高错误,因为尾数位与乘法指数位混杂在一起。

    但是,这只是4条快速指令。 预乘,从“整数”(以对数)转换为乘方,转换为“整数”(以对数为单位)。 这种SSE实施的转换速度非常快。 我们还可以将一个额外的常系数挤入第一乘法中。

    template< unsigned expnum, unsigned expden, unsigned coeffnum, unsigned coeffden >
    __m128 fastpow( __m128 arg ) {
            __m128 ret = arg;
    //      std::printf( "arg = %,vgn", ret );
            // Apply a constant pre-correction factor.
            ret = _mm_mul_ps( ret, _mm_set1_ps( exp2( 127. * expden / expnum - 127. )
                    * pow( 1. * coeffnum / coeffden, 1. * expden / expnum ) ) );
    //      std::printf( "scaled = %,vgn", ret );
            // Reinterpret arg as integer to obtain logarithm.
            asm ( "cvtdq2ps %1, %0" : "=x" (ret) : "x" (ret) );
    //      std::printf( "log = %,vgn", ret );
            // Multiply logarithm by power.
            ret = _mm_mul_ps( ret, _mm_set1_ps( 1. * expnum / expden ) );
    //      std::printf( "powered = %,vgn", ret );
            // Convert back to "integer" to exponentiate.
            asm ( "cvtps2dq %1, %0" : "=x" (ret) : "x" (ret) );
    //      std::printf( "result = %,vgn", ret );
            return ret;
    }
    

    一些指数= 2.4的试验表明,这一直高估了5%左右。 (例程总是保证高估。)你可以简单地乘以0.95,但是更多的指令会让我们精确到4位十进制数字,这对图形应该足够了。

    关键是要将过高估计与低估相匹配,并取平均值。

  • 计算x ^ 0.8:四条指令,错误〜+ 3%。
  • 计算x ^ -0.4:一个rsqrtps 。 (这非常准确,但是确实牺牲了零工作能力。)
  • 计算x ^ 0.4:一个mulps
  • 计算x ^ -0.2:一个rsqrtps
  • 计算x ^ 2:一个mulps
  • 计算x ^ 3:一个mulps
  • x ^ 2.4 = x ^ 2 * x ^ 0.4:一个mulps 。 这是高估。
  • x ^ 2.4 = x ^ 3 * x ^ -0.4 * x ^ -0.2:两个mulps 。 这是低估。
  • 上面的平均值:一个addps ,一个mulps

  • 指令计数:14,其中包括两次延迟= 5的转换和两次倒数平方根估计,吞吐量= 4。

    为了恰当地取平均值,我们希望通过预期误差对估计值进行加权。 低估将错误提高到0.6对0.4的功率,所以我们预计它是错误的1.5倍。 加权不会添加任何说明; 它可以在预设因素中完成。 调用系数a:a ^ 0.5 = 1.5 a ^ -0.75,a = 1.38316186。

    最终的误差约为0.015%,或比最初的fastpow结果好2个数量级。 对于具有volatile源变量和目标变量的繁忙循环,运行时间大约是十几个周期......虽然它与迭代重叠,但真实世界的用法也将看到指令级并行性。 考虑到SIMD,这是每3个周期一个标量结果的吞吐量!

    int main() {
            __m128 const x0 = _mm_set_ps( 0.01, 1, 5, 1234.567 );
            std::printf( "Input: %,vgn", x0 );
    
            // Approx 5% accuracy from one call. Always an overestimate.
            __m128 x1 = fastpow< 24, 10, 1, 1 >( x0 );
            std::printf( "Direct x^2.4: %,vgn", x1 );
    
            // Lower exponents provide lower initial error, but too low causes overflow.
            __m128 xf = fastpow< 8, 10, int( 1.38316186 * 1e9 ), int( 1e9 ) >( x0 );
            std::printf( "1.38 x^0.8: %,vgn", xf );
    
            // Imprecise 4-cycle sqrt is still far better than fastpow, good enough.
            __m128 xfm4 = _mm_rsqrt_ps( xf );
            __m128 xf4 = _mm_mul_ps( xf, xfm4 );
    
            // Precisely calculate x^2 and x^3
            __m128 x2 = _mm_mul_ps( x0, x0 );
            __m128 x3 = _mm_mul_ps( x2, x0 );
    
            // Overestimate of x^2 * x^0.4
            x2 = _mm_mul_ps( x2, xf4 );
    
            // Get x^-0.2 from x^0.4. Combine with x^-0.4 into x^-0.6 and x^2.4.
            __m128 xfm2 = _mm_rsqrt_ps( xf4 );
            x3 = _mm_mul_ps( x3, xfm4 );
            x3 = _mm_mul_ps( x3, xfm2 );
    
            std::printf( "x^2 * x^0.4: %,vgn", x2 );
            std::printf( "x^3 / x^0.6: %,vgn", x3 );
            x2 = _mm_mul_ps( _mm_add_ps( x2, x3 ), _mm_set1_ps( 1/ 1.960131704207789 ) );
            // Final accuracy about 0.015%, 200x better than x^0.8 calculation.
            std::printf( "average = %,vgn", x2 );
    }
    

    那么...对不起,我无法更快发布。 并将其扩展到x ^ 1 / 2.4作为练习; v)。


    更新统计信息

    我实施了一个小测试线束和两个x(5/12)的情况对应上述。

    #include <cstdio>
    #include <xmmintrin.h>
    #include <cmath>
    #include <cfloat>
    #include <algorithm>
    using namespace std;
    
    template< unsigned expnum, unsigned expden, unsigned coeffnum, unsigned coeffden >
    __m128 fastpow( __m128 arg ) {
        __m128 ret = arg;
    //  std::printf( "arg = %,vgn", ret );
        // Apply a constant pre-correction factor.
        ret = _mm_mul_ps( ret, _mm_set1_ps( exp2( 127. * expden / expnum - 127. )
            * pow( 1. * coeffnum / coeffden, 1. * expden / expnum ) ) );
    //  std::printf( "scaled = %,vgn", ret );
        // Reinterpret arg as integer to obtain logarithm.
        asm ( "cvtdq2ps %1, %0" : "=x" (ret) : "x" (ret) );
    //  std::printf( "log = %,vgn", ret );
        // Multiply logarithm by power.
        ret = _mm_mul_ps( ret, _mm_set1_ps( 1. * expnum / expden ) );
    //  std::printf( "powered = %,vgn", ret );
        // Convert back to "integer" to exponentiate.
        asm ( "cvtps2dq %1, %0" : "=x" (ret) : "x" (ret) );
    //  std::printf( "result = %,vgn", ret );
        return ret;
    }
    
    __m128 pow125_4( __m128 arg ) {
        // Lower exponents provide lower initial error, but too low causes overflow.
        __m128 xf = fastpow< 4, 5, int( 1.38316186 * 1e9 ), int( 1e9 ) >( arg );
    
        // Imprecise 4-cycle sqrt is still far better than fastpow, good enough.
        __m128 xfm4 = _mm_rsqrt_ps( xf );
        __m128 xf4 = _mm_mul_ps( xf, xfm4 );
    
        // Precisely calculate x^2 and x^3
        __m128 x2 = _mm_mul_ps( arg, arg );
        __m128 x3 = _mm_mul_ps( x2, arg );
    
        // Overestimate of x^2 * x^0.4
        x2 = _mm_mul_ps( x2, xf4 );
    
        // Get x^-0.2 from x^0.4, and square it for x^-0.4. Combine into x^-0.6.
        __m128 xfm2 = _mm_rsqrt_ps( xf4 );
        x3 = _mm_mul_ps( x3, xfm4 );
        x3 = _mm_mul_ps( x3, xfm2 );
    
        return _mm_mul_ps( _mm_add_ps( x2, x3 ), _mm_set1_ps( 1/ 1.960131704207789 * 0.9999 ) );
    }
    
    __m128 pow512_2( __m128 arg ) {
        // 5/12 is too small, so compute the sqrt of 10/12 instead.
        __m128 x = fastpow< 5, 6, int( 0.992245 * 1e9 ), int( 1e9 ) >( arg );
        return _mm_mul_ps( _mm_rsqrt_ps( x ), x );
    }
    
    __m128 pow512_4( __m128 arg ) {
        // 5/12 is too small, so compute the 4th root of 20/12 instead.
        // 20/12 = 5/3 = 1 + 2/3 = 2 - 1/3. 2/3 is a suitable argument for fastpow.
        // weighting coefficient: a^-1/2 = 2 a; a = 2^-2/3
        __m128 xf = fastpow< 2, 3, int( 0.629960524947437 * 1e9 ), int( 1e9 ) >( arg );
        __m128 xover = _mm_mul_ps( arg, xf );
    
        __m128 xfm1 = _mm_rsqrt_ps( xf );
        __m128 x2 = _mm_mul_ps( arg, arg );
        __m128 xunder = _mm_mul_ps( x2, xfm1 );
    
        // sqrt2 * over + 2 * sqrt2 * under
        __m128 xavg = _mm_mul_ps( _mm_set1_ps( 1/( 3 * 0.629960524947437 ) * 0.999852 ),
                                    _mm_add_ps( xover, xunder ) );
    
        xavg = _mm_mul_ps( xavg, _mm_rsqrt_ps( xavg ) );
        xavg = _mm_mul_ps( xavg, _mm_rsqrt_ps( xavg ) );
        return xavg;
    }
    
    __m128 mm_succ_ps( __m128 arg ) {
        return (__m128) _mm_add_epi32( (__m128i) arg, _mm_set1_epi32( 4 ) );
    }
    
    void test_pow( double p, __m128 (*f)( __m128 ) ) {
        __m128 arg;
    
        for ( arg = _mm_set1_ps( FLT_MIN / FLT_EPSILON );
                ! isfinite( _mm_cvtss_f32( f( arg ) ) );
                arg = mm_succ_ps( arg ) ) ;
    
        for ( ; _mm_cvtss_f32( f( arg ) ) == 0;
                arg = mm_succ_ps( arg ) ) ;
    
        std::printf( "Domain from %gn", _mm_cvtss_f32( arg ) );
    
        int n;
        int const bucket_size = 1 << 25;
        do {
            float max_error = 0;
            double total_error = 0, cum_error = 0;
            for ( n = 0; n != bucket_size; ++ n ) {
                float result = _mm_cvtss_f32( f( arg ) );
    
                if ( ! isfinite( result ) ) break;
    
                float actual = ::powf( _mm_cvtss_f32( arg ), p );
    
                float error = ( result - actual ) / actual;
                cum_error += error;
                error = std::abs( error );
                max_error = std::max( max_error, error );
                total_error += error;
    
                arg = mm_succ_ps( arg );
            }
    
            std::printf( "error max = %8gt" "avg = %8gt" "|avg| = %8gt" "to %8gn",
                        max_error, cum_error / n, total_error / n, _mm_cvtss_f32( arg ) );
        } while ( n == bucket_size );
    }
    
    int main() {
        std::printf( "4 insn x^12/5:n" );
        test_pow( 12./5, & fastpow< 12, 5, 1059, 1000 > );
        std::printf( "14 insn x^12/5:n" );
        test_pow( 12./5, & pow125_4 );
        std::printf( "6 insn x^5/12:n" );
        test_pow( 5./12, & pow512_2 );
        std::printf( "14 insn x^5/12:n" );
        test_pow( 5./12, & pow512_4 );
    }
    

    输出:

    4 insn x^12/5:
    Domain from 1.36909e-23
    error max =      inf    avg =      inf  |avg| =      inf    to 8.97249e-19
    error max =  2267.14    avg =  139.175  |avg| =  139.193    to 5.88021e-14
    error max = 0.123606    avg = -0.000102963  |avg| = 0.0371122   to 3.85365e-09
    error max = 0.123607    avg = -0.000108978  |avg| = 0.0368548   to 0.000252553
    error max =  0.12361    avg = 7.28909e-05   |avg| = 0.037507    to  16.5513
    error max = 0.123612    avg = -0.000258619  |avg| = 0.0365618   to 1.08471e+06
    error max = 0.123611    avg = 8.70966e-05   |avg| = 0.0374369   to 7.10874e+10
    error max =  0.12361    avg = -0.000103047  |avg| = 0.0371122   to 4.65878e+15
    error max = 0.123609    avg =      nan  |avg| =      nan    to 1.16469e+16
    14 insn x^12/5:
    Domain from 1.42795e-19
    error max =      inf    avg =      nan  |avg| =      nan    to 9.35823e-15
    error max = 0.000936462 avg = 2.0202e-05    |avg| = 0.000133764 to 6.13301e-10
    error max = 0.000792752 avg = 1.45717e-05   |avg| = 0.000129936 to 4.01933e-05
    error max = 0.000791785 avg = 7.0132e-06    |avg| = 0.000129923 to  2.63411
    error max = 0.000787589 avg = 1.20745e-05   |avg| = 0.000129347 to   172629
    error max = 0.000786553 avg = 1.62351e-05   |avg| = 0.000132397 to 1.13134e+10
    error max = 0.000785586 avg = 8.25205e-06   |avg| = 0.00013037  to 6.98147e+12
    6 insn x^5/12:
    Domain from 9.86076e-32
    error max = 0.0284339   avg = 0.000441158   |avg| = 0.00967327  to 6.46235e-27
    error max = 0.0284342   avg = -5.79938e-06  |avg| = 0.00897913  to 4.23516e-22
    error max = 0.0284341   avg = -0.000140706  |avg| = 0.00897084  to 2.77556e-17
    error max = 0.028434    avg = 0.000440504   |avg| = 0.00967325  to 1.81899e-12
    error max = 0.0284339   avg = -6.11153e-06  |avg| = 0.00897915  to 1.19209e-07
    error max = 0.0284298   avg = -0.000140597  |avg| = 0.00897084  to 0.0078125
    error max = 0.0284371   avg = 0.000439748   |avg| = 0.00967319  to      512
    error max = 0.028437    avg = -7.74294e-06  |avg| = 0.00897924  to 3.35544e+07
    error max = 0.0284369   avg = -0.000142036  |avg| = 0.00897089  to 2.19902e+12
    error max = 0.0284368   avg = 0.000439183   |avg| = 0.0096732   to 1.44115e+17
    error max = 0.0284367   avg = -7.41244e-06  |avg| = 0.00897923  to 9.44473e+21
    error max = 0.0284366   avg = -0.000141706  |avg| = 0.00897088  to 6.1897e+26
    error max = 0.485129    avg = -0.0401671    |avg| = 0.048422    to 4.05648e+31
    error max = 0.994932    avg = -0.891494 |avg| = 0.891494    to 2.65846e+36
    error max = 0.999329    avg =      nan  |avg| =      nan    to       -0
    14 insn x^5/12:
    Domain from 2.64698e-23
    error max =  0.13556    avg = 0.00125936    |avg| = 0.00354677  to 1.73472e-18
    error max = 0.000564988 avg = 2.51458e-06   |avg| = 0.000113709 to 1.13687e-13
    error max = 0.000565065 avg = -1.49258e-06  |avg| = 0.000112553 to 7.45058e-09
    error max = 0.000565143 avg = 1.5293e-06    |avg| = 0.000112864 to 0.000488281
    error max = 0.000565298 avg = 2.76457e-06   |avg| = 0.000113713 to       32
    error max = 0.000565453 avg = -1.61276e-06  |avg| = 0.000112561 to 2.09715e+06
    error max = 0.000565531 avg = 1.42628e-06   |avg| = 0.000112866 to 1.37439e+11
    error max = 0.000565686 avg = 2.71505e-06   |avg| = 0.000113715 to 9.0072e+15
    error max = 0.000565763 avg = -1.56586e-06  |avg| = 0.000112415 to 1.84467e+19
    

    我怀疑rsqrt操作限制了更准确的5/12的rsqrt


    另一个答案,因为这是与我以前的答案非常不同,这是快速的。 相对误差是3e-8。 想要更精确? 添加几个Chebychev条款。 最好保持奇数的顺序,因为这会导致2 ^ n-epsilon和2 ^ n + epsilon之间的小间断。

    #include <stdlib.h>
    #include <math.h>
    
    // Returns x^(5/12) for x in [1,2), to within 3e-8 (relative error).
    // Want more precision? Add more Chebychev polynomial coefs.
    double pow512norm (
       double x)
    {
       static const int N = 8;
    
       // Chebychev polynomial terms.
       // Non-zero terms calculated via
       //   integrate (2/pi)*ChebyshevT[n,u]/sqrt(1-u^2)*((u+3)/2)^(5/12)
       //   from -1 to 1
       // Zeroth term is similar except it uses 1/pi rather than 2/pi.
       static const double Cn[N] = { 
           1.1758200232996901923,
           0.16665763094889061230,
          -0.0083154894939042125035,
           0.00075187976780420279038,
          // Wolfram alpha doesn't want to compute the remaining terms
          // to more precision (it times out).
          -0.0000832402,
           0.0000102292,
          -1.3401e-6,
           1.83334e-7};
    
       double Tn[N];
    
       double u = 2.0*x - 3.0;
    
       Tn[0] = 1.0;
       Tn[1] = u;
       for (int ii = 2; ii < N; ++ii) {
          Tn[ii] = 2*u*Tn[ii-1] - Tn[ii-2];
       }   
    
       double y = 0.0;
       for (int ii = N-1; ii >= 0; --ii) {
          y += Cn[ii]*Tn[ii];
       }   
    
       return y;
    }
    
    
    // Returns x^(5/12) to within 3e-8 (relative error).
    double pow512 (
       double x)
    {
       static const double pow2_512[12] = {
          1.0,
          pow(2.0, 5.0/12.0),
          pow(4.0, 5.0/12.0),
          pow(8.0, 5.0/12.0),
          pow(16.0, 5.0/12.0),
          pow(32.0, 5.0/12.0),
          pow(64.0, 5.0/12.0),
          pow(128.0, 5.0/12.0),
          pow(256.0, 5.0/12.0),
          pow(512.0, 5.0/12.0),
          pow(1024.0, 5.0/12.0),
          pow(2048.0, 5.0/12.0)
       };
    
       double s;
       int iexp;
    
       s = frexp (x, &iexp);
       s *= 2.0;
       iexp -= 1;
    
       div_t qr = div (iexp, 12);
       if (qr.rem < 0) {
          qr.quot -= 1;
          qr.rem += 12;
       }
    
       return ldexp (pow512norm(s)*pow2_512[qr.rem], 5*qr.quot);
    }
    

    附录:这里发生了什么?
    根据请求,下面解释了上述代码的工作原理。

    概观
    上面的代码定义了两个函数, double pow512norm (double x)double pow512 (double x) 。 后者是套件的入口点; 这是用户代码应该调用来计算x ^(5/12)的函数。 函数pow512norm(x)使用Chebyshev多项式来逼近x ^(5/12),但仅限于范围[1,2]中的x。 (使用pow512norm(x)作为该范围之外的x的值,结果将是垃圾。)

    函数pow512(x)将传入x成一对(double s, int n)使得x = s * 2^n并且使得1≤ s <2。 进一步将n划分为(int q, unsigned int r)使得n = 12*q + rr小于12,可以将x ^(5/12)的问题分解为几部分:

  • 通过(uv)^ a =(u ^ a)(v ^ a)得到正数x^(5/12)=(s^(5/12))*((2^n)^(5/12))你,真实的。
  • s^(5/12)通过pow512norm(s)计算。
  • (2^n)^(5/12)=(2^(12*q+r))^(5/12)
  • 通过u^(a+b)=(u^a)*(u^b)计算出正数u, 2^(12*q+r)=(2^(12*q))*(2^r) A,b。
  • (2^(12*q+r))^(5/12)=(2^(5*q))*((2^r)^(5/12))
  • (2^r)^(5/12)由查找表pow2_512计算。
  • 计算pow512norm(s)*pow2_512[qr.rem] ,我们pow512norm(s)*pow2_512[qr.rem] 。 这里qr.rem是上面步骤3中计算出的r值。 所需要的就是将其乘以2^(5*q)以产生期望的结果。
  • 这正是数学库函数ldexp
  • 函数逼近
    这里的目标是想出一个易于计算的f(x)= x ^(5/12)的近似值,对于手头的问题来说,'足够好'。 我们的近似值在某种意义上应该接近f(x)。 修辞问题:“接近”是什么意思? 两种相互竞争的解释是最小化均方误差与最大绝对误差的最小化。

    我将用股市类比来描述这些差异。 假设你想为你的最终退休储蓄。 如果你二十多岁,最好的办法是投资股票或股票市场基金。 这是因为在足够长的时间内,股票市场平均击败任何其他投资计划。 然而,我们都看到过把资金投入股票是一件非常糟糕的事情。 如果你五六十岁(或四十岁,如果你想退休),你需要更保守地投资。 这些下滑可能会对你的退休投资组合产生影响。

    回到函数逼近:作为一些近似的消费者,您通常担心最糟糕的错误,而不是“平均”的性能。 使用一些近似值来构建“平均”最佳性能(例如最小平方),墨菲定律指出,如果性能远远低于平均水平,您的程序将花费大量时间使用近似值。 你想要的是一个极小极大值近似值,它可以使某个域的最大绝对误差最小化。 一个好的数学图书馆将采用极大极小的方法,而不是最小二乘法,因为这可以让数学图书馆的作者提供一定的保证性能。

    数学库通常使用多项式或有理多项式来近似某个域a≤x≤b上的某个函数f(x)。 假设函数f(x)是在这个域上解析的,并且你想用函数N的某个多项式p(x)来近似函数N.对于给定的度N,存在一些神奇的唯一多项式p(x),使得p x)-f(x)在[a,b]上有N + 2个极值,并且这些N + 2极值的绝对值都相等。 寻找这个神奇的多项式p(x)是函数逼近器的圣杯。

    我没有为你找到圣杯。 我改用Chebyshev近似。 第一类Chebyshev多项式是一个正交(但不正交)的多项式集,当涉及到函数逼近时,它具有一些非常好的特性。 切比雪夫近似常常非常接近那个魔法多项式p(x)。 (实际上,确实发现圣杯多项式的Remez交换算法通常以切比雪夫近似开始。)

    pow512norm(x)的
    该函数使用切比雪夫近似来找出近似x ^(5/12)的多项式p *(x)。 在这里,我使用p *(x)来区分这个Chebyshev近似值和上面描述的神奇多项式p(x)。 Chebyshev近似p *(x)很容易找到; 发现p(x)是一只熊。 Chebyshev逼近p *(x)是sum_i Cn [i] * Tn(i,x),其中Cn [i]是Chebyshev系数,Tn(i,x)是在x处评估的Chebyshev多项式。

    我用Wolfram alpha来为我找到Chebyshev系数Cn 。 例如,这会计算Cn[1] 。 输入框之后的第一个框具有所需的答案,在这种情况下为0.166658。 这不是我想要的那么多数字。 点击'更多数字',瞧,你会得到更多的数字。 Wolfram alpha是免费的; 它会做多少计算是有限制的。 它在更高阶的条件下达到了上限。 (如果您购买或可以访问数学,您将能够高度精确地计算这些高阶系数。)

    Chebyshev多项式Tn(x)在阵列Tn中计算。 除了给出非常接近神奇多项式p(x)的东西之外,使用切比雪夫近似的另一个原因是这些切比雪夫多项式的值很容易计算:从Tn[0]=1并且Tn[1]=x ,然后迭代地计算Tn[i]=2*x*Tn[i-1] - Tn[i-2] 。 (我在代码中使用'ii'作为索引变量而不是'i',我从不使用'i'作为变量名称。连续'我的?)

    pow512(x)的
    pow512是用户代码应该调用的函数。 我已经在上面描述了这个函数的基础。 一些更多的细节:数学库函数frexp(x)返回尾数s和指数iexp为输入x 。 (次要问题:我希望s在1和2之间用于pow512normfrexp返回一个介于0.5和1之间的值)。数学库函数div在一个swell foop中返回整数除法的商和余数。 最后,我使用数学库函数ldexp将三个部分放在一起形成最终答案。


    伊恩斯蒂芬森写这个代码,他声称优于pow() 。 他将这个想法描述如下:

    Pow基本上是用log's: pow(a,b)=x(logx(a)*b) 。 所以我们需要一个快速的日志和快速指数 - 它并不重要,所以我们使用的是2.诀窍是浮点数已经是日志格式的格式:

    a=M*2E
    

    取双方的日志给出:

    log2(a)=log2(M)+E
    

    或者更简单:

    log2(a)~=E
    

    换句话说,如果我们取一个数字的浮点数表示,并提取Exponent,那么我们就有一个很好的起点作为其日志。 事实证明,当我们通过按摩位模式来做到这一点时,尾数最终给出了一个很好的近似误差,并且它工作得很好。

    这对于简单的照明计算应该足够好,但如果您需要更好的东西,则可以提取尾数,然后使用它计算相当准确的二次修正系数。

    链接地址: http://www.djcxy.com/p/85721.html

    上一篇: Optimizations for pow() with const non

    下一篇: Accurate C/C++ clock on a multi