This site探索了几种在C中有效计算立方根的计算方法，并提供了一些源代码供您下载。
（编辑：在谷歌上搜索“快速立方根”，会出现几个看起来很有希望的结果。）
立方根是一个有趣的主题，因为它们在许多常见的公式中使用，而Microsoft Visual Studio中没有快速的立方根函数。
在没有特殊的立方根函数的情况下，一个典型的策略是通过幂函数（例如，pow（x，1.0/3.0））进行计算。这在速度方面可能存在问题，并且在负数处理不当时，在准确性方面可能存在问题。
他的网站上有一些关于所用方法的基准测试，它们都比pow()快得多。

32-bit float tests
----------------------------------------
cbrt_5f      8.8 ms    5 mbp   6.223 abp
pow        144.5 ms   23 mbp  23.000 abp
halley x 1  31.8 ms   15 mbp  18.961 abp
halley x 2  59.0 ms   23 mbp  23.000 abp
newton x 1  23.4 ms   10 mbp  12.525 abp
newton x 2  48.9 ms   20 mbp  22.764 abp
newton x 3  72.0 ms   23 mbp  23.000 abp
newton x 4  89.6 ms   23 mbp  23.000 abp

字符串
请参阅网站下载源。

下面的实现比AVX-512 CPU上的std：：pow快4倍，具有相对较高的容差（0.000001）。它由垂直自动向量化循环组成，用于每个基本操作，如乘法和除法，以便它一次计算8，16，32个元素，而不是水平向量化Newton-Raphson循环。

#include <cmath> 
/* 
   Newton-Raphson iterative solution
   f_err(x) = x*x*x - N
   f'_err(x) = 3*x*x
   x = x - (x*x*x - N)/(3*x*x)
   x = x - (x - N/(x*x))/3   <--- repeat until error < tolerance
   but with vertical-parallelization 
*/
template < typename Type, int Simd, int inverseTolerance> 
    inline 
void cubeRootFast(Type *const __restrict__ data, 
        Type *const __restrict__ result) noexcept 
{ 
    // alignment 64 required for AVX512 vectorization
    alignas(64) 
    Type xd[Simd]; 

    alignas(64) 
    Type resultData[Simd]; 

    alignas(64) 
    Type xSqr[Simd]; 

    alignas(64) 
    Type nDivXsqr[Simd]; 

    alignas(64) 
    Type diff[Simd]; 

    // cube root checking mask
    for (int i = 0; i < Simd; i++) 
    { 
        xd[i] = data[i] <= Type(0.000001);
    } 

    // skips division by zero if input is zero or close to zero
    for (int i = 0; i < Simd; i++) 
    { 
        resultData[i] = xd[i] ? Type(1.0) : data[i]; 
    } 

    // Newton-Raphson Iterations in parallel 
    bool work = true; 
    while (work) 
    { 
        // compute x*x
        for (int i = 0; i < Simd; i++) 
        { 
            xSqr[i] = resultData[i] *resultData[i]; 
        } 

        // compute N/(x*x)
        for (int i = 0; i < Simd; i++) 
        { 
            nDivXsqr[i] = data[i] / xSqr[i]; 
        } 

        // compute x - N/(x*x)
        for (int i = 0; i < Simd; i++) 
        { 
            nDivXsqr[i] = resultData[i] - nDivXsqr[i]; 
        } 

        // compute (x-N/(x*x))/3
        for (int i = 0; i < Simd; i++) 
        { 
            nDivXsqr[i] = nDivXsqr[i] / Type(3.0); 
        } 

        // compute x - (x-N/(x*x))/3
        for (int i = 0; i < Simd; i++) 
        { 
            diff[i] = resultData[i] - nDivXsqr[i]; 
        } 

        // compute error
        for (int i = 0; i < Simd; i++) 
        { 
            diff[i] = resultData[i] - diff[i]; 
        } 

        // compute absolute error
        for (int i = 0; i < Simd; i++) 
        { 
            diff[i] = std::abs(diff[i]); 
        } 

        // compute condition to stop looping (error < tolerance)?
        for (int i = 0; i < Simd; i++) 
        { 
            diff[i] = diff[i] > Type(1.0/inverseTolerance); 
        } 

        // all SIMD lanes have to have zero work left to end
        Type check = 0; 
        for (int i = 0; i < Simd; i++) 
        { 
            check += diff[i]; 
        } 
        work = (check > Type(0.0)); 

        // compute the next x guess
        for (int i = 0; i < Simd; i++) 
        { 
            resultData[i] = resultData[i] - nDivXsqr[i]; 
        } 
    } 

    // if input was close to zero, output zero
    // output result otherwise
    for (int i = 0; i < Simd; i++) 
    { 
        result[i] = xd[i] ? Type(0.0) : resultData[i]; 
    } 
} 
#include <iostream> 
 
