2018q1 Homework2 (assessment)

contributed by < TingL7 >

第 1 週 測驗 一

考慮以下程式碼:

#include <stdlib.h>
int isMultN(unsigned int n) {
    int odd_c = 0, even_c = 0; /* variables to count odd and even SET bits */
    
    if (n == 0)    // return true if difference is 0.
        return 1;
    if (n == 1)    // return false if the difference is not 0.
        return 0;
    
    while (n) {
        if (n & 1)   // odd bit is SET, increment odd_C
            odd_c++;
        n >>= 1; 
    
        if (n & 1)   // even bit is SET, increment even_c
            even_c++;
        n = n >> 1;
    }
    
    /* Recursive call till you get 0/1 */
    return(isMultN(abs(odd_c - even_c)));
}

其作用為檢查輸入整數是否為 N 的倍數，那麼 N 為多少？

---

解題想法&思考

• 程式碼解讀
  • 由這兩段程式碼，可以得知這是一個遞迴程式，結束條件是 odd_c 和 even_c 相等或相差 1 ，且相等時，為N的倍數。
    ​​​​​​​​if (n == 0)    // return true if difference is 0.
    ​​​​​​​​    return 1;
    ​​​​​​​​if (n == 1)    // return false if the difference is not 0.
    ​​​​​​​​    return 0;
    

    return(isMultN(abs(odd_c - even_c)));
  • 而由這段程式碼可以知道 odd_c 及 even_c 分別在計算奇數位及偶數位為 1 的個數
    ​​​​​​​​while (n) {
    ​​​​​​​​    if (n & 1)   // odd bit is SET, increment odd_C
    ​​​​​​​​        odd_c++;
    ​​​​​​​​    n >>= 1; 
    ​​​​
    ​​​​​​​​    if (n & 1)   // even bit is SET, increment even_c
    ​​​​​​​​        even_c++;
    ​​​​​​​​    n = n >> 1;
    ​​​​​​​​}
    

  • 綜合上述，並簡單化這個題目，在程式只有進行一次遞迴時，奇數位及偶數位總合相等時，即為 N 的倍數。將選項帶入有此特性的只有 3(11₂) 。

延伸問題: 為什麼上述程式碼可運作呢？

數字與文字間也要空白喔！
課程助教

好的，已修改
TingL7

• 程式碼可以運作，是依據 3 的倍數，在二進位的情況下的特性，即奇數位數字和及偶數位數和相等或相差 3 的倍數。因此，只要證明 3 的倍數有此特性，即可確認上述程式碼是運作良好的。

• 這個特性與十進位中的 11 很像。參考 100 以內的質數倍數判別法及其延伸中， 11 的倍數判別法，可寫出以下證明:

  設有一 (n+1) 位數 N = (a_na_n-1..a₁a₀)₂ ，其中 a_n

  $\ne$ 0
  令 k 為
  $\le n$ 的正整數
  
  $$\begin{array}{rl}{a}_n& =\sum _{i=0}^n{a}_i\times 2^i\\ & =[a_0+a_2(3+1)+a_4(15+1)+...+a_{2k}(2^{2k}-1+1)+...]\\ & +[a_1(3-1)+a_3(9-1)+a_5(33-1)+...+a_{2k+1}(2^{2k+1}+1-1)+...]\\ & =(a_0+a_2+a_4+...+a_{2k}+...)\\ & +[a_2(3)+a_4(15)+...+a_{2k}(2^{2k}+1)+...]\\ & -(a_1+a_3+a_5+...+a_{2k+1})\\ & +[a_1(3)+a_3(9)+a_5(33)+...+a_{2k+1}(2^{2k+1}+1)+...]...\mathrm{第}\mathrm{一}\mathrm{式}\end{array}$$
  利用歸納證明法，可證明 2²ⁿ-1 及 2²ⁿ⁺¹+1 為 3 的倍數:

  n=1:則 2²-1=3 為 3 的倍數
  令 n=k 成立，即 2^2k-1=3x ， x 可為任一正整數
  n=k+1:
  2^2(k+1)-1=2^2k+2-1=4*2^2k-1=4(2^2k-1)+3=4(3x)+3=3(4x+1)
  得證: 2²ⁿ-1 為 3 的倍數
  因此， [a₂(3)+a₄(15)+…+a_2k(2^2k+1)+…] 為 3 的倍數…第二式

  n=0:則 2¹-1=3 為 3 的倍數
  令 n=k 成立，即 2^2k+1+1=3x ， x 可為任一正整數
  n=k+1:
  2^2(k+1)+1+1=2^2k+3+1=4*2^2k+1+1=4(2^2k+1+1)-3=4(3x)-3=3(4x-1)
  得證: 2²ⁿ⁺¹-1 為 3 的倍數
  因此， [a₁(3)+a₃(9)+a₅(33)+…+a_2k+1(2^2k+1+1)+…] 為 3 的倍數…第三式

  由第一、二、三式可得知:
  如果 (a₀+a₂+a₄+…+a_2k+…)-(a₁+a₃+a₅+…+a_2k+1+…) 為 3 的倍數，則 N 必為 3 的倍數。
  因此，要判定二進位中，是否為 3 的倍數，可由 奇數位數字和及偶數位數和相等或相差 3 的倍數 的特性去判斷。

