2020q3 Homework (quiz5)

contributed by < shauming1020 >

tags: homework

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

假設 divop() 的第二個參數必為大於 0 的整數，而且不超過 int 型態能表達的數值上界。

１.　解釋上述程式碼運作原理

• int od = slots & 1; 判斷 slots 是否為奇數
• double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);
  • od ? (slots + D2) >> 1 : slots >> 1) 當 slots 是奇數時，必須 + 1 後再進行右移 (/2) ，如此才會滿足遞迴中止條件 (slots == 1)，因此 D2 為 1
  • orig / D1 浮點數無法透過右移 bit 來進行除 2 的動作，因此 D1 為 2
  • 將 orig 表示為分數，分子與分母同除 2 後不會影響最後結果，當分母為 1 時，分子即為商數
• if (od) result += divop(result, slots);
  • 當 slots 為奇數時，slots 加 1 後再去除 2 會產生誤差，導致求出的商數比原商數小

    e.g.

    $1÷3=0.333$ 在計算時會變成

    $1÷4=0.250$，誤差為

    $0.083$，恰近似於

    $0.250/3$
    將上述例子用數學式描述

    ${\displaystyle \frac{1}{3}}={\displaystyle \frac{1}{4}}+{\displaystyle \frac{{\displaystyle \frac{1}{4}}}{3}}={\displaystyle \frac{1}{4}}+{\displaystyle \frac{1}{12}}={\displaystyle \frac{4}{12}}$

  • 因此需要將誤差給加上去，而 slots 為偶數時，可以直接同除 2，不會有誤差的產生
  • 所以每次求誤差時會一直遞迴執行程式，e.g. 求
    ${\displaystyle \frac{{\displaystyle \frac{1}{4}}}{3}}$ 時會再去計算誤差

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

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

觀察 gcc -S -O 編譯器行為

• Recursive

int fib(int n) {
    if (n == 0) return 0;
    if (n == 1) return 1;
    return fib(n - 1) + fib (n - 2);
}
main()
{
    fib(100);
}

• Assembly code

...
fib:
.LFB0:
	.cfi_startproc
	pushq	%rbp # Push to stack
	.cfi_def_cfa_offset 16
	.cfi_offset 6, -16
	pushq	%rbx # Push to stack
	.cfi_def_cfa_offset 24
	.cfi_offset 3, -24
	subq	$8, %rsp
	.cfi_def_cfa_offset 32
	movl	%edi, %ebx
	testl	%edi, %edi
	je	.L2
	cmpl	$1, %edi
	jne	.L4
.L2:
	movl	%ebx, %eax
	addq	$8, %rsp
	.cfi_remember_state
	.cfi_def_cfa_offset 24
	popq	%rbx
	.cfi_def_cfa_offset 16
	popq	%rbp
	.cfi_def_cfa_offset 8
	ret
.L4:
	.cfi_restore_state
	leal	-1(%rdi), %edi
	call	fib
	movl	%eax, %ebp
	leal	-2(%rbx), %edi
	call	fib
	leal	0(%rbp,%rax), %ebx
	jmp	.L2
	.cfi_endproc
...
main:
.LFB1:
	.cfi_startproc
	subq	$8, %rsp
	.cfi_def_cfa_offset 16
	movl	$100, %edi
	call	fib
	movl	$0, %eax
	addq	$8, %rsp
	.cfi_def_cfa_offset 8
	ret
	.cfi_endproc
...

• Tail recursion

int fib(int n, int a, int b) {
    if (n == 0) return a;
    return fib(n - 1 , b, a + b);
}
main()
{
    fib(100, 0, 1);
}

---

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

1. 解釋上述程式碼何以運作

搭配研讀 以牛頓法求整數開二次方根的近似值 (第 19 頁)

---

測驗 3

$kx=N-{\displaystyle \frac{(k-1)k}{2}}$，由於 x 要為正整數，因此

$x=(N-{\displaystyle \frac{(k-1)k}{2}})/k\in Z^{+}$，即

$(N-{\displaystyle \frac{(k-1)k}{2}})\mathrm{\%}k=0$，而 k 的範圍可從近似解

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

1. 解釋上述程式碼何以運作

• 撰寫與推導後的數學式相同的程式碼

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

• 觀察

  $S_k=N-{\displaystyle \frac{(k-1)k}{2}},1\le k<\sqrt{2N}$

  $k$	$S_k$	$S_k-S_{k-1}$
  1	N	0
  2	N + (-1)	-1
  3	N + (-3)	-2
  4	N + (-6)	-3
  5	N + (-10)	-4
  n	N + ${\displaystyle \frac{(-n+1)n}{2}}$	1-n

• 檢驗

  $S_k\mathrm{\%}k=0$ 是否成立來判斷正整數解 x，從上表可以觀察出
  $S_k=S_{k-1}+(1-k)$，故 ZZZ 為 (1 - i)

  e.g.

  $S_5=N+(-10)=S_4+(1-5)=N+(-6)+(-4)$

• 而

  $S_k\mathrm{\%}k=0$，且 k 為大於等於 2 的正整數，因此
  $2\le k\le S_k$
  • 在 for 迴圈內的 x += (1-i) 計算出
    $S_k$ 後，進入下個迴圈前，k 會先被加一再和
    $S_k$ 比較大小，因此範圍為
    $(k+1)<S_k$，迴圈的條件判斷式才會是 i < x

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

---