在浮點值中獲得平方根的最快方法-有解無憂

我試圖找到一種最快的方法來計算 C 中任何浮點數的平方根。我在一個巨大的粒子運動計算中使用這種型別的函式，比如計算兩個粒子之間的距離，我們需要一個平方根等。所以如果有任何建議會非常有幫助。我試過了，下面是我的代碼

#include <math.h>
#include <iostream>
#include <chrono>

using namespace std;
using namespace std::chrono;

#define CHECK_RANGE 100

inline float msqrt(float a)
{
    int i;
    for (i = 0;i * i <= a;i  );
    
    float lb = i - 1; //lower bound
    if (lb * lb == a)
        return lb;
    float ub = lb   1; // upper bound
    float pub = ub; // previous upper bound
    for (int j = 0;j <= 20;j  )
    {
        float ub2 = ub * ub;
        if (ub2 > a)
        {
            pub = ub;
            ub = (lb   ub) / 2; // mid value of lower and upper bound
        }
        else
        {
            lb = ub; 
            ub = pub;
        }
    }
    return ub;
}

void check_msqrt()
{
    for (size_t i = 0; i < CHECK_RANGE; i  )
    {
        msqrt(i);
    }
}

void check_sqrt()
{
    for (size_t i = 0; i < CHECK_RANGE; i  )
    {
        sqrt(i);
    }
}

int main()
{
    auto start1 = high_resolution_clock::now();
    check_msqrt();
    auto stop1 = high_resolution_clock::now();

    auto duration1 = duration_cast<microseconds>(stop1 - start1);
    cout << "Time for check_msqrt = " << duration1.count() << " micro secs\n";


    auto start2 = high_resolution_clock::now();
    check_sqrt();
    auto stop2 = high_resolution_clock::now();

    auto duration2 = duration_cast<microseconds>(stop2 - start2);
    cout << "Time for check_sqrt = " << duration2.count() << " micro secs";
    
    //cout << msqrt(3);

    return 0;
}

上述代碼的輸出顯示實作的方法比 math.h 檔案的 sqrt 慢 4 倍。我需要比 math.h 版本更快的版本。在浮點值中獲得平方根的最快方法

uj5u.com熱心網友回復：

簡而言之，我認為不可能實作比標準庫版本更快的東西sqrt。

在實作標準庫函式時，性能是一個非常重要的引數，可以公平地假設這樣一個常用的函式盡可能sqrt地進行了優化。

擊敗標準庫函式需要特殊情況，例如：

在標準庫尚未專門針對的特定系統上提供合適的匯編指令 - 或其他專門的硬體支持。
所需范圍或精度的知識。標準庫函式必須處理標準指定的所有情況。如果應用程式只需要其中的一個子集，或者可能只需要一個近似結果，那么也許可以進行優化。
對計算進行數學簡化或以智能方式組合一些計算步驟，以便可以為該組合進行有效實施。

uj5u.com熱心網友回復：

這是二進制搜索的另一種替代方法。它可能沒有那么快std::sqrt，還沒有測驗過。但它肯定會比你的二進制搜索更快。

auto
Sqrt(float x)
{
    using F = decltype(x);
    if (x == 0 || x == INFINITY || isnan(x))
        return x;
    if (x < 0)
        return F{NAN};
    int e;
    x = std::frexp(x, &e);
    if (e % 2 != 0)
    {
          e;
        x /= 2;
    }
    auto y = (F{-160}/567*x   F{2'848}/2'835)*x   F{155}/567;
    y = (y   x/y)/2;
    y = (y   x/y)/2;
    return std::ldexp(y, e/2);
}

在排除 /-0、nan、inf 和負數之后，它通過將分解為 [ ¹ / ₄ , 1) 乘以 2 ^efloat范圍內的尾數來作業，其中是偶數。答案是 sqrt(mantissa)* 2 ^{^e}^/^₂。e

^{找到尾數的 sqrt 可以通過在 [ 1} / ₄ , 1]范圍內擬合的最小二乘二次曲線來猜測。然后通過牛頓-拉夫森的兩次迭代來完善這個好的猜測。這將使您在正確舍入結果的1 ulp范圍內。一個好的std::sqrt通常會得到最后一點正確。

uj5u.com熱心網友回復：

我也嘗試過https://en.wikipedia.org/wiki/Fast_inverse_square_root中提到的演算法，但沒有找到想要的結果，請檢查

#include <math.h>
#include <iostream>
#include <chrono>

#include <bit>
#include <limits>
#include <cstdint>

using namespace std;
using namespace std::chrono;

#define CHECK_RANGE 10000

inline float msqrt(float a)
{
    int i;
    for (i = 0;i * i <= a;i  );
    
    float lb = i - 1; //lower bound
    if (lb * lb == a)
        return lb;
    float ub = lb   1; // upper bound
    float pub = ub; // previous upper bound
    for (int j = 0;j <= 20;j  )
    {
        float ub2 = ub * ub;
        if (ub2 > a)
        {
            pub = ub;
            ub = (lb   ub) / 2; // mid value of lower and upper bound
        }
        else
        {
            lb = ub; 
            ub = pub;
        }
    }
    return ub;
}

/* mentioned here ->  https://en.wikipedia.org/wiki/Fast_inverse_square_root */
inline float Q_sqrt(float number)
{
    union Conv {
        float    f;
        uint32_t i;
    };
    Conv conv;
    conv.f= number;
    conv.i = 0x5f3759df - (conv.i >> 1);
    conv.f *= 1.5F - (number * 0.5F * conv.f * conv.f);
    return 1/conv.f;
}

void check_Qsqrt()
{
    for (size_t i = 0; i < CHECK_RANGE; i  )
    {
        Q_sqrt(i);
    }
}

void check_msqrt()
{
    for (size_t i = 0; i < CHECK_RANGE; i  )
    {
        msqrt(i);
    }
}

void check_sqrt()
{
    for (size_t i = 0; i < CHECK_RANGE; i  )
    {
        sqrt(i);
    }
}

int main()
{
    auto start1 = high_resolution_clock::now();
    check_msqrt();
    auto stop1 = high_resolution_clock::now();

    auto duration1 = duration_cast<microseconds>(stop1 - start1);
    cout << "Time for check_msqrt = " << duration1.count() << " micro secs\n";


    auto start2 = high_resolution_clock::now();
    check_sqrt();
    auto stop2 = high_resolution_clock::now();

    auto duration2 = duration_cast<microseconds>(stop2 - start2);
    cout << "Time for check_sqrt = " << duration2.count() << " micro secs\n";
    
    auto start3 = high_resolution_clock::now();
    check_Qsqrt();
    auto stop3 = high_resolution_clock::now();

    auto duration3 = duration_cast<microseconds>(stop3 - start3);
    cout << "Time for check_Qsqrt = " << duration3.count() << " micro secs\n";

    //cout << Q_sqrt(3);
    //cout << sqrt(3);
    //cout << msqrt(3);
    return 0;
}

轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/438738.html

標籤：C c 11 平方数学.sqrt

上一篇：C 11和C 17基于范圍的For回圈可以迭代到特定位置而不是地圖的整個范圍嗎？

下一篇：將std::string添加到具有記憶體中c-str緩沖區值的STL容器的有效方法是什么？