延伸問題: 將 N 改為其他質數再改寫為對應的程式碼，一樣使用遞迴

• 要將 N 改為其他質數，要先找出其他質數是否類似 3 的特殊判斷特性。
• N=5
  • 參考從二進位判斷數字是否被 5 整除， N 是否為 5 的倍數判定，可由 (N>>2)-(N&3) 是否為 5 的倍數決定。
  • 證明:
    (N>>2)-(N&3) 為 N 除以 4 的商數 (a) 減去 N 除以 4 的餘數 (b)。寫成算式為
    
    $$\begin{gathered} N÷4=a...b \\ =>N=a\times 4+b=a\times 5-a+b \\ =>\mathrm{如}\mathrm{果}b-a\mathrm{為}5\mathrm{的}\mathrm{倍}\mathrm{數}，\mathrm{則}N\mathrm{為}5\mathrm{的}\mathrm{倍}\mathrm{數} \end{gathered}$$
  • 程式碼
  ​​​​#include <stdlib.h>
  ​​​​int isMultN(unsigned int n) {
  ​​​​    if (n == 0|| n == 5)    // return true if n is 0 or 5.
  ​​​​        return 1;
  ​​​​    if (n < 5)    // return false n is smaller than 5
  ​​​​        return 0;
  ​​​​
  ​    n = (n >> 2) + (n & 3); 
  ​​​​
  ​​​​    /* Recursive call till n <= 5 */
  ​​​​    return(isMultN(n));
  ​​​​}
  

• N=7
  • 由 N=5 的證明，是推廣到 N=7 ，甚至可以說是任何數字都可以用，只要找出與 2 的次方的關係。
  • $$\begin{gathered} N=a\times 8+b=a\times 7+(a+b) \\ =>\mathrm{如}\mathrm{果}a+b\mathrm{為}7\mathrm{的}\mathrm{倍}\mathrm{數}，\mathrm{則}N\mathrm{為}7\mathrm{的}\mathrm{倍}\mathrm{數} \end{gathered}$$
  • 程式碼
  ​​​​#include <stdlib.h>
  ​​​​int isMultN(unsigned int n) {
  ​​​​    if (n == 0|| n == 7)    // return true if n is 0 or 7.
  ​​​​        return 1;
  ​​​​    if (n < 7)    // return false n is smaller than 7.
  ​​​​        return 0;
  ​​​​
  ​    n = (n >> 3) + (n & 7); 
  ​​​​
  ​​​​    /* Recursive call till n <= 5 */
  ​​​​    return(isMultN(n));
  ​​​​}
  

• N=3
  • 由上述推得的方法，可以將 N=3 改寫。
  • 改寫一:
    
    $$\begin{gathered} N=a\times 4+b=a\times 3+(a+b) \\ =>\mathrm{如}\mathrm{果}a+b\mathrm{為}3\mathrm{的}\mathrm{倍}\mathrm{數}，\mathrm{則}N\mathrm{為}3\mathrm{的}\mathrm{倍}\mathrm{數} \end{gathered}$$
  • 程式碼
  ​​​​#include <stdlib.h>
  ​​​​int isMultN(unsigned int n) {
  ​​​​    if (n == 0|| n == 3)    // return true if n is 0 or 3.
  ​​​​        return 1;
  ​​​​    if (n < 3)    // return false n is smaller than 3.
  ​​​​        return 0;
  ​​​​
  ​    n = (n >> 2) + (n & 1); 
  ​​​​
  ​​​​    /* Recursive call till n <= 5 */
  ​​​​    return(isMultN(n));
  ​​​​}
  

  • 改寫二:
    
    $$\begin{gathered} N=a\times 2+b=a\times 3+(-a+b) \\ =>\mathrm{如}\mathrm{果}b-a\mathrm{為}3\mathrm{的}\mathrm{倍}\mathrm{數}，\mathrm{則}N\mathrm{為}3\mathrm{的}\mathrm{倍}\mathrm{數} \end{gathered}$$
  • 程式碼
  ​​​​#include <stdlib.h>
  ​​​​int isMultN(unsigned int n) {
  ​​​​    if (n == 0)    // return true if n is 0.
  ​​​​        return 1;
  ​​​​    if (n == 1)    // return false n is 1.
  ​​​​        return 0;
  ​​​​
  ​    n = (n >> 1) - (n & 1); 
  ​​​​
  ​​​​    /* Recursive call till n <= 5 */
  ​​​​    return(isMultN(n));
  ​​​​}
  

• 通用
  • 其實，上面的方法是可以通用的，不論是否為質數，只要算清楚與 2 的次方差多少，再利用 shift 及 + 或 - 就可以判斷倍數。

第 2 週 測驗 一

請完成下方程式碼，依循 IEEE 754 單精度規範，輸出

$2^x$ 的浮點數表達法，而且考慮兩種狀況：

• 當
  $x$ 過小的時候，回傳
  $0.0$
• 當
  $x$ 過大的時候，回傳
  $+\infty$

注意：這裡假設 u2f 函式返回的浮點數值與其無號數輸入有著相同的位數，也就是至少 32-bit

#include <math.h>
static inline float u2f(unsigned int x) { return *(float *) &x; }

float exp2_fp(int x) {
    unsigned int exp /* exponent */, frac /* fraction */;

    /* too small */
    if (x < 2 - pow(2, Y0) - 23) {
        exp = 0;
        frac = 0;
    /* denormalize */
    } else if (x < Y1 - pow(2, Y2)) {
        exp = 0;
        frac = 1 << (unsigned) (x - (2 - pow(2, Y3) - Y4));
    /* normalized */
    } else if (x < pow(2, Y5)) {
        exp = pow(2, Y6) - 1 + x;
        frac = 0;
    /* too large */
    } else {
        exp = Y7;
        frac = 0;
    }

    /* pack exp and frac into 32 bits */
    return u2f((unsigned) exp << 23 | frac);
}

---

解題想法&思考

• 要解這題，需要充分了解浮點數的邊界判定。
  • 浮點數 = (-1)^s * 2^E * M
  • bias = 2^8-1-1 = 127
  • denormalized:(0~1)
    • s: +/-，本題中沒有負的
    • E: 1-bias= -126 = 2-2⁷
    • M: 0.f₂₂f₂₁f₂₀…f₀
  • normalized:(1~max)
    • s: +/-，本題中沒有負的
    • E: e-bias = e+1-2⁷
    • M: 1.f₂₂f₂₁f₂₀…f₀
• too small
  ​​​​if (x < 2 - pow(2, Y0) - 23) {
  ​​​​    exp = 0;
  ​​​​    frac = 0;
  

  • 判斷d enormalized 可以表達的最小數值，在單精度的情況下為
    0 00000000 00000000000000000000001 即為
    $(-1{)}^0\times 2^{2-2^7}\times 2^{-23}=2^{2-2^7-23}$
    因此， Y0=7
  • x 太小，回傳 0.0 。
• denormalized
  ​​​​} else if (x < Y1 - pow(2, Y2)) {
  ​​​​    exp = 0;
  ​​​​    frac = 1 << (unsigned) (x - (2 - pow(2, Y3) - Y4));
  

  • 判斷 denormalized 可以表達的最大數值為
    0 00000000 11111111111111111111111 即為
    $(-1{)}^0\times 2^{2-2^7}\times \sum _{i=1}^{23}{2}^{-i}=2^{2-2^7}\times (1-2^{-23})$
  • (1-2^-23) 取近似值為 1 ，因此 Y1=2 ， Y2=7 。
  • $2^x=2^{2-2^7}\times (frac\times 2^{-23})$ ，則
    $frac=2^{x-(2-2^7-23)}$
    因此， Y3=7 ， Y4=23 。
• normalized
  ​​​​} else if (x < pow(2, Y5)) {
  ​​​​    exp = pow(2, Y6) - 1 + x;
  ​​​​    frac = 0;
  

  • 判斷 normalized 可以表達的最大數值為
    0 11111110 11111111111111111111111 即為
    $(-1{)}^0\times 2^{2^8-1-1-(2^7-1)}\times (1+1-2^{-23})=2^{2^7-1}\times (2-2^{-23})=2^{2^7}-2^{104}$
  • 最大可表示的值不會超過
    $2^x$ 為
    $2^{2^7}$ ，因此 Y5=7 。
  • 由於要表達 2 的次方，不需要小數點後的表示，因此 frac=0 。
  • $2^x=2^{exp-bias}=2^{exp-(2^7-1)}$ ，則
    $exp=x+(2^7-1)$
    因此， Y6=7 。
• too large
  ​​​​} else {
  ​​​​    exp = Y7;
  ​​​​    frac = 0;
  ​​​​}
  

  • x 太大，回傳
    $+\infty$， INF 表示 exp 為 11111111 = 0xFF ， frac=0 。
    因此， Y7=0xFF 。

延伸問題: 給定一個單精度的浮點數值，輸出較大且最接近的

$2^x$ 值，需要充分考慮到邊界

• 邊界
  • denormalized 最小值:
    0 00000000 00000000000000000000001 =
    $(-1{)}^0\times 2^{2-2^7}\times 2^{-23}=2^{2-2^7-23}$
  • denormalized 最大值:
    0 00000000 11111111111111111111111 =
    $(-1{)}^0\times 2^{2-2^7}\times \sum _{i=1}^{23}{2}^{-i}=2^{2-2^7}\times (1-2^{-23})=2^{lo{g}_2(2^{-126}-2^{-149})}$
  • normalized 最大值:
    0 11111110 11111111111111111111111 =
    $(-1{)}^0\times 2^{2^8-1-1-(2^7-1)}\times (1+1-2^{-23})=2^{2^7-1}\times (2-2^{-23})=2^{lo{g}_2(2^{128}-2^{104})}$
