Try   HackMD

2024q1 Homework4 (quiz3+4)

contributed by < LindaTing0106 >

第三週測驗題

測驗 1

i_sqrt 版本一利用 log2 函式得出 N 的最高有效位元在 msb 。 隨後再迴圈內檢查該 N 為二進位表示時,將在逼近的數字加上

2a 後將其平方後,不會超過原數字。
藉由此種方式得出 N 的平方根。

  • 版本一
#include <math.h>
int i_sqrt(int N)
{
    int msb = (int) log2(N);
    int a = 1 << msb;
    int result = 0;
    while (a != 0) {
        if ((result + a) * (result + a) <= N)
            result += a;
        a >>= 1;
    }
    return result;
}
  • 版本二

將 log2 函式利用迴圈方式將其取代。

int i_sqrt(int N)
{
    int msb = 0;
    int n = N;
    while (n > 1) {
        n >>= 1;
        msb++;
    }
    int a = 1 << msb;
    int result = 0;
    while (a != 0) {
        if ((result + a) * (result + a) <= N)
            result += a;
        a >>= 1;
    }
    return result;
}

i_sqrt 版本三則是利用 Digit-by-digit calculation 的方式,其中先將要開平方根的數

N 用二進位方式表示。

假設

N=(25)2,使其變為
N=(11001)2

若將其表示位十進位則得
N=(24+23+21)2
,則要算出
x
就要將其展開,則為
(a4+a3+a2+a1+a0)2
,在此
ai=2i
or
0
,將其利用矩陣展開後可以得其規律。

an2+(an1+2an)an1+(an2+2(an1+an))an2+...+(2i=1nai+a0)a0
故此處可以假設
Pm=an+an1+...+am

an2+(an1+2Pn)an1+(an2+2Pn1)an2+...+(2P1+a0)a0

Pm=Pm+1+am
Pm
N
的平方根,故我們想求出
Pm
的話就需要去知道每個
am
為零或不為零,但如果每一輪都將
NPm2
的話將會導致運算成本過高。

  • Xm=NPm2=Xm1Ym
  • Xm1=NPm+12
  • Ym=Pm2Pm+12=Pm+12+2Pm+1am+am2Pm+12=2Pm+1am+am2

在這裡將

