2020q3 Homework5 (quiz5)

contributed by < JimmyLiu0530 >

tags: 進階電腦系統理論與實作

測驗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 型態能表達的數值上界。請補完程式碼。

程式說明

原理:

${\displaystyle \frac{N}{d}}$ =

${\displaystyle \frac{N}{d+1}}\ast {\displaystyle \frac{d+1}{d}}$ =

${\displaystyle \frac{N}{d+1}}\ast (1+\frac{1}{d})$ =

${\displaystyle \frac{N}{d+1}}+{\displaystyle \frac{N}{d(d+1)}}$

e.g.

${\displaystyle \frac{N}{5}}={\displaystyle \frac{N}{5+1}}+{\displaystyle \frac{N}{5(5+1)}}$

Line 8 的作用在於:

若除數 slots 是偶數，那麼被除數 orig 與除數 slots 可以先同時除以 2，再進行運算，因此 D1 = 2。
若 slots 為奇數，則根據上述的原理

${\displaystyle \frac{orig}{slots}}$ =

${\displaystyle \frac{orig}{slots+1}}+{\displaystyle \frac{orig}{slots(slots+1)}}$ 告訴我們兩件事:

1. 首先就是 D2 = 1。又因為 slots 原本是奇數加 1 後變成偶數，所以也可以與 orig 同除以 2。
2. 再來就是最終結果 result 除了計算
   ${\displaystyle \frac{orig}{slots+1}}$ 外還要再加上
   ${\displaystyle \frac{orig}{slots(slots+1)}}$ ，因此 line 10 需要 divop(result, slots)。

延伸問題

1. 以編譯器最佳化的角度，推測上述程式碼是否可透過 tail call optimization 進行最佳化，搭配對應的實驗;

搭配閱讀 C 語言: 編譯器和最佳化原理 及 C 語言: 遞迴呼叫

2. 比照 浮點數運算和定點數操作，改寫上述程式碼，提出效率更高的實作，同樣避免使用到 FPU 指令 (可用 gcc -S 觀察對應的組合語言輸出)，並與使用到 FPU 的實作進行效能比較

---

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

Reference:

• Floating point division and square root algorithms and implementation in the AMD-K7/sup TM/ microprocessor
• √Square Roots
• 以牛頓法求整數開二次方根的近似值 (第 19 頁)

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

程式說明

• 牛頓法求整數開根號近似值的核心概念:
  • 假設估算
    $\sqrt{N}$，則令
    $f(x)=x^2-N$，利用
    $f(x)=0$ 求得
    $x=\pm \sqrt{N}$。
  • 接著，取一個估算值 a，對
    $(a,f(a))$ 做切線，並與 x 軸相交於 b。一樣再對
    $(b,f(b))$ 做切線與 x 軸相交於 c，重複上述步驟即可逼近
    $\sqrt{N}$。

延伸問題

1. 嘗試改寫為更有效率的實作，並與 FPU 版本進行效能和 precision /accuracy 的比較

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

---

測驗3: LeetCode 829. Consecutive Numbers Sum

給定一個正整數 N，問 N 能寫成多少種連續正整數之和，例如 9 可寫成

$4+5$，或者

$2+3+4$。由於要寫成連續正整數之和，則可用等差數列來思考，且差值為 1，這個等差數列不必從 1 開始，假設從 x 開始，且個數共有 k 個，則可寫出該等差數列為：

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

令其和為 N，根據等差數列的求和公式，可寫出以下等式：

$kx+{\displaystyle \frac{(k-1)k}{2}}=N$

整理後可得：

$kx=N-{\displaystyle \frac{(k-1)k}{2}}$

對任意 k 值，倘若 x 能得到正整數解，就表示必定會有個對應的等差數列和為 N。由於 k 是等差數列的長度，首先必定大於 0，即為下限。由於 x 也必是正整數，可得:

$N-{\displaystyle \frac{(k-1)k}{2}}>0$

從而得到近似解：

$k<\sqrt{2N}$

參考實作如下:

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

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

程式說明

根據題目的敘述及提示可知