• exp&frac
  • too small
    • 為 0
    • exp = 0
    • frac = 0
  • denormalized
    • exp = 0
    • $2^x=2^{2-2^7}\times (frac\times 2^{-23})$，則
      $frac=2^{x-(2-2^7-23)}$
  • normalized
    • $2^x=2^{exp-bias}\times (1+frac\times 2^{-23})=2^{exp-127}\times (1+frac\times 2^{-23})$
    • $exp=127+x(\mathrm{整}\mathrm{數}\mathrm{部}\mathrm{分})$
    • $frac=(2^{x(\mathrm{小}\mathrm{數}\mathrm{部}\mathrm{分})}-1)\times 2^{23}$
  • too large
    • INF
    • exp = 0xFF
    • frac = 0
• 程式碼

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
static inline float u2f(unsigned int x) { return *(float *) &x; }

float exp2_fp(float x) {
    unsigned int exp /* exponent */, frac /* fraction */;
	float frac_f;

    /* too small */
    if (x < 2 - pow(2, 7) - 23) {
        exp = 0;
        frac = 0;
    /* denormalize */
    } else if (x < log(pow(2, -126) - pow(2, -149))/log(2)) {
        exp = 0;
        frac_f = pow(2.0, (x + 149));
		frac = ((int)(frac_f*2)&1)?(int)frac_f+1:(int)frac_f;//向上捨入
    /* normalized */
    } else if (x < log(pow(2, 128) - pow(2, -104))/log(2)) {
        exp = 127 + (int)((x<0)?x-0.99:x);//x為負數時，無條件進位
        frac_f = ((pow(2.0,x-(float)((int)((x<0)?x-0.99:x)))-1.0)*(1<<23));
		frac = ((int)(frac_f*2)&1)?(int)frac_f+1:(int)frac_f;//向上捨入
    /* too large */
    } else {
        exp = 0xFF;
        frac = 0;
    }

    /* pack exp and frac into 32 bits */
    return u2f((unsigned) exp << 23 | frac);
}

• 執行結果

2^-150
exp2_fp():0.000000000000000000000000000000000000000000000000000000000000
   powf():0.000000000000000000000000000000000000000000000000000000000000
2^-149.5
exp2_fp():0.000000000000000000000000000000000000000000000000000000000000
   powf():0.000000000000000000000000000000000000000000001401298464324817
2^-127
exp2_fp():0.000000000000000000000000000000000000005877471754111437500000
   powf():0.000000000000000000000000000000000000005877471754111437500000
2^-126.3
exp2_fp():0.000000000000000000000000000000000000009547961485342182800000
   powf():0.000000000000000000000000000000000000009547961485342182800000
2^-126
exp2_fp():0.000000000000000000000000000000000000011754943508222875000000
   powf():0.000000000000000000000000000000000000011754943508222875000000
2^127
exp2_fp():170141183460469230000000000000000000000.000000000000000000000000000000000000000000000000000000000000
   powf():170141183460469230000000000000000000000.000000000000000000000000000000000000000000000000000000000000
2^127.5
exp2_fp():240615965050037600000000000000000000000.000000000000000000000000000000000000000000000000000000000000
   powf():240615965050037600000000000000000000000.000000000000000000000000000000000000000000000000000000000000
2^128
exp2_fp():1.#INF00000000000000000000000000000000000000000000000000000000
   powf():1.#INF00000000000000000000000000000000000000000000000000000000
2^0.5
exp2_fp():1.414213538169860800000000000000000000000000000000000000000000
   powf():1.414213538169860800000000000000000000000000000000000000000000

• 結果探討

• 可以從結果看出，基本上與 powf() 的結果相同，除了 -149.5 ，為何 -149~-150 之間仍可以表示?
• (2^0.5)²=2，但 1.41421353816986080²=1.99999993154 ，比實際值略小，如何輸出較大且最接近的
  $2^x$ 值?

這發現很重要！你應該要能從數學的角度計算出上述程式碼和 powf(x, 0.5) 以及理論值的差異，在 IEEE 754 中，已有明確規範。

另外，請改善程式碼縮排，一率用 4 個 space，而非 tab，務必讓版面簡潔清爽。
jserv

第 3 週 測驗 一

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

#include <stdint.h>

/* A union allowing us to convert between a double and two 32-bit integers.
 * Little-endian representation
 */
typedef union {
    double value;
    struct {
        uint32_t lsw;
        uint32_t msw;
    } parts;
} ieee_double_shape_type;

/* Set a double from two 32 bit ints.  */
#define INSERT_WORDS(d, ix0, ix1)       \
    do {                                \
        ieee_double_shape_type iw_u = { \
            .parts.msw = ix0,           \
            .parts.lsw = ix1,           \
        };                              \
        (d) = iw_u.value;               \
    } while (0)

/* Get two 32 bit ints from a double.  */
#define EXTRACT_WORDS(ix0, ix1, d)   \
    do {                             \
        ieee_double_shape_type ew_u; \
        ew_u.value = (d);            \
        (ix0) = ew_u.parts.msw;      \
        (ix1) = ew_u.parts.lsw;      \
    } while (0)