int main() 
{ 
    constexpr int n = 8192; 
    constexpr int simd = 16; 
    constexpr int inverseTolerance = 1000;
    float data[n]; 
    for (int i = 0; i < n; i++) 
    { 
        data[i] = i; 
    } 
    for (int i = 0; i < n; i += simd) 
    { 
        cubeRootFast<float, simd, inverseTolerance> (data + i, data + i); 
    } 
    for (int i = 0; i < 10; i++) 
        std::cout << data[i *i *i] << std::endl; 
    return 0; 
}

字符串
它只在GCC上测试，所以它可能需要在每个循环上使用额外的MSVC-pragmas来强制自动向量化。如果你有OpenMP，那么你也可以使用#pragma omp simd safelen(Simd)来实现同样的事情。
性能仅在[0，1]范围内保持。要使用更大的值，您应该像这样使用范围缩小：

// example: max value is 1000
for(auto & input:inputs)
    input = input/1000.0f // normalize

for(..)
    cubeRootFast<float, simd, inverseTolerance> (input + i, input + i)

 for(auto & input:inputs)
    input = 10.0f*input // de-normalize  (1000 = 10 x 10 x 10)

型
如果你只需要在低范围（如[0，1000]）上有0.005的错误，并具有16倍的加速，你可以尝试下面使用多项式近似的实现（Horner Scheme应用于使用FMA指令进行计算，并且不需要显式的自动向量化，因为它不包括任何分支/循环）：

// optimized for [0,1] range: ~1 cycles on AVX512, 0.003 average error 
// polynomial approximation with Horner Scheme for FMA optimization 
template<typename T> 
T cubeRootFast(T x) 
{ 
 T xd = x-T(1.0); 
 T result = T(-55913.0/4782969.0); 
 result *= xd; 
 result += T(21505.0/1594323.0);  
 result *= xd; 
 result += T(-935.0/59049.0);  
 result *= xd; 
 result += T(374.0/19683.0);  
 result *= xd; 
 result += T(-154.0/6561.0);  
 result *= xd; 
 result += T(22.0/729.0);  
 result *= xd; 
 result += T(-10.0/243.0); 
 result *= xd; 
 result += T(5.0/81.0); 
 result *= xd; 
 result += T(-1.0/9.0);  
 result *= xd; 
 result += T(1.0/3.0); 
 result *= xd; 
 result += T(1.0); 
 return result; 
} 
 
// range reduction + dereduction: ~ 1 cycles on AVX512 
 for(int i=0;i<8192;i++) 
 {  
     float inp = input[i]; 
     // scaling + descaling for range [1,999] 
     float scaling = (inp>333.0f)?(1000.0f):(333.0f); 
     scaling = (inp>103.0f)?scaling:(103.0f); 
     scaling = (inp>29.0f)?scaling:(29.0f); 
     scaling = (inp>7.0f)?scaling:(7.0f); 
     scaling = (inp>3.0f)?scaling:(3.0f); 
      
     output[i] = powf(scaling,0.33333333333f)*cubeRootFast<float>(inp/scaling);  
 }

型

MSC 2012一定是他们的最后一个版本，没有在math.h中实现cbrt。MS系统cbrt基准测试非常糟糕，但在真实的代码中并不差。尽管它不是最准确的。
我发现在大多数编译器上既准确又快速的cbrt公共代码的最佳组合是布鲁斯埃文斯的Kahan的魔术常数算法的Sun's implementation。英特尔在他们的2023编译器中的系统cbrt通过在最终改进中小心使用截断和FMA而令人难以置信地准确。
如果没有更好的选择，一个简单的选择是：

double cbrt(double x)
{  // fast and accurate. NR isn't sufficient to tidy it up
   // could be improved by precision (y*y)*y-x using FMA
   double t, y;
   if (x) y = exp(log(abs(x)) / 3); else return x;
   if (x<0) y = -y;
   t = y*y*y;
   return y - y * (t - x) / (2 * t + x);  // Halley refinement
}

字符串
根据CPU架构的不同，这两种方法都值得一试。