Ym 拆解為
cm
dm
,其中

  • cm=2Pm+1am=Pm+12m+1
  • dm=am2=(2m)2
  • Ym={cm+dm  ,if am=2m0      ,if am=0

以此類推

  • cm1=2Pmam1=Pm2m=(Pm+1+am)2m=Pm+12m+am2m={cm/2+dm ,if am=2mcm/2    ,if am=0
  • dm1=am12=(2m1)2=(2m)222=dm4

則在程式進行時,會由

2n 往下查看
N
中有無
2i

  • Xn1=N

  • cn=0

  • P0=an+an1+...+a0=c1
    故可以看出,求出
    c1
    則等於求出
    N
    的平方根為何。

  • 版本三

int i_sqrt(int x)
{
    if (x <= 1) /* Assume x is always positive */
        return x;

    int z = 0;
    for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
        int b = z + m;
        z >>= 1;
        if (x >= b)
            x -= b, z += m;               
    }
    return z;
}

確認 x 不為 0 或是其他小於 0 的數。

 if (x <= 1) /* Assume x is always positive */
        return x;

int m 為

dm ,故第一次進來迴圈時, int m =
dn=(2n)2
。 因為
(2n)2=2n2
n2
必不可能為奇數,因此要 & ~1UL
因為
dm1=dm4
故每次迴圈往右位移 2 。

for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {

Ym={cm+dm  ,if am=2m0      ,if am=0

並且每次得出

cm1={cm/2+dm ,if am=2mcm/2    ,if am=0

且若

Xm>=Ym 則減去
Ym
Xm1

int b = z + m;
z >>= 1;
if (x >= b)
    x -= b, z += m;  

ffs 取代 __builtin_clz 實作

在每一次迴圈中,使用 ffs 得到 temp 中最低有效位元 var 後,將其右移 var ,一直重複此操作,最終可以得到其最高有效位元的位置。但還需要將其減 1 才會與 (31 - __builtin_clz(x)) 相等。

int msb = 0;
for (int temp = x; temp; ){
    int var = ffs(temp);
    msb += var;
    temp >>= var;
}

for (int m = 1UL << ((msb - 1) & ~1UL); m; m >>= 2) 

在迴圈內使用 mask 取代原先條件句。

int mask = -(x >= b);
x -= mask & b;
z += mask & m;

commit aa272b4

測驗 2

一般在計算商數和餘數時,我們都會直接使用 / 、 % 對其進行運算,如此造成運算成本較高的後果。

static void str_add(char *b, char *a, char *res, size_t size)
{
    int carry = 0;

    for (int i = 0; i < size; i++) {
        int tmp = (b[i] - '0') + (a[i] - '0') + carry;
        carry = tmp / 10;
        tmp = tmp % 10;
        res[i] = tmp + '0';
    }
}

因此採用 bitwise operation 來實作除法,但可以發現當我們希望將 n 除以 10 時,無法直接使用 shift 指令完成,因為 10 並不能直接使用

2i 的方式來表示。因此我們需要利用 n / x 來完成除以 10 的動作。
但我們應該要先證明其除以 x 後,應該在小數點後第一位時的精確度都是正確的。
故設

  • n=ab
  • l=cd
  • nl

ab10nxab.910 ,確保 n 除以 x 後其跟 n 除以 10 直到小數點第一位時的精確度都是正確的。

ab10nxab.910a.bnxa.b9na.b9xna.b

,由於

n=ab ,故
na.b9x10

cd10lxcd.910cd10la.b9ncd.910cd10cda.b9abcd.910

a.b9ab=a.b+0.09ab ,故

cda.b9ab=cd(a.bab+0.09ab)=cd110+cd0.09ab=c.d+cd0.09ab

cd10c.d+cd0.09ab , 而因為
cd<ab
,故
c.d+cd0.09abcd.910

故可以保證

n/x 其在小數點後第一位的精確度。

因為被除數的大小不會超過 19 ,所以我們接著直接用 19 來做討論。

故從上面式子,我們可以帶入

n=19

191019x19.910191.99x191.99.54x10

故我們只要找到一

x 其在此範圍內,便可保證將某數除以
x
時,其到小數點後第一位的精確度都會是正確的。

在這邊老師用了

128139.85 這個數字來當作 x 。

故將 n 除以

12813 ,即可取代除以 10 ,而 n 除以
12813
同時也等同於 n *
13128
。因此我們要想辦法用 bitwise operation 操作出
13128

((((tmp >> 3) + (tmp >> 1) + tmp)

得出

tmp8+tmp2+tmp=13tmp8

((((tmp >> 3) + (tmp >> 1) + tmp) << 3

往左位移 3 讓

13tmp88=13tmp

為了避免在向右位移時導致前面 3 個位元丟失,先用 mask 的方式將位元暫存下來。

d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;
q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2)

但這裡不太懂的是 d2 應該是為了避免 (tmp >> 3) 丟失的位元, d0 則是應該是為了避免 (tmp >> 1) 丟失的位元,那為什麼需要加 d1 呢?

接著再將

13tmp/128 ,故往右位移 7 。

q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7

最後為了要算出 mod 10 ,將 tmp 減 ( div 10 的結果乘以 10 )。

r = tmp - (((q << 2) + q) << 1);

最終可將程式包裝為

#include <stdint.h>
void divmod_10(uint32_t in, uint32_t *div, uint32_t *mod)
{
    uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */
    uint32_t q = (x >> 4) + x;
    x = q;
    q = (q >> 8) + x;
    q = (q >> 8) + x;
    q = (q >> 8) + x;
    q = (q >> 8) + x;

    *div = (q >> CCCC);
    *mod = in - ((q & ~0x7) + (*div << DDDD));   
}

在這邊我們會先把

in|1 ,之後再將其減掉
in>>2
,如此我們可以得到
x=in34

uint32_t x = (in | 1) - (in >> 2);

in | 1 是在 in 為偶數時,會使 in + 1 ,其用意為?

q=in364+in4864=in102128
1021280.8=810

uint32_t q = (x >> 4) + x;

此處用來使

102128 更接近
0.8

x = q;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;

如果不進行逼近的話,會導致結果有錯嗎?

8in10 變為
in10

*div = (q >> 3);

最後為了要算出 mod 10 ,將 in 減 ( div 10 的結果乘以 10 )。
由於 (q & ~0x7) = q & ~0..0111 = q & 1..1000 。也就等於將 *div << 3 。

故為要湊出 *div

10 。故將 *div << 3 + *div << 1 。

*mod = in - ((q & ~0x7) + (*div << 1));

測驗 3

ilog2 版本一是用來計算以 2 為底的對數,原先的版本為利用將二進位的數字往右位移 1 ,則表示其可以除以 2 一次,重複以上算式,直到數字被除至 0 。

int ilog2(int i)
{
    int log = -1;
    while (i) {
        i >>= 1;
        log++;
    }
    return log;
}

但此函式的缺點為,此數字的最高 set bit 的位置在 n 處,該迴圈就應該要執行 n 次。
故將其改寫為 ilog2 版本二,
若有一數為

232 ,則在原本的程式中需要執行 32 次,但在改寫完的程式中只須執行 2 次。因此可省下原本需要的時間。

while (i >= 2^16) {
    result += 16;
    i >>= 16;
}
while (i >= 2^8) {
    result += 8;
    i >>= 8;
}
while (i >= 2^4) {
    result += 4;
    i >>= 4;
}
while (i >= 2) {
    result += 1;
    i >>= 1;
}

亦可寫為 ilog2 版本三, builtin_clz 用以得出一二進位數字其最高有效位元前,共有多少個 0 ,在利用 31 減去此數後,則可以得到其為

2i 。而需要傳入 v | 1 則是為了避免 v 為 0 時,由於 builtin_clz 對於 0 未定義導致錯誤的情況發生。

int ilog32(uint32_t v)
{
    return (31 - __builtin_clz(v | 1));
}

測驗 4

指數加權移動平均為一種在取平均的手段,其特點為可以調整歷史資料加權降低的程度。其公式為

St={Y0,t=0αYt+(1α)St1,t>0 ,其中
α
表示歷史資料加權降低的程度,
Yt
表示在
t
時獲得的新資料,
St
表示在
t
時算出的指數加權移動平均。將以上式子展開後可得
St=αYt+(1α)αYt1+(1α)2αYt2+...+(1α)iαYti+...+(1α)tαY0

在程式中,利用 factor 去縮放新傳進來的資料, weight 則是表示歷史資料加權降低程度。

測驗 5

實作出 ceil_ilog2 ,在最一開始處先將 x - 1 ,此目的為最終可以直接無條件進位,下面 r = (x > 0xFFFF) << 4; 則是若

x>216 則將
1<<4
並傳給 r ,而後再讓 x 去向右位移 r ,最後是為了判斷 x 是否是一個大於 0 的數字,若 x 大於 0 則其 ceil_ilog2 必定要加一。

改進程式碼,使其得以處理 x = 0 的狀況,並仍是 branchless

利用 x 減去 (x > 0) ,使當 x = 0 的情況下不會減到 1 , x 也就不會變為 0xFFFF_FFFF ,最後在原本加 1 處改為加 (原本的 x >= 1) ,讓 x 為 0 時不會加 1 。

commit 1fd5a21


第四週測驗題

測驗 1

在 popcount_naive 函式中的缺點為,如果 v 中有多少個 set bit ,則程式就需要執行多少次。故引入了 popcount_branchless 的寫法,利用

xx2x4x8 可以知道在 x[3:0] 範圍內 1 的數量。
推導為:設
x=23b3+22b2+21b1+20b0
,其中
bi=0
or
1

xx2x4x8 =
23b3+22b2+21b1+20b022b3+21b2+20b121b3+20b220b3

整理後為

(23222120)b3+(222120)b2+(2120)b1+(20)b0

故可以發現在

b3 為 1 的時候
(23222120)b3
也會為 1 ,以此類推。

故可以得知確實可以利用

xx2x4x8 得知 4 bits 範圍內的 1 的 bits 數。

    unsigned n;
    n = (v >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;

利用此段程式碼我們可以得到此圖。

b_31 b_30 b_29 b_28 ... b7 b6 b5 b4 b3 b2 b1 b0  // v
   0 b_31 b_30 b_29 ...  0 b7 b6 b5  0 b3 b2 b1  // (v >> 1) & 0x77777777
   0    0 b_31 b_30 ...  0  0 b7 b6  0  0 b3 b2  // (n >> 1) & 0x77777777
   0    0    0 b_31 ...  0  0  0 b7  0  0  0 b3  // (n >> 1) & 0x77777777

也可以得到此算式

vv2v4v8

  • 假設 v = 32'b0100_1111_0000_0100_0111_0001_0001_1100
  • 則運算完的結果會為 v = 32'b0001_0100_0000_0001_0011_0001_0001_0010 ,但可以看出我們應該要將每個 nibble 的位置調整好,相加後才是 1 的總數。
v = (v + (v >> 4)) & 0x0F0F0F0F

假設

Bn 代表第 n 個 nibble (4 位元) 中的數值

B7 B6 B5 B4 B3 B2 B1 B0  // v
 0 B7 B6 B5 B4 B3 B2 B1  // (v >> 4)
B7 (B7+B6) (B6+B5) (B5+B4) (B4+B3) (B3+B2) (B2+B1) (B1+B0)  // (v + (v >> 4))
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)  // (v + (v >> 4)) & 0x0F0F0F0F
  • 則繼續使用剛才的例子 v = 32'b0001_0100_0000_0001_0011_0001_0001_0010
  • 經過此式子後v = 32'b0000_0101_0000_0001_0000_0100_0000_0011
v *= 0x01010101

此處的

Ai 為上面剛計算完的
Bi+1+Bi

                         0 A6  0 A4  0 A2  0 A0
                      x  0  1  0  1  0  1  0  1
---------------------------------------------------
                         0 A6  0 A4  0 A2  0 A0
                   0 A6  0 A4  0 A2  0 A0  0
             0 A6  0 A4  0 A2  0 A0  0
       0 A6  0 A4  0 A2  0 A0  0
---------------------------------------------------
                            ↑_______________________A6+A4+A2+A0
  • v = 32'b0000_0101_0000_0001_0000_0100_0000_0011
                         0 0101  0 0001  0 0100  0 0011
                      x  0    1  0    1  0    1  0    1
-------------------------------------------------------
                         0 0101  0 0001  0 0100  0 0011
                 0 0101  0 0001  0 0100  0 0011
         0 0101  0 0001  0 0100  0 0011
 0 0101  0 0001  0 0100  0 0011
-------------------------------------------------------
                            ↑_______________________每個 nibble 加總起來的地方
  • v = 32'b0000_1101_0000_xxxx_0000_xxxx_0000_xxxx

可以看到我們所求的數字就位於 v[31:24] 處,故最後將其向右為移即可結束。

return v >> 24
  • return 32'b..._0000_1101
unsigned popcount_branchless(unsigned v)
{
    unsigned n;
    n = (v >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;

    v = (v + (v >> 4)) & 0x0F0F0F0F;
    v *= 0x01010101;                                     

    return v >> 24;
}

在與算出 hammingDistance 有關的題目中,可以利用 popcount 計算兩數有幾個不相同的值,而最後回傳回去時需要右移 1 ,是因為此處會重複運算到一次,所以需要除以 2 也就是右移 1 。

int totalHammingDistance(int* nums, int numsSize)
{
    int total = 0;;
    for (int i = 0;i < numsSize;i++)
        for (int j = 0; j < numsSize;j++)
            total += __builtin_popcount(nums[i] ^ nums[j]); 
    return total >> 1;
}

使用位元方式找尋 Total Hamming Distance

用位元的位置去看,在 x 位元時, nums 中有幾個 1 ,從此處也可以得出 0 的數量,把 1 數量乘上 0 的數量,即是當前位元位置的 Hamming Distance ,最後把所有位置的 Hamming Distance 都加總起來後即是 Total Hamming Distance 。

commit 840a409

測驗 2

由於

2k={1(mod3),k even1(mod3),k odd , 故
n=bn1bn2...b2b1b0
,而
bi=1
,若 i 為奇數,則 n mod 3 的結果會 -1 。

n(mod3)=0n1bi(1)i(mod3)

可以看出我們只要將

bi 在偶數的位元數減去
bi
在奇數的位元數即可得到 n 。

故得出此式 n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)

tic-tac-toe 程式碼中,進行一個 3x3 大小的圈圈叉叉,由於大小為 3x3 故總操作數只會有 9 次。此時的操作為每次從陣列中選出一個可以操作的位置,選擇完後將本次選擇的位置,從陣列中移除。並將其給 board | move_masks[move] ,隨後去檢查該名玩家是否勝利。

static const uint32_t move_masks[9] = {
    0x40040040, 0x20004000, 0x10000404,
    0x04020000, 0x02002022, 0x01000200,
    0x00410001, 0x00201000, 0x00100110,
};

若左邊有一條直線達成連線則 player_board = 0x44470041
若中間有一條直線達成連線則 player_board = 0x22207022
若右邊有一條直線達成連線則 player_board = 0x11100714

若上面有一條橫線達成連線則 player_board = 0x70044444
若中間有一條橫線達成連線則 player_board = 0x07022222
若右邊有一條橫線達成連線則 player_board = 0x00711111

若左上有一條斜線達成連線則 player_board = 0x42142172
若右上有一條橫線達成連線則 player_board = 0x12412427

可以看出只要 player_board 在 16 進位的表示下,只要有其中一個位元為 7 則表示勝利。
故在此段程式碼中將每個位元都加 1 ,再將每個位元都 & 8 ,如果回傳不為 0 則表示此為玩家勝利。

static inline uint32_t is_win(uint32_t player_board)
{
    return (player_board + 0x11111111) & 0x88888888;
}

測驗 3

在 Makefile 中可以看到當執行 make check 命令時,其會去對 treeint 檔進行編譯,可以從其中看到 insert, find, remove 等操作的執行時間。而 treeint 呼叫的函式是在 treeint_xt 定義。

  • treeint_xt_destroy 函式中會呼叫 xtree 中的 xt_destroy
    • xt_root 存在時,會再進入 __xt_destroy,使用遞迴方式將樹的左右子樹刪除,最終再刪除自己本身。
  • treeint_xt_insert 函式中會呼叫 xtree 中的 xt_insert ,在函式中會呼叫 __xt_find,去尋找應該 insert 至何處,如果找到此值已經在樹中了則不用新增,反之找到該新增位置後,就將節點插入進去。
  • treeint_xt_remove 函式中會呼叫 xtree 中的 xt_remove , 在函式中會呼叫 xt_find ,但 xt_find 不同於 __xt_find, xt_find 將不會回傳位置,其目的只是在於確認欲刪除的節點確實存在。
    • 如果欲刪除的節點有右子樹:將要被刪除的節點用 xt_replace_right 替換為右子樹中第一個遇到的節點,也就是右子樹中最左邊的節點。
    • 如果欲刪除的節點有左子樹:將要被刪除的節點用 xt_replace_left 替換為左子樹中最後一個遇到的節點,也就是左子樹中最右邊的節點。
  • xt_replace_right
  • xt_replace_left