Try   HackMD

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 型態能表達的數值上界。

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

  • 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

      將上述例子用數學式描述
      13=14+143=14+112=412

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

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(k1)k2,由於 x 要為正整數,因此
x=(N(k1)k2)/kZ+
,即
(N(k1)k2)%k=0
,而 k 的範圍可從近似解
k<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; }

  • 觀察

    Sk=N(k1)k2,1k<2N

    k
    Sk
    SkSk1
    1 N 0
    2 N + (-1) -1
    3 N + (-3) -2
    4 N + (-6) -3
    5 N + (-10) -4
    n N +
    (n+1)n2
    1-n
  • 檢驗

    Sk%k=0 是否成立來判斷正整數解 x,從上表可以觀察出
    Sk=Sk1+(1k)
    ,故 ZZZ 為 (1 - i)

    e.g.

    S5=N+(10)=S4+(15)=N+(6)+(4)

  • Sk%k=0,且 k 為大於等於 2 的正整數,因此
    2kSk

    • 在 for 迴圈內的 x += (1-i) 計算出
      Sk
      後,進入下個迴圈前,k 會先被加一再和
      Sk
      比較大小
      ,因此範圍為
      (k+1)<Sk
      ,迴圈的條件判斷式才會是 i < x

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