2020q3 Homework5 (quiz5)

contributed by < OscarShiang >

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

我們可以從計算 od 與 divop 的 slots 進行的位移運算知道

D1 應該會是

D1 = 2

因為 >> 與 << 這些 bitwise 的運算子只適用在 int 的資料型態上，所以如果我們要利用 2 化簡除數與被除數的時候就會需要使用 orig / 2 的方式來進行

而當我們將兩者化簡到不能再以 2 來化簡的時候，就需要進行逼近的方式來進算，因此在 D2 的位置就嘗試透過將 slots 加一使其改成偶數的方式來逼近 orig / slots 的浮點數數值。

最後因為 od == 1 的時候我們採取的手法是用逼近的方式，所以最後我們需要利用 result += divop(result, slots) 這行指令來做數值的修正。

最後以 float 的除法與 divop 來比較兩者在計算

$1÷3$ 的表現

測試程式的行為如下所示

int main(void)
{
    printf("%.30f\n", divop(1, 3));
    printf("%.30f\n", (float) 1 / 3);
    return 0;
}

$ ./precise
divop:     0.333333333333333314829616256247
float div: 0.333333343267440795898437500000

在精準度的表現上，divop 優於浮點數的除法

為了測試兩者在執行時間上的差異，我將測試程式修改為下列

int main(void) 
{
    struct timespec start, end;

    float tmp;

    clock_gettime(CLOCK_MONOTONIC, &start);
    tmp = divop(1, 3);
    clock_gettime(CLOCK_MONOTONIC, &end);

    printf("divop(1, 3): %ld\n", end.tv_nsec - start.tv_nsec);

    clock_gettime(CLOCK_MONOTONIC, &start);
    tmp = (float) 1 / 3;
    clock_gettime(CLOCK_MONOTONIC, &end);

    printf("float div: %ld\n", end.tv_nsec - start.tv_nsec);
    return 0;
}

測量的結果如下，

$ ./precise
divop(1, 3): 22807
float div: 37

從結果中可以看到雖然使用 divop 計算得精準度較浮點數除法好，但是在執行時間上因為 divop 需要透過函式不斷遞迴呼叫的方式來逼近計算數值，所以在執行時間上大幅落後浮點數運算。

測驗 2

本題嘗試以位元與整數運算來計算出平方根

在這個實作中可以分為幾個階段

1. 去除負數與 NaN 等特殊值
2. normalize
3. 利用牛頓法計算平方根

為了能夠使用位元運算來判斷該變數的各個位元的狀況，在實作中使用 macro 來實作這樣的功能，將 double 變數的 bit pattern 轉換到 int 變數中

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

其原理是使用 union 來產生可以以 int 或是 float 的方式來存取的變數，即可將 float 的 bit pattern 取出來存在 int 變數中。

QQ2

本題所在位置在

float ieee754_sqrt(float x)
{
    
    <...>
    
    ix0 = (q >> 1) + QQ1;
    ix0 += (m << QQ2);
	
    INSERT_WORDS(z, ix0);

    return z;
}

這個部分主要就是將我們先前在

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

所計算出來的 exponent 部分寫入到新的 float 變數中。而 exponenet 的部分對應到 IEEE 754 floating point number 的定義

所以我們需要將 m 的數值放在上圖綠色的位置中，因此 QQ2 的答案就是

QQ2 = 23

測驗 3

這題的用意是在實作一個可以計算給定數值能夠表示成幾種連續整數和的程式

例如 9 可以表示成

$$9=9=4+5=2+3+4$$

總共三種可能

而 15 則可以表示成

$$15=15=8+7=4+5+6=1+2+3+4+5$$

因此總共有 4 種表示方法

參考資料

tags: sysprog2020