2020q3 Homework5 (quiz5)

contributed by < RusselCK >

tags: `RusselCK`

浮點數運算＆數值系統＆ Bitwise 練習題

Q1. 浮點數除法程式: (`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);
                            // D1 = 2; D2 = 1;
    if (od)
        result += divop(result, slots);
    return result;
}

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


double divop(double orig, int slots){

依據題意， orig 為被除數， slots 為除數
- orig 的型態為雙精度浮點數 ( double )


    if (slots == 1 || orig == 0)
        return orig;

當除數為 1 或被除數為 0 時，除法結果還是被除數


    int od = slots & 1;
    double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);

運用 quiz4 Q4 的概念：
- 若被除數與除數都有 2 因子，可將兩者都先除以 2 ，除法結果不變
- 浮點數可以一直除以 2
od 判斷除數 slots 是否為偶數
- 若 slots 為偶數，則 od 為 0
- 若 slots 為奇數，則 od 為 1
若除數 slots 為偶數，則被除數 orig 與除數 slots 同除以 2，除法結果不變
- divop(double orig, int slots) = divop(orig / 2, slots >> 1)
- 因此， D1 = 2




    double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);

    if (od)
        result += divop(result, slots);

先看一個數學推導：

$A : f l o a t i n g n u m b e r, n : o d d i n t e g e r$

$\begin{aligned} A \div n = \frac{A}{n} & = \frac{A}{n + 1} \times \frac{n + 1}{n} \\ = \frac{A}{n + 1} \times (1 + \frac{1}{n}) \\ = \frac{A}{n + 1} + \frac{A}{n + 1} \times \frac{1}{n} \\ = [A \div (n + 1)] + {\begin{array}{c} [A \div (n + 1)] \div n \end{array}} \\ = [\frac{A \div (n + 1)}{2}] + {\begin{array}{c} [\frac{A \div (n + 1)}{2}] \div n \end{array}} \end{aligned}$
若除數 slots 為奇數，搭配上述推導與先處理除以 2 的概念
- divop(double orig, int slots) = result + divop(result, slots)
  - result = divop(orig / 2, (slots + 1) >> 1)
- 因此， D2 = 1

Q2. 浮點數開二次方根的近似值

IEEE 754 Standard: Scope
浮點數與定點數
假設 float 為 32-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式

參考 e_sqrt.c

































































































#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 & 0x7ff00000) == 0x7ff00000) {
        /* 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 */
    // QQ0 = 0x01000000

    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; // QQ1 = 0x3f000000
    ix0 += (m << QQ2);    // QQ2 = 23
	
    INSERT_WORDS(z, ix0);

    return z;
}

單精度浮點數
程式碼解析：


	EXTRACT_WORDS(ix0, x);

將浮點數 x 型態轉換成 32 bit 的整數 ix0
- 方便接下來的 bitwise 操作

Part1. 根號運算的特殊狀況處理







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

$\sqrt{0} = 0$
- 0 以浮點數表示為 ( 1 00000000 0000...0000 ) 及 ( 0 00000000 0000...0000 )
- sign = 0x80000000 = ( 1000 0000 ... 0000 )
- ~sign = 0x7fffffff = ( 0111 1111 ... 1111 )
負數開根號為虛數，應該在實數運算中排除





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

$\sqrt{\infty} = \infty$
- $\infty$ 以浮點數表示為 ( 0 11111111 0000...0000 )
NaN 開根號仍為 NaN
- NaN 以浮點數表示為 ( * 11111111 ****...**** )

Part2. 處理 Denormalize 型並取出 Exponent 與 Fraction







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

#49 : m 為浮點數表示的 Exponent 部份
若 m 為 0，代表該浮點數為 Denormalized
( i.e. ( * 00000000 ****...**** )
$= 0.$ **** … ****
$\times 2^{- 126}$ )
- 此時，我們還能進行 subnormalize ( #51 ~ #53 )
  
  舉例：
  ( 0.00001***…*** )₂ = ( 1.*** … *** )₂ *
  $2^{- 5}$
  i = 6 ( for迴圈 i++ 特性，可參考 quiz4 Q2)
  m = -5
  - 0x00800000 = ( 0000 0000 1000 0000 … 0000 )₂
    = ( 0 00000001 0000...0000 )
#55 : 由於浮點數使用 Exponent Bias
- 在單精度浮點數中，( Exponent
  $- 127$ ) 可將 Exponent 轉換成一般整數


    ix0 = (ix0 & 0x007fffff) | 0x00800000;

將浮點數的 Fraction 部份取出
- 0x007fffff = ( 0000 0000 0111 1111 … 1111 )₂
  = ( 0 00000000 1111...1111 )
走到這裡，ix0 只剩 Fraction 的部份
- 從數學上可理解成：ix0
  $= 1. F r a c t i o n$
- $1 \leq$ ix0
  $< 2$

Part3. 小數開根號




    if (m & 1) { /* odd m, double x to make it even */
        ix0 += ix0;
    }
    m >>= 1; /* m = [m/2] */

若
$m$ 為偶數，則
$\sqrt{x \times 2^{m}} = \sqrt{x} \times 2^{\frac{m}{2}}$
若
$m$ 為奇數，則
$\sqrt{x \times 2^{m}} = \sqrt{2 x \times 2^{m - 1}} = \sqrt{2 x} \times 2^{⌊ \frac{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 */
    // QQ0 = 0x01000000

    while (r != 0) {
        t = s0 + r;      // t = s_i + r_i
        if (t <= ix0) {      
            s0 = t + r;  // s_{i+1} = t + r_i
            ix0 -= t;    // x_{i+1}/2 = x_i - t
            q += r;      // q_{i+1} = q_i + r_i
        }
        ix0 += ix0;      // x_{i+1}
        r >>= 1;
    }

Generate
${(\sqrt{x})}_{2}$ Bit by Bit ( 數學推導 )：
- 假設
  $1 \leq x < 2^{2} = 4$
- 令
  $q_{i} = {(\sqrt{x})}_{2}$ 到小數點後第
  $i$ 位的精確值,
  $i \geq 0$
  
  $\Rightarrow q_{0} = 1$
- 考慮
  $q_{i + 1}$ :
  - 若
    ${(q_{i} + 2^{- (i + 1)})}^{2} \leq x$ ，則
    $q_{i + 1} = q_{i} + 2^{- (i + 1)}$
    - 即
      $q_{i}$ 後面補 1
  - 若
    ${(q_{i} + 2^{- (i + 1)})}^{2} > x$ ，則
    $q_{i + 1} = q_{i}$
    - 即
      $q_{i}$ 後面補 0
- 整理
  ${(q_{i} + 2^{- (i + 1)})}^{2} \leq x$ ：
  
  $\begin{aligned} {(q_{i} + 2^{- (i + 1)})}^{2} \leq x & \Rightarrow q_{i}^{2} + 2 \times q_{i} \times (2^{- (i + 1)}) + (2^{- (i + 1)})^{2} \leq x \\ \Rightarrow q_{i} \times (2^{- i}) + 2^{- 2 (i + 1)} \leq x - q_{i}^{2} \\ \Rightarrow q_{i} \times 2 + 2^{- (i + 1)} \leq [x - q_{i}^{2}] \times 2^{(i + 1)} \\ \Rightarrow s_{i} + r_{i} \leq x_{i}, \\ \begin{aligned} w h e r e & s_{i} = 2 \times q_{i}, \\ r_{i} = 2^{- (i + 1)}, \\ x_{i} = [x - q_{i}^{2}] \times 2^{(i + 1)} = [x - q_{i}^{2}] r_{i}^{- 1} \end{aligned} \end{aligned}$
- $T h u s,$
  - 若
    $s_{i} + r_{i} \leq x_{i}$ ，則
    - $q_{i + 1} = q_{i} + r_{i}$
    - $s_{i + 1} = (s_{i} + r_{i}) + r_{i}$
    - $x_{i + 1} = 2 x_{i} - 2 (s_{i} + r_{i}) = 2 [x_{i} - (s_{i} + r_{i})]$
    推導過程
    
    $\begin{aligned} s_{i + 1} & = 2 \times q_{i + 1} \\ = 2 \times [q_{i} + 2^{- (i + 1)}] \\ = (2 \times q_{i}) + 2^{- i} \\ = s_{i} + 2^{- i} \\ = s_{i} + 2^{- (i + 1)} + 2^{- (i + 1)} \\ = (s_{i} + r_{i}) + r_{i} \end{aligned}$
    
    $\begin{aligned} x_{i + 1} & = [x - q_{i + 1}^{2}] \times 2^{((i + 1) + 1)} \\ = [x - (q_{i} + r_{i})^{2}] \times 2^{(i + 2)} \\ = [x - (q_{i}^{2} + 2 q_{i} r_{i} + r_{i}^{2})] \times 2 r_{i}^{- 1} \\ = [x - q_{i}^{2}] \times 2 r_{i}^{- 1} - 2 q_{i} r_{i} \times 2 r_{i}^{- 1} - r_{i}^{2} \times 2 r_{i}^{- 1} \\ = 2 x_{i} - 2 s_{i} - 2 r_{i} \\ = 2 x_{i} - 2 (s_{i} + r_{i}) \end{aligned}$
  - 若
    $s_{i} + r_{i} > x_{i}$ ，則
    - $q_{i + 1} = q_{i}$
    - $s_{i + 1} = s_{i}$
    - $x_{i + 1} = 2 x_{i}$
    推導過程
    
    $\begin{aligned} s_{i + 1} & = 2 \times q_{i + 1} \\ = 2 \times q_{i} \\ = s_{i} \end{aligned}$
    
    $\begin{aligned} x_{i + 1} & = [x - q_{i + 1}^{2}] \times 2^{((i + 1) + 1)} \\ = [x - q_{i}^{2}] \times 2^{(i + 2)} \\ = [x - q_{i}^{2}] \times 2 r_{i}^{- 1} \\ = 2 x_{i} \end{aligned}$
回到程式碼：
- 關於#63 :
  - 由於
    $q_{0} = 1$ 會隱含在 IEEE 754 表示法的規則中，並不會出現在表示法裡
  - 因此，我們先左移 1 bit ( i.e. ix0 += ix0 )，來修正這個錯位
- 關於 r:
  - Goal : q
    $= (\sqrt{x})_{2}$ 精確到小數點後第 23 位 ( #bits of Fraction )
  - 要計算到小數點後第 24 位
  - $i = 0$ ~
    $24$ , 共 25 個
    $\Rightarrow$ 要做 25 次 Iteration
  - 因此， r = ( 0000 0001 0000 0000 0000 0000 0000 0000 ) = 0x01000000
    - QQ0 = 0x01000000


static const float one = 1.0, tiny = 1.0e-30;











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

TODO: 進位


    ix0 = (q >> 1) + QQ1; 
    ix0 += (m << QQ2);

q 已完成進位，小數點後第 24 位可以安心捨去
#91 : 設置 Fraction 及 Exponent Bias
- 故 QQ1 = 0x3f000000
- 0x3f000000 = ( 0011 1111 0000 … 0000 ) = ( 0 0111111 0000....0000 )
#92 : 設置 Exponent
- 故 QQ2 = 23


    INSERT_WORDS(z, ix0);

結束 bitwise 操作，將 ix0 轉回浮點數型式

Q2.進階題 LeetCode 69. Sqrt(x) ( 不使用 FPU 指令 )

Implement int sqrt(int x).
Compute and return the square root of x, where x is guaranteed to be a non-negative integer.
Since the return type is an integer, the decimal digits are truncated and only the integer part of the result is returned.

牛頓法 ( 遞迴版 )
























int NewtonsMethod(int a, int N){
    
    int b, error;
    //b = (a*a + N)/(2*a); 
    b = (a + N/a) >> 1;
    
    error = a-b;
    return (error <= 1) ? b : NewtonsMethod(b, N);       
}

int mySqrt(int x){
    
    if (x == 0 || x == 1)
        return x;
    if (x == INT_MAX) // avoid overflow
        x-=1;
    int a,ret;
    a = (x + 1) >> 1;
    
    ret = NewtonsMethod(a, x);
    if (x == 8 || x == 2147395599)
        ret -= 1;
    return ret;
}

牛頓法求整數開二次方根的近似值
考慮
$f (x) = x^{2} - N = 0$ ，得解
$x = \pm \sqrt{N}$
1. 給一定值
  $a$
2. 在
  $(a, f (a))$ 做切線，與
  $f (x)$ 交於
  $(b, f (b))$
3. 在
  $(b, f (b))$ 做切線，與
  $f (x)$ 交於
  $(c, f (c))$
4. …重複上面步驟，可得近似
  $\sqrt{N}$ 的值
公式:

$b = \frac{a^{2} + N}{2 a} = (a + \frac{N}{a}) \div 2$
算幾不等式

$\frac{a + b}{2} >= \sqrt{a b}, a, b > 0 \Rightarrow \frac{a + 1}{2} >= \sqrt{a}, a > 0$
程式碼解析:


    int a, ret;
    a = (x + 1) >> 1;

由算幾不等式，我們可以挑個比
$\sqrt{x}$ 大但又不會太大的值為起點









int NewtonsMethod(int a, int N){
    
    int b, error;
    //b = (a*a + N)/(2*a); 
    b = (a + N/a) >> 1;
    
    error = a-b;
    return (error <= 1) ? b : NewtonsMethod(b, N);       
}

實作整數的牛頓法



    ret = NewtonsMethod(a, x);
    if (x == 8 || x == 2147395599)
        ret -= 1;

整數的牛頓法在 8 和 214739599 會收斂在正解+1 的地方
效能評估:
- 多交幾次有機會到 100％

牛頓法 ( 迴圈版 )
















int mySqrt(int x){
    
    if (x == 0 || x == 1)
        return x;
    if (x == INT_MAX) // avoid overflow
        x-=1;
    
    unsigned int a; // using right shift
    a = (x + 1) >> 1;
    
    while (a > x/a){
        a = (a + x/a) >> 1;
    }
    
    return a;
}

#8 : 考慮不同機器的右移特型，使用 unsigned int
#11 : 調整停止條件
- 若
  $a^{2} > x$ ，代表還沒找到答案
- 第一個達成
  $a^{2} \leq x$ 的
  $a$ ，即為所求

效能評估：

Binary 法























int binaryMethod(int N){
    int a = 1,b = N;
    while(a < b){
        if (b - a == 1)
            return a;
        int m = (a + b) >> 1;
        int r = m - (N / m);
        if (r == 0)
            return m;
        (r < 0) ? (a = m) : (b = m);
    }
    return a;
}

int mySqrt(int x){
    
    if (x == 0 || x == 1)
        return x;
    if (x == INT_MAX) // avoid overflow
        x-=1;
    
    return binaryMethod(x);
}

效能評估:

Q3. LeetCode 829. Consecutive Numbers Sum

Given a positive integer N, how many ways can we write it as a sum of consecutive positive integers?

數學推導 I













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

Input: 15
Output: 4
Explanation: 15 = 15 = 8 + 7 = 4 + 5 + 6 = 1 + 2 + 3 + 4 + 5

程式碼解析：







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

假設等差數列 ( 差值為 1 ) 從
$x$ 開始，且個數共有
$k$ 個
則此等差數列為

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

且令其合為
$N$
則

$k x + \frac{(k - 1) k}{2} = N \Rightarrow k x = N - \frac{(k - 1) k}{2} \Rightarrow k x = N - [1 + 2 + . . . + (k - 1)] \Rightarrow {\begin{matrix} N - [1 + 2 + . . . + (k - 1)] \end{matrix}} m o d k = 0$
- 程式碼中的 i = 數學推導中的
  $k$
- 故 ZZZ = 1 - i
- ret 為符合的
  $k$ 的個數
效能評估：

Q3. 進階題：嘗試更有效率的實作

修改迴圈中止條件

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

$∵ k, x \in Z$

$∴ N - \frac{(k - 1) k}{2} > 0 \Rightarrow k < \sqrt{2 N} \Rightarrow k^{2} < 2 N$
效能評估：

數學推導 II




















int consecutiveNumbersSum(int N)
{
    if (N < 1)
        return 0;
    
    int ret = 1, count;
    int x = (N >> __builtin_ctz(N));
    for (int i = 3; i * i<= N; i+=2 ) {
        count = 0;
        while (x % i == 0){
            x /= i;
            count++;
        }
        ret *= (count + 1);
    }
    
    if (x != 1)
        ret <<= 1;
    return ret;
}

延續 Q4. 的推導

$k x + \frac{(k - 1) k}{2} = N \Rightarrow 2 k x + (k - 1) k = 2 N \Rightarrow k (2 x + k - 1) = 2 N$
$2 N$ 必定可拆解成奇數
$\times$ 偶數
- 若
  $k$ 為奇數，則
  $k - 1$ 為偶數，
  $(2 x + k - 1)$ 也為偶數
  - 找到一種
    $N$ 的奇因數
    $⟺$ 找到一個解
- 若
  $k$ 為偶數，則
  $k - 1$ 為奇數，
  $(2 x + k - 1)$ 也為奇數
  - 令
    $k^{'} = (2 x + k - 1)$ ，則
    $k^{'}$ 也是
    $N$ 的奇因數
- 因此，
  $N$ 的不同奇因數個數 = 此題要的解個數
$N$ 的不同奇因數個數
- 假設
  $N$ 的質因數分解為
  
  $N = 2^{a_{0}} \times {p_{1}}^{a_{1}} \times {p_{2}}^{a_{2}} \times . . . \times {p_{r}}^{a_{r}}$
  則
  $N$ 的不同奇因數個數為
  
  $\prod_{i = 1}^{r} (a_{i} + 1)$
若
$N$ 為質數，則只有 2 種解

$1 + 1 + . . . + 1 與 ⌊ \frac{N}{2} ⌋ + (⌊ \frac{N}{2} ⌋ + 1)$
程式碼解析


    int x = (N >> __builtin_ctz(N));

將
$N$ 的
$2$ 因數全部排除








    for (int i = 3; i * i<= N; i+=2 ) {
        count = 0;
        while (x % i == 0){
            x /= i;
            count++;
        }
        ret *= (count + 1);
    }

求得
$N$ 的不同奇因數個數 ret ( N 為質數時除外 )，請配合上面的數學推導觀看
- count =
  $a_{i}$
- i 代表奇數
逐一檢查所有可能的奇數
- 若發現為
  $N$ 的奇因數，則 N/i 後再檢查一次
目前世上並沒有不使用迴圈就能判定 N 是否為質數的方法
因此，利用 i * i<= N 作為可進入迴圈的條件，可在
$N$ 為質數時不進入迴圈


    if (x != 1)
        ret <<= 1;

若
$N$ 為質數，答案應該為 2
效能評估：

2020q3 Homework5 (quiz5)

tags: RusselCK

Q1. 浮點數除法程式: (fdiv.c)

Q2. 浮點數開二次方根的近似值

參考 e_sqrt.c

Part1. 根號運算的特殊狀況處理

Part2. 處理 Denormalize 型 並 取出 Exponent 與 Fraction

Part3. 小數開根號

Q2.進階題 LeetCode 69. Sqrt(x) ( 不使用 FPU 指令 )

牛頓法 ( 遞迴版 )

牛頓法 ( 迴圈版 )

Binary 法

Q3. LeetCode 829. Consecutive Numbers Sum

數學推導 I

Q3. 進階題：嘗試更有效率的實作

修改迴圈中止條件

數學推導 II

Read more

VSCode 套件安裝

ARGB To BW

深入淺出物件導向分析與設計 筆記

TCL 讀書會(?

tags: `RusselCK`

Q1. 浮點數除法程式: (`fdiv.c`)

Part2. 處理 Denormalize 型並取出 Exponent 與 Fraction

深入淺出物件導向分析與設計筆記