2020q3 Homework5 (quiz5)

contributed by < Holycung >
2020q3 Homework5 (quiz5) 題目

測驗 1

題目

考慮到以下浮點數除法程式: (fdiv.c)












#include <stdio.h>
#include <stdlib.h>

double divop(double orig, int slots) {
    if (slots == 1 || orig == 0)
        return orig;
    int od = slots & 1;
    double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);
    if (od)
        result += divop(result, slots);
    return result;
}

假設 divop() 的第二個參數必為大於 0 的整數，而且不超過 int 型態能表達的數值上界。請補完程式碼。

作答

回答這題要先了解一下在 IEEE 754 下浮點數的定義。

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

這題的觀念是把原本的浮點數除法都替換成除以 2 進行約分，因為整數除以二可以用往右 shift 一位來加速，浮點數除以二只需要把 exponent 的部分減一就好。那就可以把狀況分成下面兩種，假設

X \div Y

，已知 Y 為大於零的整數，第一種狀況就是

Y

為偶數，那我們可以直接對

X Y

都向右 shift 一位進行約分，也就是

X \div Y = \frac{X}{2} \div \frac{Y}{2}

。第二種狀況就是

Y

為奇數，這就是我們所要探討的。

Y

為奇數了話我們就沒有辦法直接除二，那這邊的做法就是先將

Y + 1

，把原本的

X \div Y

變成

X \div (Y + 1)

，這樣

Y + 1

就是偶數就可以用第一種狀況的方法，不過我們知道

X \div (Y + 1)

跟原本的值會不同，且

X \div Y > X \div (Y + 1)

，所以我們會需要補差值回去，那個差值就是

X \div Y - X \div (Y + 1)

，稍微化減一下如下。

\begin{aligned} 差 值 & = \frac{X}{Y} - \frac{X}{Y + 1} \\ = \frac{X (Y + 1)}{Y (Y + 1)} - \frac{X Y}{Y (Y + 1)} \\ = \frac{X}{Y (Y + 1)} \end{aligned}

所以這邊我們得到了差值

= \frac{X}{Y (Y + 1)}

，可以發現這邊

\frac{X}{(Y + 1)}

跟我們前面第二種狀況是一樣的，接下來我們就可以來看程式。












#include <stdio.h>
#include <stdlib.h>

double divop(double orig, int slots) {
    if (slots == 1 || orig == 0)
        return orig;
    int od = slots & 1;
    double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);
    if (od)
        result += divop(result, slots);
    return result;
}

第 5 行這邊是判斷如果，除數 slot 是 1 或是被除數 orig 是 0 就直接回傳被除數 orig。
第 7 行 od 是判斷除數是否為奇數，只要拿最後一個位元跟 1 做 bitwise and 就可以得到。
第 8 行這邊就是我們剛剛上面推論的方法，把浮點數的除法全部換成除以二約分遞迴下去，divop 第一個參數是 orig / D1，就是對被除數除以二的動作，所以 D1 = 2這邊是浮點數除二所以能用 shift。
第二個參數會先判斷除數 slots 是否為奇數，如果是偶數就直接進行 shift，如果是奇數就加一後在進行 shift，所以 D2 = 1。最後把遞迴的結果存回到 result 變數，所以 result =

\frac{X}{(Y + 1)}

。
第 9 10 行最後就是要判斷除數是否為奇數，如果是奇數了話我們需要把差值補回去，經過剛剛的推導差值

= \frac{X}{Y (Y + 1)}

，也就是

\frac{r e s u l t}{Y}

。

延伸問題

TODO

以編譯器最佳化的角度，推測上述程式碼是否可透過 tail call optimization 進行最佳化，搭配對應的實驗;
比照浮點數運算和定點數操作，改寫上述程式碼，提出效率更高的實作，同樣避免使用到 FPU 指令 (可用 gcc -S 觀察對應的組合語言輸出)，並與使用到 FPU 的實作進行效能比較

