5/30 一對一討論 問題回答

contributed by < willy-liu >

Q1: 不使用 FPU 指令進行浮點數開平方根

Assume x is a single precision floating point number and always >= 1. Calculate square root (casting to int for easing complexity) without FPU, as precise as FPU operations.

int my_sqrtf(float x){}

首先注意到為了讓題目可以在十分鐘內回答出來，老師簡化了回傳值，並且限定輸入>=1，而我當下太在意小數點的精確度，而回答了 IEEE 754 的浮點數規格，然後打算用 LUT 的方式確保小數點，但這題若只要回傳"整數"，我們就需要確認精確度的問題。

二分搜尋法

先考慮計算無條件捨去的情況
在無法使用 FPU 的情況，我只能使用int和位元運算來實現，最基本的方式是透過二分搜尋法來逼近。

int my_sqrtf(float x) {
    const int num = (int)x;

    if (num == 1)
        return num;

    int low = 1, high = num >> 2 + 1, ans = 1;
    while (low <= high) {
        int mid = low + ((high - low) >> 1);
        if (mid <= num / mid) {
            ans = mid;
            low = mid + 1;
        } else {
            high = mid - 1;
        }
    }
    return ans;
}

p.s. 原本條件為mid * mid <= num，避免溢位改成mid <= num / mid

因為是二分搜尋法，所以迴圈要跑

$\log (x)$次，執行時間會隨著值上升變久，而且中間需要用到除法（如果要小數後一位四捨五入可以在最後判斷

$(ans\ast 10+5{)}^2$是否

$\le x\ast 100$
if ((ans*10+5) * (ans*10+5) <= num * 100) ans++;）

直式開方法

考慮更快的方法，直式開方法（又叫做digit-by-digit法），計算平方根時我們就是要找到一個最大的

$y$使

$y^2\le x$，假設a是b的前一位，我們可以將

$(y^2)$視為

$(a+b{)}^2$，改寫式子為

$a^2+(2a+b)b$，至此，就可以從高位數往低位數去逼近平方根。

10進位範例：
x = 6242 =

$a^2+(2a+b)b$
由於y的每位數會貢獻x的2位數，所以一次要考慮兩位數。
首先，找到最大的a使

$a^2\le 62$，得出a = 7
故x = 6242 =

$7{?}^2$ =

${70}^2+(2\ast 70+b)b$
首兩位"62"由於a的貢獻只會剩下

$62-7^2=13$
接著再往後考慮兩位，1342會是

$(2\ast 70+b)b$的貢獻，得出b=9
1342經過b的貢獻會剩下1342-

$(2\ast 70+9)9=1342-1341=1$
根據這個公式可以再無限往下擴展，例如c為b的後一位數
一樣再往下取兩位數100，會是

$[2\times 79\times 10+c]\times c=(1580+c)\times c$的貢獻，得出c=0
…

$N_k=10N_{k-1}+d_k$，

$R_k=R_k^{\prime }-(20N_{k-1}+d_k)d_k$
其中

$N_{k-1}$ 是前

$k-1$ 位組合成的數字（不含小數點，視為整數）

$R_{k-1}$ 是當前被減數

$d_k$ 是本輪要決定的新數字

參考資料：https://zh.wikipedia.org/wiki/平方根#長除式算法

套用到二進制，每一輪就只要決定新數字要是0/1就好

int my_sqrtf_bitwise(float x) {
    int num = (int)x; 
    int res = 0;
    int bit = 1 << 30; // The second-to-top bit is set

    // "bit" starts at the highest power of four <= x
    while (bit > num) bit >>= 2;

    while (bit != 0) {
        if (num >= res + bit) { // decide b = 1 or 0
            num -= res + bit;
            res = (res >> 1) + bit;
        } else {
            res >>= 1; 
        }
        bit >>= 2;
    }
    return res;
}

