2020q3 Homework5 (quiz5)

contributed by < quantabase13 >

程式運作原理

測驗一：

程式碼運作原理如下：

$$\begin{array}{rl}{\displaystyle \frac{a}{b}}-{\displaystyle \frac{a}{b+1}}& ={\displaystyle \frac{a}{b\cdot (b+1)}}\\ & ={\displaystyle \frac{{\displaystyle \frac{a}{b+1}}}{b}}\end{array}$$
換句話說，要計算

${\displaystyle \frac{a}{b}}$ 時，可以換成計算

${\displaystyle \frac{a}{b+1}}+{\displaystyle \frac{{\displaystyle \frac{a}{b+1}}}{b}}$ 。

• 需要討論的事項：

1. 比較該除法與一般除法的效能（實驗）

測驗二：

為了集中討論程式碼的核心，我將題目的程式碼縮減成如下：

#include <stdint.h>

typedef union {
    float value;
    uint32_t v_int;
} ieee_float_shape_type;

/* Set a float from a 32 bit int. */
#define  INSERT_WORDS(d, ix0)	        \
    do {                                \
        ieee_float_shape_type iw_u;     \
            iw_u.v_int = (ix0);         \
        (d) = iw_u.value;               \
    } while (0)

/* Get a 32 bit int from a float. */
#define EXTRACT_WORDS(ix0, d)        \
    do {                             \
        ieee_float_shape_type ew_u;  \
        ew_u.value = (d);            \
        (ix0) = ew_u.v_int;          \
    } while (0)

    float ieee754_sqrt(float x)
    {
        float z;
        uint32_t r;
        int32_t ix0, s0,q,m,t,i;

        EXTRACT_WORDS(ix0, x);
        m = ix0 >> 23;
        m -= 127;
        ix0 = (ix0 & 0x007fffff) | 0x00800000;
        if (m&1) {
            ix0 += ix0;
        }
        m >>= 1;
        q = s0  = 0;
        r = 0x01000000;
        ix0 += ix0;
        while (r != 0) {
            t = s0 + r;
            if (t <= ix0) {
                s0 = t + r;
                ix0 -= t;
                q += r;
            }
            ix0 += ix0;
            r >>= 1;
        }
        ix0 = (q >> 1) + 0x3f000000;
        ix0 += (m << 23);

        INSERT_WORDS(z, ix0);
        return z;

    }

根據 IEEE-754 ，浮點數格式轉換成科學記號的關係如下：

$$(-1{)}^{sign}\cdot (1.fraction{)}_2\cdot 2^{exponent-bias}$$
在 sign bit 為 0 (這裡我們只考慮正數)時，平方根可以用

$$\sqrt{(1.fraction{)}_2}\cdot 2^{{\displaystyle \frac{exponent-bias}{2}}}$$
表示。於是，我們可以將計算拆成 fraction 和 exponent 兩部份。
以下是取得 fraction 和 exponent 的步驟：

• 第一步：取得指數

        m = ix0 >> 23;
        m -= 127;

這裡是根據 IEEE-754 的規範，將指數的部份取出，並減去 bias = 127 。
另外，指數有可能為奇數，之後除以 2 時要把多出的 2 乘回 1.fraction，這個想法對應到程式碼

        if (m&1) {
            ix0 += ix0;
        }

• 第二步：取得 fraction

ix0 = (ix0 & 0x007fffff) | 0x00800000;

這裡取得 fraction 的部份只有 (ix0 & 0x007fffff) 這裡，真正的目的是令 ix0 值為

$(1.\stackrel{23\mathrm{個}bit}{\overbrace{fraction}}{)}_2$ ，binary point 在從 LSB 數來第 23 個 bit 和第 24 個 bit 之間。
到此取得了 fraction 、 exponent。平方根的指數部份能簡單求出，重點是

$\sqrt{(1.fraction{)}_2}$ 的求法。參考 YLowy 的筆記 後，我發現我需要了解開 2 次方根的運算原理。
我查到了開 n 次方根的直式計算與原理 這篇。
以求 根號2 為例：