$k$ 為正整數且

$k<\sqrt{2N}$，有了 k 值的範圍後，在這範圍內一一去檢驗

$kx=N-{\displaystyle \frac{(k-1)k}{2}}$ 是否成立就可以知道 N 能不能寫成 k 個連續整數相加。

因為任何正整數一定能寫成 1 個數相加(即自己)，所以直接從 k=2 開始看。

接著用一個例子來說明如何將數學式轉換成程式碼: (N=9為例)

k	$kx=N-{\displaystyle \frac{(k-1)k}{2}}$	x	是否能寫成 k 個連續整數相加
2	$2x=9-{\displaystyle \frac{1\ast 2}{2}}=8$	4	:O:
3	$3x=9-{\displaystyle \frac{2\ast 3}{2}}=6$	2	:O:
4	$4x=9-{\displaystyle \frac{3\ast 4}{2}}=3$	${\displaystyle \frac{3}{4}}$	:X:

根據上面例子會發現

1. 每一次迭代，
   $N-{\displaystyle \frac{(k-1)k}{2}}$ 的值都會比上一次等差地減少。
   • e.g. k=2 時，值為 8 會是上一輪的 9 減掉 1; k=3 時，值為 6 會是上一輪的 8 減掉 2, 以此類推，因此 ZZZ = 1-i。
2. $2x=8$ 有整數解釋因為 8 被 2 整除，然而
   $4x=3$ 無整數解是因為 3 不被 4 整除，因此可以推廣到
   $kx=m$ 有整數解就代表 m 要能被 k 整除，呼應 line 9 的判斷式。

• 執行效能:
  
  

延伸問題

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

int consecutiveNumbersSum(int N){
    int res = 1, count;
        while (N % 2 == 0) N /= 2;
        for (int i = 3; i * i <= N; i += 2) {
            count = 0;
            while (N % i == 0) {
                N /= i;
                count++;
            }
            res *= count + 1;
        }
        return N == 1 ? res : res * 2;
}

TODO: 透過 shift 改寫

jserv

• 此方法的原理：
  由上面的推論可得

  $kx=N-{\displaystyle \frac{(k-1)k}{2}}$ ，變換一下得到
  $2N=k(2x+k-1)$，因為
  $2N$ 是偶數，代表著
  $k$ 或是
  $(2x+k-1)$ 要是偶數，而當 k 偶數，
  $(2x+k-1)$ 就會是奇數; 當 k 奇數，
  $(2x+k-1)$ 就會是偶數。

  也就是說，N 的每一個 odd factor 都會對應一組答案。因此，此問題就變成了 N 有幾個 odd factor。

• 計算 odd factor 數量:

  若

  $N=3^a\ast 5^b\ast 7^c\ast ...$，則 N 共有 (a+1)(b+1)(c+1)… 個因數
  ( 因為在 line 3 會將所有2因數移除，所以才可以如此假設 N )

• 程式碼說明:
  上面這個新的程式碼就是用來找 N 有多少個 odd factor。首先，因為我們只關心 odd factor，所以先移除 N 的所有 2 因數，再對所有小於 N 的奇數去看是否為 N 的因數。若是，則去計算該因數有多少次方(即 count)，加 1 後再乘入 res 中。

  最後，for 迴圈結束檢查若 N 為 1，那麼直接回傳 res; 否則代表 N 中還剩一個最大質因數 P，且其次方必為1，因此回傳 res * (1+1)。

  因為 line 4 for 迴圈的中止條件為 i * i <= N，所以執行完 for 迴圈後 N 有可能剩下的就只有其最大質因數 P (因為 (1) 若為合數，那麼早在 for 迴圈就被除掉了。(2) 假設 N 還存在另一個質因數 Q，Q < P，則代表 for 迴圈沒有執行到 i = Q，即 Q * Q > N。又

  $N\ge P\ast Q>Q\ast Q$，很明顯矛盾，因此 N 中不存在另一個質數)，使得 P * P > N，也因為 P * P > N，所以 P 的次方必為 1。

• 執行效能: