2020q3 Homework5 (quiz5)

contributed by < Rsysz >.

tags: `sysprog`

測驗題

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

將浮點數 orig 除以大於 0 的整數 slots
若 slots 為偶數

\frac{o r i g / 2}{s l o t s / 2} = \frac{o r i g}{s l o t s}

若 slots為奇數，透過 result += divop(result, slots) 補償

\frac{o r i g / 2}{(s l o t s + 1) / 2} + \frac{\frac{o r i g / 2}{(s l o t s + 1) / 2}}{s l o t s} = \frac{o r i g / 2}{(s l o t s + 1) / 2} \times (1 + \frac{1}{s l o t s}) = \frac{o r i g / 2}{(s l o t s + 1) / 2} \times \frac{(s l o t s + 1)}{s l o t s} = \frac{o r i g}{s l o t s}

對奇數 slots 加 1 後右移與補償，巧妙地加速計算速度，D1 = (c) 2, D2 = (d) 1

2.

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

首先看到 do while(0) 封裝的 marco，這種用法是為了確保替換後程式正確編譯與減少內部 Bug 的發生，再來看到各 marco 的功能。

INSERT_WORDS 將 int ix0 轉為 float d
EXTRACT_WORDS 將 float d 轉為 int ix0
因無論 unsigned 與否，內部數值皆為相同的，因此皆以 int 表示

接著看到整體程式碼，首先將 float x 轉為 int ix0

ix0 <= 0 也就是 float x 為負數時有兩個結果
- (ix0 & (~sign)) == 0 表示除 sign 以外皆是 0，代表 float x為 +-0 開根號還是 0，回傳 x
- 先前條件不滿足且 ix0 < 0，表示 float x 為負浮點數，無法開根號，因此回傳 NaN

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

接著處理 INF 與 NaN，先前條件不滿足且 (ix0 & 0x7f800000) == 0x7f800000 表示 float x 的 Exponent 部分皆為 1

+INF 表示為
$1_{s i g n} 11111111_{e x p o n e n t} . . ._{f r a c t i o n}$ fraction 皆為 0
NaN 表示為
$0_{s i g n} 11111111_{e x p o n e n t} . . ._{f r a c t i o n}$ fraction 部分皆可，此外先前 (ix0 < 0) 已排除 sign 為 1 的結果

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

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

對浮點數做正規化，首先取出 ix0 的前九個bit 也就是 float x 的 sign 與 exponent 部分
- 若 m == 0 則對 ix0 做左移，直到 exponent 部分的
  $2^{0}$ 為 1
- 因 m 表示 exponent 部分，根據上圖對浮點數的正規化，m = m + 1 - i
接著將 m - 127 得到對應 exponent 的數字，並以 ix0 = (ix0 & 0x007fffff) | 0x00800000 擷取 ix0 的後 23 bit並將第 24 bit 設為 1，得到
$1. f r a c t i o n$
最終得到
$2^{m} \times i x 0 ， 1 \leq i x 0 < 2$ 接著為了將
$\sqrt{2^{m}} = 2^{\frac{m}{2}}$ ，若 m 為奇數，將多的 2 乘至 ix0 部分，再將 m 除以 2 得到
$2^{\frac{m}{2}} \times (2 \times 1. f r a c t i o n) ， 1. f r a c t i o n \leq i x 0 \leq (2 \times 1. f r a c t i o n)$
又
$1 \leq 1. f r a c t i o n < 2 ，所以 1 \leq i x 0 < 4$ ，整理如下
$2^{\frac{m}{2}} \times i x 0 ， 1 \leq i x 0 < 4 ， m \in e v e n$

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

接著因
$1 \leq i x 0 < 4$ 表示
$1 \leq \sqrt{i x 0} < 2$ ，使得可以從 1 開始逼近
$\sqrt{i x 0}$ 的數值。
此外前面將 fraction 部分擴大至 24 bit(新增 1 的部分)，這裡又將 ix0 += ix0，因此範圍變成
$2 \leq i x 0 < 8$ ，所以為了保證 ix0 最小值 1.fraction 能正確逼近，這裡要從 25 位開始做判斷， QQ0 = (a) 0x01000000

想請問這裡一下，這裡 ix0 += ix0 與 ix0 << 1 在各種方面上有差異嗎?

令
$q_{i} = \sqrt{i x 0}$ 取小數點後
$i$ 位的精度
$q_{0} = 1$ ，
$s_{i} = 2 q_{i}$ 且
$i x 0_{i} = 2^{i + 1} (i x 0 - q_{i}^{2})$
- 為了透過
  $q_{i}$ 計算
  $q_{i + 1}$ ，若
  $(q_{i} + 2^{- (i + 1)})^{2} \leq i x 0$ 成立則
  $q_{i + 1} = q i + 2^{- (i + 1)} (補 1) 否則 q_{i + 1} = q_{i} (補 0)$