static const double one = 1.0, tiny = 1.0e-300;

double ieee754_sqrt(double x)
{
    double z;
    int32_t sign = 0x80000000;
    uint32_t r, t1, s1, ix1, q1;
    int32_t ix0, s0, q, m, t, i;

    EXTRACT_WORDS(ix0, ix1, x);

    /* take care of INF and NaN */
    if ((ix0 & KK1) == KK2) {
        /* sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN */
        return x * x + x;
    }
    /* take care of zero */
    if (ix0 <= 0) {
        if (((ix0 & (~sign)) | ix1) == 0)
            return x; /* sqrt(+-0) = +-0 */
        if (ix0 < 0)
            return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
    }
    /* normalize x */
    m = (ix0 >> 20);
    if (m == 0) { /* subnormal x */
        while (ix0 == 0) {
            m -= 21;
            ix0 |= (ix1 >> 11);
            ix1 <<= 21;
        }
        for (i = 0; (ix0 & 0x00100000) == 0; i++)
            ix0 <<= 1;
        m -= i - 1;
        ix0 |= (ix1 >> (32 - i));
        ix1 <<= i;
    }
    m -= KK3; /* unbias exponent */
    ix0 = (ix0 & 0x000fffff) | 0x00100000;
    if (m & 1) { /* odd m, double x to make it even */
        ix0 += ix0 + ((ix1 & sign) >> 31);
        ix1 += ix1;
    }
    m >>= 1; /* m = [m/2] */

    /* generate sqrt(x) bit by bit */
    ix0 += ix0 + ((ix1 & sign) >> 31);
    ix1 += ix1;
    q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
    r = 0x00200000;       /* 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 + ((ix1 & sign) >> 31);
        ix1 += ix1;
        r >>= 1;
    }

    r = sign;
    while (r != 0) {
        t1 = s1 + r;
        t = s0;
        if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
            s1 = t1 + r;
            if (((t1 & sign) == sign) && (s1 & sign) == 0)
                s0 += 1;
            ix0 -= t;
            if (ix1 < t1)
                ix0 -= 1;
            ix1 -= t1;
            q1 += r;
        }
        ix0 += ix0 + ((ix1 & sign) >> 31);
        ix1 += ix1;
        r >>= 1;
    }

    /* use floating add to find out rounding direction */
    if ((ix0 | ix1) != 0) {
        z = one - tiny; /* trigger inexact flag */
        if (z >= one) {
            z = one + tiny;
            if (q1 == (uint32_t) 0xffffffff) {
                q1 = 0;
                q += 1;
            } else if (z > one) {
                if (q1 == (uint32_t) KK4)
                    q += 1;
                q1 += 2;
            } else
                q1 += (q1 & 1);
        }
    }
    ix0 = (q >> 1) + 0x3fe00000;
    ix1 = q1 >> 1;
    if ((q & 1) == 1)
        ix1 |= sign;
    ix0 += (m << KK5);
    INSERT_WORDS(z, ix0, ix1);
    return z;
}

---

解題想法與思考

