2020q3 Homework5 (quiz5)

contributed by < ccs100203 >

tags: linux2020

第 5 週測驗題

測驗 1

浮點數除法程式: (fdiv.c)
假設 divop() 的第二個參數必為大於 0 的整數，而且不超過 int 型態能表達的數值上界。

#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 / 2, od ? (slots + 1) >> 1 : slots >> 1);
    if (od)
        result += divop(result, slots);
    return result;
}

• 利用遞迴的方式做除法。簡單來說程式會不斷的對被除數與除數同除 2，直到符合終止條件就會把 result 傳回去。但遇到 slots (除數) 是奇數時，因為會在第 8 行視為偶數直接運算，所以藉由第 10 行補上誤差。

• 第 5 行的 if 就是作為遞迴的終止條件，當 slots 等於 1 時就會把當前的 orig (也就是當前運算的 result) return 回去上一層。而 orig == 0 的判斷是用來處理 orig (被除數)是 0，或者是數字已經小到精度不夠，程式判定他為 0 時，也會結束並回傳。

• 第 8 行會直接做

  $÷2$ 並作為參數進入 recursive 的下一層。注意到他有對 od 做特別處理，也就是對 odd (奇數) 會特別做處理。當 slots 是奇數時會先做 +1 才除 2，也就是無條件進位。

• 但這樣子在 slots 是奇數時，數值會比實際答案小，因為除數無條件進位了，所以在第 9~10 行就是用來彌補奇數時的值，也就是會取近似值，不斷的把值透過遞迴呼叫加到 result，維持除法精準度。

假設

$10.0/3$，程式會把 slots 為 1 的結果全部相加起來，得以近似 3.3333…

orig: 2.500000,	slots: 1
orig: 0.625000,	slots: 1
orig: 0.156250,	slots: 1
orig: 0.039062,	slots: 1
orig: 0.009766,	slots: 1
orig: 0.002441,	slots: 1
orig: 0.000610,	slots: 1
...

寫成數學式的話會變成:

$$\sum _{i=1}^{\mathrm{\infty }}\frac{10}{4^i}=\frac{10}{3}$$

換句話說，就是利用等比級數的方式去算出誤差值，會計算出

$S_n$，將這些結果累加求得最後的 result

$$S_n=\frac{a(1-r^n)}{1-r}$$

數學推導，假設 A 除以 n

$$\begin{gathered} \frac{A}{n}=\frac{A(n+1)}{n(n+1)}=\frac{An+A}{n(n+1)} \\ =\frac{A}{n+1}+\frac{A}{n(n+1)} \\ =\frac{A}{n+1}+\frac{A}{n}\ast \frac{1}{n+1} \\ =\frac{A}{n+1}+\frac{1}{n+1}(\frac{A}{n+1}+\frac{A}{n}\ast \frac{1}{n+1}) \end{gathered}$$

這段式子可以不斷延伸下去，就會發現符合等比級數的數學式

其中在 n 屬於偶數時會使用

$÷2$ 的方式去處理，式子就會類似下列情況，還是符合推導

$$\frac{A÷2}{n÷2}\vee \frac{A÷2}{(n+1)÷2}$$

測驗 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 = 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;
}

關於 do…while 的解釋，請參考 為什麼需要 do {…} while(0)

程式運作原理

• EXTRACT_WORDS 輸入 float，藉由 union 的方式得到 int
• INSERT_WORDS 輸入 int，藉由 union 的方式得到 float
  關於 union 的說明參考這裡

x is less than or equal 0

根據 Division by zero

In IEEE 754 arithmetic, a ÷ +0 is positive infinity when a is positive, negative infinity when a is negative, and NaN when a = ±0. The infinity signs change when dividing by −0 instead.

如果 x 是 -0 的話就直接 return，其餘則是 -NaN

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

x is NaN or +INF

根據 IEEE 754 的規範，Exponent 的部分全為 1 時會是 INF 或 NaN。

程式遇到 NaN 或 +INF 時會直接 return

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

x is denormalized

根據 IEEE 754 的規範，Exponent 為 0 時，可能遇到 subnormal 的情況

m 代表的就是 Exponent 的部分，如果 m 是 0，那程式會強制 normalize，會從

$fraction$ 找最高位的 1，不斷對

$fraction$ 往左 shift，一直到補進了 normalized 的 1，再將 shift 多少位數從 Exponent 的地方扣掉。

/* normalize x */
m = (ix0 >> 23);
if (m == 0) { /* subnormal x */
    for (i = 0; (ix0 & 0x00800000) == 0; i++)
        ix0 <<= 1;
    m -= i - 1;
}

x is normalized and Greater than 0

這時候的 x 一定是 normalized 的數字並且大於 0。
- m -= 127 是把 bias 的部分扣除，得到真正的 Exponent
- ix0 存放的是

$significand$ 也等同

$fraction+1$
- 再來會做 m >>= 1，意即