- 接著透過代換得到
  $s_{i} + 2^{- (i + 1)} \leq i x 0_{i}$
  - 若成立，
    $s_{i + 1} = s_{i} + 2^{- i} ， i x 0_{i + 1} = i x 0_{i} - s_{i} - 2^{- (i + 1)}$
  - 否則，
    $s_{i + 1} = s_{i} ， i x 0_{i + 1} = i x 0_{i}$
這裡以圖解說明逼近過程，首先假設 ix0 = 1.10...

ix0 += ix0

q = s0 = 0
r = 0x01000000
t = s0 + r 也就是

s_{i} + 2^{- (i + 1)} \leq i x 0_{i}

首項 s 為 0 是因為還沒進入小數位
t <= ix0 t = 0x01000000 < ix0 = 0x018XXXXX
s0 = t + r 也就是

s_{i + 1} = s_{i} + 2^{- i}

s0 = 0x02000000
ix0 -= t 也就是

i x 0_{i + 1} = i x 0_{i} - s_{i} - 2^{- (i + 1)}

ix0 = 0x008XXXXX

q += r q = 0x01000000
ix0 += ix0 ix0 = 0x01XXXXXX

r >>= 1 r = 0x00800000

透過 ix0 += ix0 對 ix0 左移，透過 r >> 1 對 r 右移，不斷逼近
$\sqrt{i x 0}$ 並將逼近後每位數值存於 q
剩下的數值將會根據不同系統上的精度做 rounding，四捨五入

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

最後將計算結果還原成浮點數表示
- q 表示的是小數點落在 24~25 bit間的編碼，因此右移一次回到正常範圍
- QQ1 = (d) 0x3f000000，照理來講應該加上 0x3f800000 ，但因 q = 1.fraction 其最高位的 1 右移時補上了缺口
m 表示 exponent 部分，QQ2 = (g) 23 左移回正確位置

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

    INSERT_WORDS(z, ix0);
    return z;

延伸問題

LeetCode 69. Sqrt(x)

typedef union {
    float value;
    uint32_t v_int;
} ieee_float_shape_type;

#define EXTRACT_WORDS(ix0, d)        \
    do {                             \
        ieee_float_shape_type ew_u;  \
        ew_u.value = (d);            \
        (ix0) = ew_u.v_int;          \
    } while (0)

int mySqrt(int x)
{
    int32_t sign = 0x80000000;
    int32_t ix0, m, 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 */
    }   
    /* 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] */
    
    unsigned int lo = 1 << m, hi = 1 << (m + 1);
    unsigned int mid = (lo + hi) >> 1;
    while(1){
        if(x > (mid + 1) * (mid + 1)){
            lo = mid;
            mid = (lo + hi) >> 1;
        }
        else if (x < mid * mid){
            hi = mid;
            mid = (lo + hi) >> 1;
        }
        else
            return mid;
    }
}

因本題無須深入求小數點位的精度要求，因此借用本題方式轉換為

2^{\frac{m}{2}} \times 1. f r a c t i o n

接著利用

2^{\frac{m}{2}}

做 Binary search 以此快速縮小範圍，最終回傳 mid。

3.

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

給定正整數N，求N能寫成多少種連續正整數的和，假設從 x 開始，且個數共有 k 個
根據等差數列的求和公式

k x = N - \frac{(k - 1) k}{2}

求有幾種組合k 與 x 能滿足 N - {0 ~ k}的和 = kx
- 首先看到 x = N 對應著公式中的 N，代表 ZZZ 對應著 {0~k}的和，因此 i 對應 k
- for 迴圈則是用來計算不同終點的等差數列，所以ZZZ = (d) 1 - i
舉例來說，若 N = 3
- 第一個迴圈N - {0~1}的和 = 2 也就是 x += 1 - i，因此檢查 x % i 也就是
  $x = \frac{N - \frac{(k - 1) k}{2}}{k}$ 就能知道是否有可行的解
- 第二個迴圈N - {0~2}的和 = 0 因 x 於第一個迴圈已被更新，因此算出結果仍然等於 x += 1 - i，且 x % i == 0 因此得到 x = 0 的整數解一個，以此類推，直到 k = x。

2020q3 Homework5 (quiz5)

tags: sysprog

1.

2.

延伸問題

3.

Read more

MTK / Bluetooth

ONVM-IPS / Hyperscan

VPP & DPDK

2020q3 Gameboy-emu

tags: `sysprog`