這個方法僅使用位元運算和加減法，而且無論輸入是多少都只要進行16次的迴圈，如果要考慮小數後一位，一樣可以在最後透過位移，再多算一次小數後一位的情況，如果小數後一位的位元是"1"，代表小數部分

$\ge$ 0.5，所以回傳的int就要加一。

    num <<= 2;
    res <<= 2;
    bit = 1;
    
    if (num >= res + bit) 
        res = (res>>1) + bit;
    else 
        res >>= 1;    

這段程式碼其實也等價於res += (num >= res + 1 ? 1 : 0);，因為僅考慮下一位帶來的影響，做位移僅僅是為了好讀。

時間分析

Binary-search   mean time: 75.09 ns
Bitwise         mean time: 45.73 ns (1.64× faster)

誤差分析

test range: [1.0, 1000.0], step: 0.2
without rounding：  mean error = 3.003131%, max error = 48.360222% (x = 3.8)
with rounding：     mean error = 1.552714%, max error = 39.697735% (x = 2.8)

在x較小時誤差會比較多，最大高達40~50%。

Q2: following Q1, the return type is float. 誤差 GPU

由於要讓回傳值也是float，首先要了解IEEE 754中如何定義一個single precision floating point

由定義可以知道一個 IEEE 754 的單精度浮點數的值表示法如下

$$\text{value}=(-1{)}^{\text{sign}}\times 2^{(E-127)}\times (1+\sum _{i=1}^{23}{b}_{23-i}{2}^{-i})$$

直式開方法

在完成「四捨五入到小數點第一位」的處理後，我想到可以把這個方法向更高精度推進，並且將小數部分也一併納入計算，

$\mathrm{以}\sqrt{12300}\text{為例}\sqrt{12300}=\sqrt{1.23\ast 10000}=\sqrt{1.23}\ast 100$也就是說，只要針對有效數字（即 1.23）計算開方，然後再將結果乘以對應的指數位移，就能得出最終值。若將這個思路應用到浮點運算中，只需對 fraction 部分進行開方計算，計算完成後再把 exponent 按原先偏移量調整即可。

$$\begin{array}{rl}& \text{令 }x=(-1{)}^{\text{sign}}{2}^{E-127}(1+\sum _{i=1}^{23}{b}_{23-i}{2}^{-i}),\\ & e=E-127=2k+\pi ,\pi \in \{0,1\},k=\lfloor \frac{e}{2}\rfloor ,\\ & \text{定義有效小數部分 }A=2^{\pi }(1+\sum _{i=1}^{23}{b}_{23-i}{2}^{-i}),\\ & \sqrt{x}=(-1{)}^{\text{sign}}{2}^k\sqrt{A},\\ & \text{令整數尾數 }{M}_{\mathrm{i}\mathrm{n}\mathrm{t}}=A\times 2^{23},\\ & B=M_{\mathrm{i}\mathrm{n}\mathrm{t}}\times 2^7=A,B\in [2^{30},2^{32}],\\ & S=\lfloor \sqrt{B}\rfloor =\lfloor \sqrt{A}{2}^{15}\rfloor ,\\ & \text{取 23 位： }\tilde{M}=S\times 2^8,\tilde{m}=\frac{\tilde{M}}{2^{23}}=\frac{S}{2^{15}}=1+\sum _{i=1}^{23}{d}_{23-i}{2}^{-i},\\ & E^{\prime }=k+127,\\ & \hat{\sqrt{x}}=(-1{)}^{\text{sign}}{2}^{E^{\prime }-127}(1+\sum _{i=1}^{23}{d}_{23-i}{2}^{-i}).\end{array}$$