$m=\lfloor m/2\rfloor$，若 m 是奇數的話則會將多除掉的部分補到 ix0 內，所以才會做 ix0 += ix0。

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

再來的部分就是做 sqrt，關於直式開 2 次方根，參考 Binary numeral system (base 2) 以及 直式開二次方根(page 2~5) 的做法。

         q (r)
       ----------
  s0  / ix0
  (t)
  
 t: 用來協助 s0 的紀錄，以及判斷當前的 ix0 是否大於 t
 r: 用來協助 q 的紀錄，為當前正在處理的 bit。
    從 0x01000000 開始是因為 ix0 最高位為 1 的 bit 會在第 25 (可能因為 m 為奇數而在第 26)

解說 while loop 運作原理:

1. 在直式運算中，會先從右到左，將 ix0 每兩個 bit 劃分為一組
2. 每一次會先將 r (當前處理的 bit) 與 s0 (直式左邊的數字) 相加存入 t
3. 從直式運算的規則可以得知若 ix0 >= t 則從 ix0 扣掉 t，q (商) 則會增加 r。而根據直式運算 2 次方根的規則，s0 應更新為上次 s0 + 2 * r，由於 t = s0 + r，故程式內的寫法成立。
4. 若 ix0 < t 或執行完上述第 3 點，在直式規則中應一次從 ix0 新增兩個位數來處理，但這裡利用做 r >>= 1 以及 ix0 += ix0，恰巧符合運算規則，能夠接著運算 ix0 下兩個 bit。這樣的好處是不需要在 ix0 < t 時對 s0 做額外補 0 的動作，並且也不用移動 s0 來對齊下一步的 ix0。在 q 上也會因為 r 一次只往右 shift 1，符合其移動規則，不須額外調整。

/* 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;
}

Rounding

參考 Rounding、Rounding modes、Rounding rules
我想程式是利用 one 與 tiny 的運算去判斷現在的系統屬於哪種 Rounding Mode，根據當前的 rounding mode 做出符合規範的 rounding。

• one - tiny < 1.0，代表當前的 mode 不會做進位，所以不需要特別處理 q，因為之後往右 shift 就會自動捨去
• one - tiny >= 1.0，代表當前的 mode 會進位，但要判斷屬於無條件進位或是四捨五入
  • 藉由 one + tiny > 1.0 可以知道系統屬於無條件進位，所以直接進位，做 q += 2
  • 藉由 one + tiny <= 1.0 可以知道系統屬於四捨五入，那只會在 q & 1 為 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);
    }
}

將數值轉換回 float

• 由於 r 是從第 25 位開始算，所以 q 需要往右 shift 1 才能夠對齊 IEEE 754 規範中
  $significand$ 的位置。
• 然後再加上 126 來符合 bias。 (原本是 127，但已經從
  $significand$ 得到 1，所以只需 126)
• 把先前的 m (Exponent) 往左 shift 23，將他放在 IEEE 754 規範的位置。

ix0 = (q >> 1) + 0x3f000000;
ix0 += (m << 23);

LeetCode 69. Sqrt(x)

TODO

測驗 3

LeetCode 829. Consecutive Numbers Sum
給定一個正整數

$N$，問

$N$ 能寫成多少種連續正整數之和，例如

$9$ 可寫成

$4+5$，或

$2+3+4$。

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

我看不懂題目上的推導敘述，所以重新從等差數列的方法去理解程式。
程式的作法會在每一次的迴圈做 x += 1 - i 然後再去判斷 x 能不能整除 i

• 用下表說明 (程式內的 i 是從 2 開始的)
  i 即為程式內的 i，代表輸入的數字 N 可以相等於被 i 個連續數字累加
  n 則是代表連續數值累加中的第一個數字。

1-i	-0	-1	-2	-3	-4	-5
i	1	2	3	4	5	6
n	n	n+1	n+2	n+3	n+4	n+5

1. 當 i 為 2 的時候，代表

   $n+(n+1)=N$，而
   $+1$ 的部分藉由 1 - i 扣掉了，所以實際上是判斷
   $(N-1)\mathrm{mod}2$ 是否為 0 (可否整除)，這樣便可以判別是否能由 2 個連續數字累加為
   $N$。

2. 再來當 i 為 3，代表

   $n+(n+1)+(n+2)=N$，新的
   $+2$ 也在這次的 1 - i 扣掉，
   $+1$ 則是在上一次的迴圈就扣掉了，所以判斷式變成
   $(N-1-2)\mathrm{mod}3$ 是否為 0。從這邊我們可以往後推論每個 i 的情況，而迴圈一直執行到 i >= x 會結束，因為 x 此時無法繼續扣除 1 - i，就算兩者相等扣完之後 x 也只剩下 1，不可能整除 i。

3. 程式內的 i 是從 2 開始的，原因在於 i 等於 1 是絕對會成立的情況，所以節省了一次計算。