• 程式碼解讀

  ​​​​typedef union {
  ​​​​    double value;
  ​​​​    struct {
  ​​​​        uint32_t lsw;
  ​​​​        uint32_t msw;
  ​​​​    } parts;
  ​​​​} ieee_double_shape_type;
  

  • 建立union，為了能夠分別操作前後32bit，如下所示。

    bit	63	62 ~~~ 52	51 ~~~~~~ 32	31 ~~~~~~~~~~ 0
    value	s	exp	frac	frac
    parts	msw	msw	msw	lsw
    	31	30 ~~~ 20	19 ~~~~~~ 0	31 ~~~~~~~~~~ 0

  ​​​​/* Set a double from two 32 bit ints.  */
  ​​​​#define INSERT_WORDS(d, ix0, ix1)       \
  ​​​​    do {                                \
  ​​​​        ieee_double_shape_type iw_u = { \
  ​​​​            .parts.msw = ix0,           \
  ​​​​            .parts.lsw = ix1,           \
  ​​​​        };                              \
  ​​​​        (d) = iw_u.value;               \
  ​​​​    } while (0)
  ​​​​    
  ​​​​/* Get two 32 bit ints from a double.  */  
  ​​​​#define EXTRACT_WORDS(ix0, ix1, d)   \
  ​​​​    do {                             \
  ​​​​        ieee_double_shape_type ew_u; \
  ​​​​        ew_u.value = (d);            \
  ​​​​        (ix0) = ew_u.parts.msw;      \
  ​​​​        (ix1) = ew_u.parts.lsw;      \
  ​​​​    } while (0)
  

  • 這兩段 define 主要功能是使 union 中的 double 和 parts 互相轉換。
  • #define在做的事，基本上就是代換，將程式其他地方出現 INSERT_WORDS(d, ix0, ix1) ，代換成下面一段程式碼，不做任何修改。通常在使用會以 INSERT_WORDS(d, ix0, ix1); 的樣子出現，因此在代換的程式碼中，最後不會在加 ;。
  • 因為直接代換的關係，可能會產生 dangling else ，所以需要加上 do{...}while(0)。
    • dangling else : 在使用 if-then(-else) 時，產生語意不清的情況，如:
      if A then if B then C else D ，會有兩種解讀:
      if A then (if B then C)else D 或
      if A then (if B then C else D)
    • Linux Kernel 巨集 do{…}while(0) 的撰寫中，舉了一個很好的例子，將其簡化來看:
      ​​​​​​​​​​​​#define func(A, B) \
      ​​​​​​​​​​​​    {              \
      ​​​​​​​​​​​​        if(A)      \
      ​​​​​​​​​​​​            B = 1; \
      ​​​​​​​​​​​​        else       \
      ​​​​​​​​​​​​            B = 0; \
      ​​​​​​​​​​​​    }              \
      
      ​​​​​​​​​​​​int funcA(int A,int B){  
      ​​​​​​​​​​​​    if(B)  
      ​​​​​​​​​​​​        func(A, B);  
      ​​​​​​​​​​​​    else  
      ​​​​​​​​​​​​        printf("Error B");  
      ​​​​​​​​​​​​    return B;  
      ​​​​​​​​​​​​}
      

      將 func(A,B) 用define的內容代換:
      ​​​​​​​​​​​​int funcA(int A,int B){  
      ​​​​​​​​​​​​    if(B)  
      ​​​​​​​​​​​​    {              
      ​​​​​​​​​​​​        if(A)      
      ​​​​​​​​​​​​            B = 1; 
      ​​​​​​​​​​​​        else       
      ​​​​​​​​​​​​            B = 0; 
      ​​​​​​​​​​​​    };
      ​​​​​​​​​​​​    else  
      ​​​​​​​​​​​​        printf("Error B");  
      ​​​​​​​​​​​​    return B;  
      ​​​​​​​​​​​​}
      

      可以發現一個很尷尬的地方，else 前面有 ; ，表示if在else之前就做完了， else 是多的，語意與原本要表達的意思完全不同，這就是dangling else的情況。
      若加上 do{...}while(0) 即可避免這樣的情況。
      ​​​​​​​​​​​​int funcA(int A,int B){  
      ​​​​​​​​​​​​    if(B)  
      ​​​​​​​​​​​​        do{              
      ​​​​​​​​​​​​            if(A)      
      ​​​​​​​​​​​​                B = 1; 
      ​​​​​​​​​​​​            else       
      ​​​​​​​​​​​​                B = 0; 
      ​​​​​​​​​​​​        }while(0);
      ​​​​​​​​​​​​    else  
      ​​​​​​​​​​​​        printf("Error B");  
      ​​​​​​​​​​​​    return B;  
      ​​​​​​​​​​​​}
      

  ​​​​/* take care of INF and NaN */
  ​​​​if ((ix0 & KK1) == KK2) {
  ​​​​    /* sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN */
  ​​​​    return x * x + x;
  ​​​​}
  

  • 考慮 INF 及 NaN 的情況，要判斷是否為這兩者，只要看exp是否全為 1 ，對 ix0 來說， bit 20~30 的位置為 exp ，因此如果是 INF 或 NaN ， ix0 必為
    x 11111111111 xxxxxxxxxxxxxxxxxxx ，需要在意的只有 exp 那一段，因此 KK1 = 0 11111111111 0000000000000000000 = 0x7FF00000
  • 且， exp 必須全為 1 ，因此， KK2 = 0x7FF00000 。
  • 負數開根號為 sNaN，因此，回傳值不能直接回傳 x 。要用-INF得到sNaN，要先了解INF的運算:
    • +INF + +INF = +INF
    • -INF + -INF = -INF
    • +INF + -INF = +INF - +INF = -INF - -INF = sNaN
    • +INF * +INF = +INF
    • -INF * -INF = +INF
    • +INF * -INF = -INF
      基本上，要得到sNaN，就必須使 +INF 和 -INF 相加，而 -INF 得到 +INF 的方法只能乘以自己，因此，回傳值為 x * x + x 才符合要求。
  • 那麼，如果想要不用乘法做回傳值，且完成sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN，有什麼辦法?
    • 先判斷負號，讓 -INF 在進到這個if之前，就因為負號被驅逐了，這樣就可以直接回傳 x 。
  ​​​​/* take care of zero */
  ​​​​if (ix0 <= 0) {
  ​​​​    if (((ix0 & (~sign)) | ix1) == 0)
  ​​​​        return x; /* sqrt(+-0) = +-0 */
  ​​​​    if (ix0 < 0)
  ​​​​        return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
  ​​​​}
  

  • 判斷為 0 或負數，±0 的開根號為 ±0 ，負數的開根號為 sNaN 。可由 ix0 判斷，這也是為什麼 ix0 的 type 要為 int32_t 而非 uint32_t 。
  • x 為 0 ，除了 sign bit 無所謂以外，其他 bit 都為 0 。拆成 msw 及 lsw 兩部份來看: ix0 & (~sign) 即是判斷 msw 中除了 sign bit 以外，其他 bit 是否全為 0 ， ix1 就是lsw全部是否為0。
  • x 為負數，須回傳sNaN。
  ​​​​/* normalize x */
  ​​​​m = (ix0 >> 20);
  ​​​​if (m == 0) { /* subnormal x */
  ​​​​    while (ix0 == 0) {
  ​​​​        m -= 21;
  ​​​​        ix0 |= (ix1 >> 11);
  ​​​​        ix1 <<= 21;
  ​​​​    }
  ​​​​    for (i = 0; (ix0 & 0x00100000) == 0; i++)
  ​​​​        ix0 <<= 1;
  ​​​​    m -= i - 1;
  ​​​​    ix0 |= (ix1 >> (32 - i));//x << i
  ​​​​    ix1 <<= i;
  ​​​​}     
  

  • m = (ix0 >> 20);，表示 m 為 sign + exp 。不過，透過前面的if負數判斷， sign 必定為 0 。
  • 在非正規化(exp = 0)的情況下，先將它變成正規化，透過 shift 來操作。先找出 frac 中，第一個 1 的位置j，再將其移到 exp 的最小位，並修改 m 。也就是說，原本 x = 2^1-bias * f ，現在 x = 2^e-bias * (f * 2^j) ，因此，實際上 x = 2^m * 2^-bias * (f * 2^j) ， m = -j+1，1代表e
  • 假設原本 f 為 0001010…0，視為 0.0001010…0， shift 後， f 為 (1) 010…0 ，視為 1.010…0。不需要再考慮是否加一，這大概是非正規化數 exp 和 frac 設計巧妙的地方。
  ​​​​m -= KK3; /* unbias exponent */
  ​​​​ix0 = (ix0 & 0x000fffff) | 0x00100000;
  

  • 之後，又進一步將 2^-bias 也包進 2^m 之中，同時，修改 ix0 exp 的部分為 0x001。此時， x = x = 2^m * f' (f' = f * 2^j ， m = e - j - bias)
    因此，KK3 = bias = 2^11-1-1 = 1023。
  ​​​​if (m & 1) { /* odd m, double x to make it even */
  ​​​​    ix0 += ix0 + ((ix1 & sign) >> 31);
  ​​​​    ix1 += ix1;
  ​​​​}
  ​​​​m >>= 1; /* m = [m/2] */
  

  • 這段要做 exp 部份的開根號，基本上就是次方除以二。但要考慮到 m 可能是奇數的情況。
  • m 如果是奇數的話，要將 exp+frac 的部分乘以二，以表示會被忽略的部份。
  ​​​​/* generate sqrt(x) bit by bit */
  ​​​​ix0 += ix0 + ((ix1 & sign) >> 31);
  ​​​​ix1 += ix1;
  ​​​​q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
  ​​​​r = 0x00200000;       /* 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 + ((ix1 & sign) >> 31);
  ​​​​    ix1 += ix1;
  ​​​​    r >>= 1;
  ​​​​}
  
  ​​​​r = sign;
  ​​​​while (r != 0) {
  ​​​​    t1 = s1 + r;
  ​​​​    t = s0;
  ​​​​    if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
  ​​​​        s1 = t1 + r;
  ​​​​        if (((t1 & sign) == sign) && (s1 & sign) == 0)
  ​​​​            s0 += 1;
  ​​​​        ix0 -= t;
  ​​​​        if (ix1 < t1)
  ​​​​            ix0 -= 1;
  ​​​​        ix1 -= t1;
  ​​​​        q1 += r;
  ​​​​    }
  ​​​​    ix0 += ix0 + ((ix1 & sign) >> 31);
  ​​​​    ix1 += ix1;
  ​​​​    r >>= 1;
  ​​​​}
  

  • 參考二進制根式法、wikipedia、√Square Roots，了解直式開根號的方法，基本上就對應了這段程式碼。兩個while迴圈，基本上就是對ix0及ix1做開根號的動作。
  ​​​​/* use floating add to find out rounding direction */
  ​​​​if ((ix0 | ix1) != 0) {
  ​​​​    z = one - tiny; /* trigger inexact flag */
  ​​​​    if (z >= one) {
  ​​​​        z = one + tiny;
  ​​​​        if (q1 == (uint32_t) 0xffffffff) {
  ​​​​            q1 = 0;
  ​​​​            q += 1;
  ​​​​        } else if (z > one) {
  ​​​​            if (q1 == (uint32_t) KK4)
  ​​​​                q += 1;
  ​​​​            q1 += 2;
  ​​​​        } else
  ​​​​            q1 += (q1 & 1);
  ​​​​    }
  ​​​​}
  

  • 這一段在做捨入， IEEE 754 預設捨入方法是向偶數捨入，但我看不懂 z 是在做什麼。
  • 在需要進位的情況下，因為[q,q1]實際為一組數字，如果 q1 全為 1 的話，進位之後，會變成 1 0000…00。
  • q1 + 2 之後，q 需要 +1 ，這就表示 KK4 = 1 000…00-2 = 111…10 = 0xfffffffe 。
  • q1 += (q1 & 1); 就比較單純，表示向偶數捨入。
  • 直接跑過程式，看能不能幫助理解。結果，都如我所預測，z = one - tiny 由於 tiny 很小，會被忽略，所以 z == one ，進入第二層 if ，在 q1沒有剛好為 0xffffffff 時，會直接進入第二個 else ， 因為 z = one + tiny; 的 tiny 也會被忽略，所以 z == one 。
  • tiny = 1.0e-300，double可表達的範圍約1.7e-308~1.7e-308
  • z = one - tiny; /* trigger inexact flag */