float my_sqrtf_bitwise(float x) {
    unsigned int i;
    memcpy(&i, &x, sizeof(i));
    if (x == 0.0f) return 0.0f;
    if (i & 0x80000000u) return -1.0f;

    unsigned int e = ((i >> 23) & 0xFF) - 127;
    unsigned int m = (i & 0x7FFFFFu) | 0x800000u;

    /* 合併奇偶指數處理並放大供整數開方 */
    unsigned int scaled = ((e & 1) ? m << 1 : m) << 7;
    e = (e - (e & 1)) / 2;

    unsigned int s = my_sqrtf_bitwise_int(scaled);
    unsigned int nm = s << 8;

    unsigned int ne = e + 127;
    if (ne >= 0xFF) {
        i = 0x7F800000u;        // overflow
    } else if (ne <= 0) {
        i = 0x00000000u;        // underflow
    } else {
        i = ((unsigned)ne << 23) | (nm & 0x7FFFFFu);
    }

    memcpy(&x, &i, sizeof(x));
    return x;
}

這邊要注意的是 exponent 要分為奇數和偶數考慮，奇數時要借位 mantissa （改進後使用位元運算避免分支），此時　mantissa　可能是24 或 25 bits(包含IEEE 754省略的1），為了提高計算的精度，將其<<7再進行直式開方（要記得把my_sqrtf_bitwise輸入、回傳改成unsigned int），才能有效利用完整的unsigned int，經過直式開方的值會剩下16bit，再<<8變回 IEEE 754 的 fraction 格式，最後再進行float的重組。

如果要再改善精度，完整利用float的fraction部分，我們可以考慮將直式開方的輸入、輸出改為long long，然後從<<46的 bit 開始計算， mantissa 傳入前也先轉型並<<23，就可以完整利用 fraction 的23個 bits ，精確度也會提升到跟sqrtf幾乎一致，但在測試後發現 mantissa 的最後一個 bit 可能差1。
為了找出原因，我查閱了sqrtf的 man page，雖然其中沒有說明當數值無法完整表示時的處理方式，但可以確定它會遵循 IEEE 754 的規範。 rounding 的相關摘要如下：

4.3.1 Rounding-direction attributes to nearest

In the following two rounding-direction attributes, an infinitely precise result with magnitude at least

$b^{emax}\times (b–½b^{1–p})$ shall round to ±∞ with no change in sign; here e_max and p are determined by the destination format (see 3.3). With:

• roundTiesToEven, The floating-point number nearest to the infinitely > precise result shall be delivered. If the two nearest floating-point numbers bracketing an unrepresentable infinitely precise result are equally near, the one with an even least-significant digit shall be delivered. If that is not possible, the one larger in magnitude shall be delivered. (This can happen for one-digit precision—for example, when rounding 9.5 to one digit both 9 and 1×10^1 have odd significands.)
• roundTiesToAway, The floating-point number nearest to the infinitely precise result shall be delivered. If the two nearest floating-point numbers bracketing an unrepresentable infinitely precise result are equally near, the one with larger magnitude shall be delivered.

很顯然，在sqrtf的情況需要再額外考慮下一位的值進行 rounding ，因此，在加入res += (num >= res + 1 ? 1 : 0);後得到的結果就與sqrtf完全一致。

-    unsigned int scaled_input_for_isqrt = input_for_isqrt_base << 7; //放大輸入以提高整數開方的精度
-    unsigned int sqrt_intermediate_result = my_sqrtf_bitwise(scaled_input_for_isqrt);
-    unsigned int new_mantissa_int_val = sqrt_intermediate_result << 8;
+    long long scaled_input_for_isqrt = (long long)input_for_isqrt_base << 23;
+    long long sqrt_intermediate_result = integer_sqrt_longlong(scaled_input_for_isqrt);
+    unsigned int new_mantissa_int_val = sqrt_intermediate_result;

Comparison with math.h sqrtf:
sqrt(2.00000000e+00):
  Custom: 1.41421354e+00(0x3fb504f3)
  math.h: 1.41421354e+00(0x3fb504f3)
sqrt(3.00000000e+00):
  Custom: 1.73205078e+00(0x3fddb3d7)
  math.h: 1.73205078e+00(0x3fddb3d7)
sqrt(9.00000000e+00):
  Custom: 3.00000000e+00(0x40400000)
  math.h: 3.00000000e+00(0x40400000)
sqrt(5.00000000e-01):
  Custom: 7.07106769e-01(0x3f3504f3)
  math.h: 7.07106769e-01(0x3f3504f3)
sqrt(2.50000000e-01):
  Custom: 5.00000000e-01(0x3f000000)
  math.h: 5.00000000e-01(0x3f000000)
sqrt(1.60000000e+01):
  Custom: 4.00000000e+00(0x40800000)
  math.h: 4.00000000e+00(0x40800000)
sqrt(1.00000000e+00):
  Custom: 1.00000000e+00(0x3f800000)
  math.h: 1.00000000e+00(0x3f800000)
sqrt(1.23456777e+04):
  Custom: 1.11111107e+02(0x42de38e3)
  math.h: 1.11111107e+02(0x42de38e3)
sqrt(1.23399996e-05):
  Custom: 3.51283350e-03(0x3b663791)
  math.h: 3.51283350e-03(0x3b663791)
sqrt(1.17549435e-38):
  Custom: 1.08420217e-19(0x20000000)
  math.h: 1.08420217e-19(0x20000000)
sqrt(3.40282347e+38):
  Custom: 1.84467430e+19(0x5f7fffff)
  math.h: 1.84467430e+19(0x5f7fffff)

牛頓法(有fpu)

牛頓法的概念就是透過微分找到切線

$f^{\prime }(x_n)=0$的

$x$值，然後用新的

$x_{n+1}$繼續迭代去逼近這個可微函數的解。
可以從函數在當前點

$x_n$處的一階泰勒展開（切線近似）出發：

$$f(x_{n+1})\approx f(x_n)+f^{\prime }(x_n)(x_{n+1}-x_n).$$
令上述切線在 (y=0) 處求交點：

$$\begin{array}{r}0=f(x_n)+f^{\prime }(x_n)(x_{n+1}-x_n)\\ f^{\prime }(x_n)(x_{n+1}-x_n)=-f(x_n)\\ x_{n+1}-x_n=-\frac{f(x_n)}{f^{\prime }(x_n)}\\ x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}.\end{array}$$