測驗 2

題目

延續從

\sqrt 2

的運算談浮點數，假設 float 為 32-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式，請補上對應數值，使其得以運作：
































































































#include <stdint.h>

/* A union allowing us to convert between a float and a 32-bit integers.*/
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)

static const float one = 1.0, tiny = 1.0e-30;

float ieee754_sqrt(float x)
{
    float z;
    int32_t sign = 0x80000000;
    uint32_t r;
    int32_t ix0, s0, q, m, t, i;
	
    EXTRACT_WORDS(ix0, x);
	
    /* take care of zero */
    if (ix0 <= 0) {
        if ((ix0 & (~sign)) == 0)
            return x; /* sqrt(+-0) = +-0 */
        if (ix0 < 0)
            return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
    }
    /* take care of +INF and NaN */
    if ((ix0 & 0x7f800000) == 0x7f800000) {
        /* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
        return x;
    }
    /* normalize x */
    m = (ix0 >> 23);
    if (m == 0) { /* subnormal x */
        for (i = 0; (ix0 & 0x00800000) == 0; i++)
            ix0 <<= 1;
        m -= i - 1;
    }
    m -= 127; /* unbias exponent */
    ix0 = (ix0 & 0x007fffff) | 0x00800000;
    if (m & 1) { /* odd m, double x to make it even */
        ix0 += ix0;
    }
    m >>= 1; /* m = [m/2] */
	
    /* generate sqrt(x) bit by bit */
    ix0 += ix0;
    q = s0 = 0; /* [q] = sqrt(x) */
    r = QQ0;       /* r = moving bit from right to left */

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

    /* use floating add to find out rounding direction */
    if (ix0 != 0) {
        z = one - tiny; /* trigger inexact flag */
        if (z >= one) {
            z = one + tiny;
            if (z > one)
                q += 2;
            else
                q += (q & 1);
        }
    }
	
    ix0 = (q >> 1) + QQ1;
    ix0 += (m << QQ2);
	
    INSERT_WORDS(z, ix0);

    return z;
}

請補完程式碼，使上述程式碼得以運作。

作答

再回答這題前可以先看過浮點數運算和定點數操作，熟悉一下 IEEE 754 對浮點數的規範，再來直接看到程式碼。





/* A union allowing us to convert between a float and a 32-bit integers.*/
typedef union {
    float value;
    uint32_t v_int;
} ieee_float_shape_type;

在 C 當中的 union 用法可以參考 Union type wiki。

a union is a value that may have any of several representations or formats within the same position in memory;















/* 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)

這邊有用到你所不知道的C語言：技巧篇提到的 do { ... } while(0) 巨集 是為了避免 dangling else 問題，也可以參考這篇 macro do{}while(0) 更詳細的解釋。
這兩個 macro INSERT_WORDS、EXTRACT_WORDS，是運用 union 來轉換 float 跟 int，參考 RinHizakura 的例子如下:

SIGN  EXPONENT  FRACTION
0     00000010  01000000000000000000000

如果用 uint32_t 的角度來解釋，就是

2^{24} + 2^{21}

如果從 float (IEEE 754) 的角度來解釋，則是

(1 + 1 * 2^{- 2}) \times 2^{(2 - 127)}









    EXTRACT_WORDS(ix0, x);
	
    /* take care of zero */
    if (ix0 <= 0) {
        if ((ix0 & (~sign)) == 0)
            return x; /* sqrt(+-0) = +-0 */
        if (ix0 < 0)
            return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
    }

第 34 行先把浮點數 x 轉成整數 ix0。
第 37 行這邊是先處理小魚等於 0 的數字，第 38 行是判斷是否為正零或負零(除了 sign bit 其他都是 0)，如果是 0 的話就回傳 0。
第 40 行是判斷是否小於 0，對小於零的數取平方根的話要回傳 NaN，第 41 行是透過 0.0 / 0.0 產生 NaN，這邊可以從 IEEE 754-2019 規格書當中 6.2.1 先了解到 NaN 的編碼格式，然後又分程 quiet NaN 和 signal NaN，再看到 7.2 有明確規範當浮點數 0.0 / 0.0 的時候會產生 quiet NaN。

6.2.1 NaN encodings in binary interchange formats
All binary NaN bit strings have the sign bit S set to 0 or 1 and all the bits of the biased exponent field E set to 1 (see 3.4). A quiet NaN bit string should be encoded with the first bit (d1) of the trailing significand field T being 1. A signaling NaN bit string should be encoded with the first bit of the trailing significand field being 0.

7.2 Invalid operation
For operations producing results in floating-point format, the default result of an operation that signals the invalid operation exception shall be a quiet NaN that should provide some diagnostic information.
e) division: division(0, 0) or division(∞, ∞)





    /* take care of +INF and NaN */
    if ((ix0 & 0x7f800000) == 0x7f800000) {
        /* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
        return x;
    }

第 44 行是判斷 +INF 和 NaN 如果是，ix0 & 0x7f800000 == 0x7f800000 就可以判斷 Exponent 是否等於 255，也就是 INF、NaN 的狀況。













    /* normalize x */
    m = (ix0 >> 23);
    if (m == 0) { /* subnormal x */
        for (i = 0; (ix0 & 0x00800000) == 0; i++)
            ix0 <<= 1;
        m -= i - 1;
    }
    m -= 127; /* unbias exponent */
    ix0 = (ix0 & 0x007fffff) | 0x00800000;
    if (m & 1) { /* odd m, double x to make it even */
        ix0 += ix0;
    }
    m >>= 1; /* m = [m/2] */

第 48 行把 ix0 向右 shift 23 個得到 m 就是 exponent 的值，如果 m == 0 代表是 denormalized number，這邊就需要先把他變成 normalized number 的形式一起處理，所以這邊的想法是把 ix0 往左 shift 也就是把 fraction 的部分往左移動，直到第 24 個位元是 1，然後記錄下移動了幾次 i，每移動一次就相當指數 m 要減一，所以最後要把 m -= i - 1，就變成了 normalized 的形式。

因為剛剛處理完 denormalized，這邊就可以直接用 normalized 的形式處理，也就是

1. F \times 2^{E - 127}

，所以第 55 行就把 m 減掉 127。

第 56 行是把 ix0 只留下後面 fraction 23 bits 的部分，把 24 位元變成 1，因為前面已經把 exponent 的部分保留下來到 m，然後也把 denormalized 的部分變成 normalized 的形式，所以第 24 位元本來就是 1，因此這邊就只留下後面 24 位元也就是 1.F。

第 57 行是判斷 m 是否為偶數，因為對

2^{E}

取平方根了話就會是

2^{E / 2}

，那如果 m 不是偶數了話，就把 ix0 乘以二，這樣就可以把 m 除以二。















    /* generate sqrt(x) bit by bit */
    ix0 += ix0;
    q = s0 = 0; /* [q] = sqrt(x) */
    r = QQ0;       /* r = moving bit from right to left */

    while (r != 0) {
        t = s0 + r;      // t = s_i + r_i
        if (t <= ix0) {      
            s0 = t + r;  // s_{i+1} = t + r_i
            ix0 -= t;    // x_{i+1}/2 = x_i - t
            q += r;      // q_{i+1} = q_i + r_i
        }
        ix0 += ix0;      // x_{i+1}
        r >>= 1;
    }

接著為計算 square root 的部分，這邊的方法可以參考 How to calculate a square root without a calculator，其方法的證明可以參考 Why the square root algorithm works，應用到二進位 Methods of computing square roots 也是一樣的，我們則可以利用這個方法把平方根換成較快的 bit shift operations，接下來進行推導，推導這部分是參考 RusselCK。

這邊我們要計算
$x$ 的平方根
$\sqrt{x}$ bit by bit 的數學推導。
- 令
  $q_{i} = \sqrt{x}, q_{i}$ 代表到小數點後面
  $i$ 位精確值
  $i \geq 0$
- 那假設今天
  $1 \leq x < 4$ ，那
  $\sqrt{x}$ 的整數位必定是 1，所以
  $q_{0} = 1$ ，就會後面會探討為甚麼要這樣假設。
- 考慮
  $q_{i + 1}$ ，也就是 bit by bit 的逼近下去:
  1. 若
    $(q_{i} + 2^{- (i + 1)})^{2} \leq x$ ，則
    $q_{i + 1} = q_{i} + 2^{- (i + 1)}$
    - 即
      $q_{i}$ 下一位補 1
  2. 若
    $(q_{i} + 2^{- (i + 1)})^{2} > x$ ，則
    $q_{i + 1} = q_{i}$
    - 即
      $q_{i}$ 下一位補 0
- 整理上面第一種情況
  $(q_{i} + 2^{- (i + 1)})^{2} \leq x$ ：
  
  $\begin{aligned} {(q_{i} + 2^{- (i + 1)})}^{2} \leq x & \Rightarrow q_{i}^{2} + 2 \times q_{i} \times (2^{- (i + 1)}) + (2^{- (i + 1)})^{2} \leq x \\ \Rightarrow q_{i} \times (2^{- i}) + 2^{- 2 (i + 1)} \leq x - q_{i}^{2} \\ \Rightarrow q_{i} \times 2 + 2^{- (i + 1)} \leq (x - q_{i}^{2}) \times 2^{(i + 1)} \\ \Rightarrow s_{i} + r_{i} \leq x_{i}, \\ \begin{aligned} w h e r e & s_{i} = 2 \times q_{i}, \\ r_{i} = 2^{- (i + 1)}, \\ x_{i} = (x - q_{i}^{2}) \times 2^{(i + 1)} = (x - q_{i}^{2}) \times r_{i}^{- 1} \end{aligned} \end{aligned}$
- 因此
  - 若
    $s_{i} + r_{i} \leq x_{i}$ ，則
    - $q_{i + 1} = q_{i} + r_{i}$
    - $s_{i + 1} = (s_{i} + r_{i}) + r_{i}$
    - $x_{i + 1} = 2 x_{i} - 2 (s_{i} + r_{i}) = 2 [x_{i} - (s_{i} + r_{i})]$
    推導過程
    
    $\begin{aligned} s_{i + 1} & = 2 \times q_{i + 1} \\ = 2 \times [q_{i} + 2^{- (i + 1)}] \\ = (2 \times q_{i}) + 2^{- i} \\ = s_{i} + 2^{- i} \\ = s_{i} + 2^{- (i + 1)} + 2^{- (i + 1)} \\ = (s_{i} + r_{i}) + r_{i} \end{aligned}$
    
    $\begin{aligned} x_{i + 1} & = [x - q_{i + 1}^{2}] \times 2^{((i + 1) + 1)} \\ = [x - (q_{i} + r_{i})^{2}] \times 2^{(i + 2)} \\ = [x - (q_{i}^{2} + 2 q_{i} r_{i} + r_{i}^{2})] \times 2 r_{i}^{- 1} \\ = [x - q_{i}^{2}] \times 2 r_{i}^{- 1} - 2 q_{i} r_{i} \times 2 r_{i}^{- 1} - r_{i}^{2} \times 2 r_{i}^{- 1} \\ = 2 x_{i} - 2 s_{i} - 2 r_{i} \\ = 2 x_{i} - 2 (s_{i} + r_{i}) \end{aligned}$
  - 若
    $s_{i} + r_{i} > x_{i}$ ，則
    - $q_{i + 1} = q_{i}$
    - $s_{i + 1} = s_{i}$
    - $x_{i + 1} = 2 x_{i}$
    推導過程
    
    $\begin{aligned} s_{i + 1} & = 2 \times q_{i + 1} \\ = 2 \times q_{i} \\ = s_{i} \end{aligned}$
    
    $\begin{aligned} x_{i + 1} & = [x - q_{i + 1}^{2}] \times 2^{((i + 1) + 1)} \\ = [x - q_{i}^{2}] \times 2^{(i + 2)} \\ = [x - q_{i}^{2}] \times 2 r_{i}^{- 1} \\ = 2 x_{i} \end{aligned}$

接下來我們可以回到程式碼。















    /* generate sqrt(x) bit by bit */
    ix0 += ix0;
    q = s0 = 0; /* [q] = sqrt(x) */
    r = QQ0;       /* r = moving bit from right to left */

    while (r != 0) {
        t = s0 + r;      // t = s_i + r_i
        if (t <= ix0) {      
            s0 = t + r;  // s_{i+1} = t + r_i
            ix0 -= t;    // x_{i+1}/2 = x_i - t
            q += r;      // q_{i+1} = q_i + r_i
        }
        ix0 += ix0;      // x_{i+1}
        r >>= 1;
    }

第 63 行，把 ix0 往左 shift 一位。這樣 ix0 的寬度可能是 25 或 26 個位元，所以 QQ0 = 0x01000000，第 25 個位元為一，等於把小數點點在第 24 25 個位元之間，從第 25 個 bit 開始，這樣代表我們會精確到小數點後面第 24 位，那在 IEEE 754 的規範中 fraction 部分有 23 位，我們就可以把最後一位做 rounding，精確到小數點後 23 位。小數點點在這邊也可以確保我們前面的假設

1 \leq x < 4

是正確的。


static const float one = 1.0, tiny = 1.0e-30;











    /* use floating add to find out rounding direction */
    if (ix0 != 0) {
        z = one - tiny; /* trigger inexact flag */
        if (z >= one) {
            z = one + tiny;
            if (z > one)
                q += 2;
            else
                q += (q & 1);
        }
    }

第 79 行這邊是判斷 rounding，首先要先判斷系統的 rounding mode，這邊透過 1 加減一個很小的數 tiny，讓系統進行 rounding 來判斷其 rounding mode。

先用 one - tiny 如果小於 1，代表是 rounding down，那不管最低位是 0 或 1 都會被捨棄掉，所以不需要額外動作。
接下來判斷 one + tiny 如果大於 1，代表是 rounding up，那 23 bit 就直接直接進位，所以這邊是 q += 2。
最後如果 one + tiny 如果等於 1，代表是 round to nearest，此時就要判斷最低位是 0 或 1，如果是 1 就要進位，所以是 q += (q & 1)。






    ix0 = (q >> 1) + QQ1;
    ix0 += (m << QQ2);
	
    INSERT_WORDS(z, ix0);

    return z;

最後要把 exponent 跟 fraction 的部分合併回來，所以第 90 行先把 q 往右 shift 1 位，不過要注意 q 現在還剩下 24 位，因為剛剛整理得到的是

1. F

所以我們這邊的 QQ1 = 0x3f000000，讓中間 exponent 最低 7 位 bit 都補上一，因為 IEEE 754 normalized 的格式是

1. F \times 2^{E - 127}

，也就是在這邊 exponent 先加上 127。第 91 行，在把 m 往左 shift 23 個，加回去 ix0。

最後在透過 INSERT_WORDS 轉回浮點數回傳。

延伸問題

練習 LeetCode 69. Sqrt(x)，提交不需要 FPU 指令的實作。

作答

從

\sqrt 2

的運算談浮點數了解到開根號可以透過二分逼近法或是十分逼近法來找到，再研讀以牛頓法求整數開二次方根的近似值 (第 19 頁)，決定使用牛頓法來實作。

\sqrt[n]{N} ≒ [(n - 1) \times a + \frac{N}{a^{n - 1}}] \div n = a_{1}

透過多次逼近可以取得更精準的值，今天我們要找到平方根，所以 n=2:

\sqrt{N} ≒ [(2 - 1) \times a + \frac{N}{a^{2 - 1}}] \div 2 = \frac{a + \frac{N}{a}}{2}

接下來看到程式。














int mySqrt(int x){
    if (x == 0 || x == 1)
        return x;
    if (x == INT_MAX) // avoid overflow
        x-=1;
    
    unsigned int a; // using right shift
    a = (x + 1) >> 1;
    
    while (a > x/a){
        a = (a + x/a) >> 1;
    }
    return a;
}

那最開始的值要從多少開始逼近呢，這邊可以從 INT_MAX 最大值開始逼近，但是我們可以透過算幾不等式，找到比

\sqrt{x}

大但又不會太大的值開始，轉換成程式碼就是 (x + 1) >> 1，不過這邊要小心因為會對 x 做加一的動作，所以如果 x == INT_MAX 要避免 overflow 就要先把他減一。

\begin{aligned} \frac{a + b}{2} \geq & \sqrt{a b}, w h e r e a, b > 0 \\ => \frac{a + 1}{2} \geq & \sqrt{a}, w h e r e a > 0, b = 1 \end{aligned}

最後 while 迴圈的終止條件則是找到第一個

a^{2} \leq x

的

a

。

另外還有找到一個快速整數平方根算法，有時間了話再來研究QQ












int mySqrt(int v)
{
    unsigned long temp, nHat = 0, b = 0x8000, bshft = 15;
    do{
        if (v >= (temp = (((nHat << 1) + b) << bshft--)))
        {
            nHat += b;
            v -= temp;
        }
    } while (b >>= 1);
    return nHat;
}

測驗 3

題目

LeetCode 829. Consecutive Numbers Sum 給定一個正整數 N，問 N 能寫成多少種連續正整數之和，例如 9 可寫成 4+5，或者　2+3+4。由於要寫成連續正整數之和，則可用等差數列來思考，且差值為 1，這個等差數列不必從 1 開始。
參考實作如下:













int consecutiveNumbersSum(int N)
{
    if (N < 1)
        return 0;
    int ret = 1;
    int x = N;
    for (int i = 2; i < x; i++) {
        x += ZZZ;
        if (x % i == 0)
            ret++;
    }
    return ret;
}

請補完程式碼，使上述程式碼得以運作。

作答

假設從 x 開始，且個數共有 k 個，則可寫出該等差數列為：

x, x + 1, x + 2 +, . . ., x + (k - 1)

經過化減後：

\begin{aligned} k x & + \frac{(k - 1) k}{2} = N \\ k x & = N - \frac{(k - 1) k}{2} \\ k x & = N - [1 + 2 + . . . + (k - 1)] \end{aligned}

因此我們就可以透過

{N - [1 + 2 + . . . + (k - 1)]} m o d k = 0

來判斷這個 k 長度的連續整數是否是一組解。

第 5 行先初始化回傳值為 1，因為任何大於 0 的整數都必定會有自己的一個解。
第 7 行的迴圈就是從 k=2 的長度開始找起，第 8 行就是

{N - [1 + 2 + . . . + (k - 1)]}

，所以 ZZZ = x += (1 - k)，第 9 行就是判斷能不能整除，可以的話就代表可以找到一個解 ret++。

延伸問題

TODO

嘗試改寫為更有效率的實作

2020q3 Homework5 (quiz5)

目錄

測驗 1

題目

作答

延伸問題

測驗 2

題目

作答

延伸問題

作答

測驗 3

題目

作答

延伸問題

參考資料

2020q3 Homework5 (quiz5)

目錄

測驗 1

題目

作答

延伸問題

測驗 2

題目

作答

延伸問題

作答

測驗 3

題目

作答

延伸問題

參考資料

Read more

2020q3 Homework5 (render)

研讀筆記: Raycasting

研讀筆記: Achieving polymorphism in C

研讀筆記 linenoise