if (z >= one)的目的是什麼？在什麼情況下會出現 z > 0 ?
  • 在什麼情況下，tiny 不會被忽略？
  ​​​​ix0 = (q >> 1) + 0x3fe00000;
  ​​​​ix1 = q1 >> 1;
  ​​​​if ((q & 1) == 1)
  ​​​​    ix1 |= sign;
  ​​​​ix0 += (m << KK5);
  ​​​​INSERT_WORDS(z, ix0, ix1);
  ​​​​return z;
  

  • exp 要將 1-bias 加回去，因此加上 0x3fe00000 ，此時為 2⁰。
  • m 為開根號後的 exp ，要加回 ix0 ，exp 在20~30的位置，因此，要做 shift ，且 KK5 = 20 。
  • 最後再將 ix0 和 ix1 合成 z 就大功告成了~~~

延伸題目: 改為 float (單精度) 版本，注意應該用更短的程式碼來實作

• 想法與思考
  • 要注意 ix0 和 ix1 的運算，float不需要拆成兩部分。
  • 不需要拆成兩部份，是否還需要union的轉換？
    float 無法做 bitwise 的操作，需要轉換成int才行。
  • 先做負號判斷，再做 +INF and NaN ，這樣可以直接回傳 x 。 -INF 的部分，在負號判斷時，已回傳 sNaN 。
  • bias = 2^8-1 -1 = 127
  • tiny 改為1e-30，因為 float 可表達的範圍約3.4e-38~e.4e-38