另外，從上圖可以發現

$x_0$的選擇很重要，如果距離解太遠的話會需要迭代很多次，如果一開始

$x_0$就離解很近，也許只需要迭代1~2次就可以與解的誤差很小，通常會使用

$a$或是

$\frac{a}{2}$出發。

由於我要計算的是

$\sqrt{a}$，實際上要找的就是

$f(x)=x^2-a$的解，帶入公式後就變成

$$x_{n+1}=x_n-\frac{(x_n^2-a)}{2x_n}=x_n-\frac{1}{2}(x_n-\frac{a}{x_n})$$

float Q_sqrt(float number , int iter)
{
    float y = number;
    while(iter-- > 0) {
        y  =  y - 0.5f * (y - number / y);
    }
    return y;
}

接著我們來思考如何去找到一個更好的

$x_0$，由前面知道float的表示法為(簡化 mantissa 部分以及省略 sign bit)

$$\begin{array}{rl}a& =2^{(E-127)}\times (1+m),0\le m<1\\ \sqrt{a}& =2^{\frac{1}{2}(E-127)}\times (1+m{)}^{\frac{1}{2}},1\le (1+m{)}^{\frac{1}{2}}<1.414...\end{array}$$

$$\mathrm{設}\sqrt{a}=2^{(E^{\prime }-127)}\times (1+m^{\prime }),1\le (1+m^{\prime })<2$$

mantissa 的部分會很接近，而對於指數部分我們取近似

