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
$\neq$ 0
令 k 為
$\leq n$ 的正整數

$\begin{aligned} a_{n} & = \sum_{i = 0}^{n} a_{i} \times 2^{i} \\ = [a_{0} + a_{2} (3 + 1) + a_{4} (15 + 1) + . . . + a_{2 k} (2^{2 k} - 1 + 1) + . . .] \\ + [a_{1} (3 - 1) + a_{3} (9 - 1) + a_{5} (33 - 1) + . . . + a_{2 k + 1} (2^{2 k + 1} + 1 - 1) + . . .] \\ = (a_{0} + a_{2} + a_{4} + . . . + a_{2 k} + . . .) \\ + [a_{2} (3) + a_{4} (15) + . . . + a_{2 k} (2^{2 k} + 1) + . . .] \\ - (a_{1} + a_{3} + a_{5} + . . . + a_{2 k + 1}) \\ + [a_{1} (3) + a_{3} (9) + a_{5} (33) + . . . + a_{2 k + 1} (2^{2 k + 1} + 1) + . . .] . . . 第一式 \end{aligned}$
利用歸納證明法，可證明 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)。寫成算式為

$N \div 4 = a . . . b => N = a \times 4 + b = a \times 5 - a + b => 如果 b - a 為 5 的倍數，則 N 為 5 的倍數$
程式碼

#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 的次方的關係。
$N = a \times 8 + b = a \times 7 + (a + b) => 如果 a + b 為 7 的倍數，則 N 為 7 的倍數$
程式碼

#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 改寫。
改寫一:

$N = a \times 4 + b = a \times 3 + (a + b) => 如果 a + b 為 3 的倍數，則 N 為 3 的倍數$
程式碼

#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));
}

改寫二:

$N = a \times 2 + b = a \times 3 + (- a + b) => 如果 b - a 為 3 的倍數，則 N 為 3 的倍數$
程式碼

#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 (f r a c \times 2^{- 23})$ ，則
  $f r a c = 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^{e x p - b i a s} = 2^{e x p - (2^{7} - 1)}$ ，則
  $e x p = 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^{l o 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^{l o 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 (f r a c \times 2^{- 23})$ ，則
    $f r a c = 2^{x - (2 - 2^{7} - 23)}$
- normalized
  - $2^{x} = 2^{e x p - b i a s} \times (1 + f r a c \times 2^{- 23}) = 2^{e x p - 127} \times (1 + f r a c \times 2^{- 23})$
  - $e x p = 127 + x (整數部分)$
  - $f r a c = (2^{x (小數部分)} - 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，務必讓版面簡潔清爽。
:notes: 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

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個小時在上面，不過，還是很多教材沒看完。但，這好像比較接近我想像中的大學課程，有點半自學的性質，也可以說很幸運的在最後一個學期修到這門課。
程式寫得好的人，數學都很厲害。
這三、四週下來，深深的體悟到什麼是真正在用的程式。其實很多事情的原理並不難，只是我沒有好好理解過，隨堂考的題目要解出來其實不難，但是在沒有細細看過程式前，很多都是靠感覺猜的，因為那些背後的理論我尚未融會貫通。但是，要完全看懂程式碼在做什麼，為什麼要那樣做就完全是另外一回事了，這大概就是理論跟實做的差距吧。

2018q1 Homework2 (assessment)

第 1 週 測驗 一

解題想法&思考

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

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

第 2 週 測驗 一

解題想法&思考

延伸問題: 給定一個單精度的浮點數值，輸出較大且最接近的 2x 值，需要充分考慮到邊界

第 3 週 測驗 一

解題想法與思考

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

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

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

第 1 週測驗 `一`

第 2 週測驗 `一`

延伸問題: 給定一個單精度的浮點數值，輸出較大且最接近的
$2^{x}$ 值，需要充分考慮到邊界

第 3 週測驗 `一`