• 程式碼

#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 & 0x7ff00000) == 0x7ff00000) {
        /* 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 = 0x01000000;       /* 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) + 0x3f000000;
    ix0 += (m << 23);
	
    INSERT_WORDS(z, ix0);

    return z;
}

• 執行結果

sqrt 1e+78
1.#INF00

sqrt -1e+78
-1.#IND00

sqrt 9
3.000000

sqrt 1.44
1.200000

sqrt 17453723534
132112.546875

因為自動飲料機而延畢的那一年及學習本課程 3 週之後的感想與啟發

• 現實與想像的落差。
  現實其實有更多的層面需要考量，我們時常將事情簡單化，無論是面對自以為龐大的專案，還是認為很簡單的幾個流程，實際上都是需要經過很多的研究時間、成本。我沒想過，一隻十塊錢不到的原子筆、坐上去舒服的椅子，這些生活用品是怎麼做出來的。

• 學用落差。
  我曾經以為，土木系學用落差嚴重只是因為我們不可能真的找一個地方給房子、蓋橋，只能紙上談兵，那牽涉到太多人力、成本。而資訊系這方面的問題，比較沒那麼嚴重，畢竟寫程式每個人的電腦都可以做，而且馬上就可以看到結果。實際上不是，寫寫小程式不難，就像土木系在做模型時一樣，但面對大一點的就容易出現問題了。

• 全心投入。
  看完這片文章，讓我想起一年前，為了做一個小房子的模型，而做過無數的實驗跟研究。那時是要蓋一個混凝土的模型，原本的想像也是很簡單，設計尺寸、做模、灌水泥、拆模。但實際下去做之後，才發現現實世界有太多問題，那時為了做這件事，整天都在想、在試。做完之後，成果也不怎麼樣，最後也是收到角落。而這之中，最寶貴的，我想，不是結果多厲害，而是在那段時間全心投入的自己。在實習過後，體會更深的，工作之餘，想要再將心力投入到其他事情上，是沒有那麼容易的，尤其是需要耗費腦力的事。希望自己能一直保有那份熱情。

• 關於課程
  一開始覺得壓力很大，這種上課模式沒有遇過 : 每週都小考，還有很多課後學習的教材要看，作業還會被同學批評、觀摩 (害羞)，一週花費超過16個小時在上面，不過，還是很多教材沒看完。但，這好像比較接近我想像中的大學課程，有點半自學的性質，也可以說很幸運的在最後一個學期修到這門課。

• 程式寫得好的人，數學都很厲害。
  這三、四週下來，深深的體悟到什麼是真正在用的程式。其實很多事情的原理並不難，只是我沒有好好理解過，隨堂考的題目要解出來其實不難，但是在沒有細細看過程式前，很多都是靠感覺猜的，因為那些背後的理論我尚未融會貫通。但是，要完全看懂程式碼在做什麼，為什麼要那樣做就完全是另外一回事了，這大概就是理論跟實做的差距吧。

「課後學習」教材和 CS:APP 3/e心得和提問