1. 根號2
   2 的
   $(1.fraction{)}_2$ 左移兩位。總共26 bit。此時 binary point 定在從 LSB 數來第 24 和第 25 個 bit 之間。

10|00|00|00|00|00|00|00|00|00|00|00|00

為了討論方便，在做開平方根的直式計算時，暫且將

10|00|00|00|00|00|00|00|00|00|00|00|00

視為整數來計算。
以下是計算流程：

  q                        1|  0| 1|
                        --------------------------------------------------
 ix0              √ 10|00|00|00|00|00|00|00|00|00|00|00|00
t = 1                    1
s0 = 10           
                        --------------------------------------------------
t = 101                100
s0 = 100              000
                        -------------------------------------------------
 t = 1001              10000
 s0 = 1010              1001

q 欄位代表平方根 ， ix0 欄位代表代表要被取平方根的數，t 欄位代表將要從 ix0 扣除的數， s0 欄位代表 q 左移一位的結果。
這種直式計算，用到的原理其實就是國中學的乘法公式：
若

$$x={(a+\beta )}^2=a^2+(2a+\beta )\beta$$
則

$$x-a^2=(2a+\beta )\beta$$
若

$$\beta =b+\gamma$$
則

$$x-(a+b{)}^2=(2a+2b+\gamma )\gamma$$
試著把上述兩個式子套到直式計算的流程：

$$a^2=(1\stackrel{13\mathrm{個}bit}{\overbrace{(00...00)}}{)}^2<=ix0<(10\stackrel{13\mathrm{個}bit}{\overbrace{(00...00)}}{)}^2$$ 所以

$$ix0=(1\stackrel{13\mathrm{個}bit}{\overbrace{(\beta )}}{)}^2=(a+\beta {)}^2$$

$$(a+b{)}^2=(1\stackrel{13\mathrm{個}bit}{\overbrace{(00...00)}}{)}^2<=ix0<(1\stackrel{13\mathrm{個}bit}{\overbrace{(10...00)}}{)}^2$$

$a^2+2ab+b^2=a^2+(2a+b)b$
所以

$$ix0=(1\stackrel{13\mathrm{個}bit}{\overbrace{(\beta )}}{)}^2=(a+b+\gamma {)}^2=(10\stackrel{12\mathrm{個}bit}{\overbrace{(\gamma )}}{)}^2$$
由於

$$a^2+2ab+b^2=a^2+(2a+b)b$$
比較

$(a+b{)}^2$是否大於 ix0 時，只要比較

$(2a+b)b$ 是否大於

$ix0-a^2$即可。
我們求平方根的求法是由左往右 bit by bit 的求，因此 b 永遠是當前 q 的 LSB，非 0 即 1。換句話說，只須比較

$(2a+1)$ 是否大於

$ix0-a^2$，就可以知道 b 是否為 1 。這就是

        while (r != 0) {
            t = s0 + r;
            if (t <= ix0) {
                s0 = t + r;
                ix0 -= t;
                q += r;
            }
            ix0 += ix0;
            r >>= 1;
        }

這段程式碼作的工作。另外的 ix0 += ix0 是為了讓 s0 可以不用 shift。
順帶一提，如果要更符合數學式的話，可以寫成

        while (r != 0) {
            t = s0 + r;
            if (t <= ix0) {
                ix0 -= t;
                q += r;
                s0 = q << 1;
            }
            ix0 += ix0;
            r >>= 1;
        }

測驗三：

從 x 開始， k 個連續正整數合為 N 的話，可以表示成如下關係式：

$kx=N-{\displaystyle \frac{(k-1)k}{2}}$
於是我們要檢查特定 k 是否能使上述關係式成立時，我們只須確認

$N-{\displaystyle \frac{(k-1)k}{2}}$ 是否為 k 的倍數即可。
接下來關於 k 的範圍，我們知道

$N-{\displaystyle \frac{(k-1)k}{2}}>0$ ，故將其設為迴圈中止條件即可。