$$\begin{array}{r}\frac{1}{2}(E-127)\approxeq (E^{\prime }-127)\\ E^{\prime }\approxeq 63.5+\frac{1}{2}E\end{array}$$
在 IEEE-754 二進制浮點裡，加上 63 的偏移量就對應到常數位元

• 63 → 十六進位 0x1F800000
• 63.5 → 十六進位 0x1FC00000，其中那個「.5」正好是把隱藏位之後的第一個 fraction bit 設成 1

然而在測試後注意到由於我們同時會對尾數做右移

$\lfloor \frac{(E-127)}{2}\rfloor$ 的操作，實際上那個「.5」加到 exponent 上之後，不會和原本的 mantissa bits 重疊覆蓋

+    uint32_t i;
    float y = number;
+    memcpy(&i, &y, sizeof(i));
+    i  = 0x1F800000 + ( i >> 1 ); // or use 0x1FC00000
+    memcpy(&y, &i, sizeof(y));

誤差分析（ magic number 使用0x1fc00000）：

在0.01~100取1000個點測試，每個點測試1000次。

Error percentages for each iteration:
Iteration 1:
  Plain: Max Error = 405.000000%, Mean Error = 243.526647%
  Magic: Max Error = 0.172270%, Mean Error = 0.023569%
Iteration 2:
  Plain: Max Error = 162.401000%, Mean Error = 88.811777%
  Magic: Max Error = 0.000159%, Mean Error = 0.000012%
Iteration 3:
  Plain: Max Error = 50.255300%, Mean Error = 23.026584%
  Magic: Max Error = 0.000036%, Mean Error = 0.000002%
Iteration 4:
  Plain: Max Error = 8.404350%, Mean Error = 2.845740%
  Magic: Max Error = 0.000036%, Mean Error = 0.000002%

可以看到加了初始值的預估後只做一次牛頓法就可以非常接近正確值，甚至不使用牛頓法，單純預估值的平均誤差大約在1.6%，最大誤差僅6%左右（Max Error = 6.044607%, Mean Error = 1.555732%），在一些場景下這樣的精確度也許就已經足夠。

相較之下若將

$x_0=a$的話做了超過四次迭代，平均誤差值仍然高達2%以上，比計算的預估值還差，由此可知這樣的分析可以大幅降低計算量。

Average time for plain sqrt: 140.21 ns
Average time for magic sqrt: 152.75 ns

然而，看完使用的時間後我發現，雖然預估初始值可以大幅降低牛頓法的迭代次數，但是花費的時間卻比較長，就算 magic sqrt 只做一次牛頓法，速度還是略慢於 plain sqrt 做10次，且可以達到一樣的精確度，就算用了union也是相同的結果，猜測原因是數值在記憶體複製的花費時間可能比想像中的長。

參考資料

牛頓法(無fpu)

和Linux 核心的 PELT的關係，看過原始碼後發現它使用定點數計算EWMA。
所以可以將原本用浮點數迭代的部分改用Q16.16定點數進行運算，並且換算時也要手動提取 IEEE 754 的 float field 進行轉換(參考quiz5)。
輔助函式：

static inline fix16_t float_to_fix16(float a) {
    int32_t i;
    memcpy(&i, &a, sizeof(i));
    int e = ((i >> 23) & 0xFF) - 127;
    int m = (i & 0x7FFFFF) | 0x800000;
    int shift = e - 7;
    if (shift >= 0) {
        return (fix16_t)(m << shift);
    } else {
        return (fix16_t)(m >> -shift);
    }
}

static inline float fix16_to_float(fix16_t a) {
    int pos = 31 - __builtin_clz(a);
    int e = pos - 16;
    int biased = e + 127;
    int m = (a << (23 - pos)) & 0x7FFFFF;
    int32_t ib = (biased << 23) | m;
    float f;
    memcpy(&f, &ib, sizeof(f));
    return f;
}

static inline fix16_t fix16_mul(fix16_t x, fix16_t y)
{
    int64_t res = (int64_t) x * y;
    return (fix16_t) (res >> 16);
}

static inline fix16_t fix16_div(fix16_t a, fix16_t b)
{
    if (b == 0) /* Avoid division by zero */
        return 0;
    fix16_t result = (((int64_t) a) << 16) / ((int64_t) b);
    return result;
}

以下是實作程式碼：

float Q_sqrt_fix16(float number , int iter)
{
    uint32_t i;
    float y;

    y  = number;
    i  = * ( uint32_t * ) &y;                    
    i  = 0x1FC00000 + ( i >> 1 );               
    y  = * ( float * ) &i;

    fix16_t y_fix16 = float_to_fix16(y);
    fix16_t number_fix16 = float_to_fix16(number);

    while(iter-- > 0) { 
        y_fix16 = y_fix16 - fix16_mul(
            float_to_fix16(0.5f),
            y_fix16 - fix16_div(number_fix16, y_fix16)
        );
    }
    y = fix16_to_float(y_fix16);
    
    return y;
}

由於選擇的是 unsigned Q16.16 定點數，所以整數部分可以表示的範圍只有

$[0,2^{15}-1]$ 而小數部分的精確度只能表示到

$2^{-16}\approx 0.00001526$。

誤差分析：

Error percentages for Q16.16 iterations:
Iteration 1: Max Error = 0.001721(0.172147%), Mean Error = 0.000235(0.023524%)
Iteration 2: Max Error = 0.000240(0.024000%), Mean Error = 0.000001(0.000109%)
Iteration 3: Max Error = 0.000240(0.024000%), Mean Error = 0.000001(0.000107%)
Iteration 4: Max Error = 0.000240(0.024000%), Mean Error = 0.000001(0.000107%)

由於前面提過，經過計算的預估值已經非常準確，所以誤差都很小，然而最大誤差卻和Q16.16的精確度差了一個數量級，

分析原因：在 Q16.16 定點格式中，最小可表達增量為

$2^{-16}\approx 1.5\times {10}^{-5}$，若當前迭代誤差

${\mathrm{\Delta }}_k$已接近此量化步長時，下一次更新

${\mathrm{\Delta }}_{k+1}\approx \frac{{\mathrm{\Delta }}_k^2}{2y_k}$將遠小於

$2^{-16}$，因此在定點運算時會被捨入為

$0$，無法再繼續縮小。

令定點最小量化步長為

$q=2^{-16}\approx 1.5259\times {10}^{-5}$，每次牛頓法迭代中，主要有兩次定點除法與一次定點乘法會各自捨入：

$$y_{k+1}=y_k-\underset{\le q}{\underbrace{{\mathrm{round}}_{Q16.16}(\frac{1}{2})}}\times \underset{\le y_k+q}{\underbrace{(y_k-{\mathrm{round}}_{Q16.16}(n/y_k))}}$$

因此，每迭代一次，最多會產生約

$3q$的累積捨入誤差。若迭代

$K$次，則量化誤差上界為

$E_{max}\le K\cdot 3q$，帶入

$K=4,q=1.5259\times {10}^{-5}⟹E_{max}\lesssim 4\times 3\times 1.5259\times {10}^{-5}\approx 1.83\times {10}^{-4}.$
再加上最後一次從定點轉回浮點時，尾數欄位也會捨入一個

$q$,\

$\therefore E_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}\lesssim 1.83\times {10}^{-4}+1.5259\times {10}^{-5}\approx 1.98\times {10}^{-4}$
與實測的

$2.4\times {10}^{-4}$誤差量級相符（因精確誤差會累積，故略大於理論上界）。

Average time for lib_sqrtf: 28.70 ns
Average time for Q_fix16 iter=1: 65.38 ns
Average time for Q_fix16 iter=2: 83.83 ns
Average time for Q_fix16 iter=3: 101.31 ns
Average time for Q_fix16 iter=4: 118.70 ns