---
tags: Linux2020
---
# 2020q3 Homework5 (quiz5)
contributed by < `YLowy` >
## 測驗 1 fdiv.c
考慮到以下浮點數除法程式: (fdiv.c)
```c=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 / 2, od ? (slots + 1) >> 1 : slots >> 1);
if (od)
result += divop(result, slots);
return result;
}
```
假設 `divop()` 的第二個參數必為大於 0 的整數，而且不超過 ` int ` 型態能表達的數值上界。
### 在欣賞程式碼前可以先看看本程式欲實作概念 :
當我想要進行除法時
以浮點數 $X$ 除以整數 $Y$
1. 在二進位中我們可以很簡單的觀察整數 $Y$ 有無 $2$ 的因數，觀察最後一位即可，且只要右移即可完成除以 2 操作。
2. 浮點數 $X$ 則可以透過將 $exp$ $-1$ 很快速的將其除 2 。
3. 如果除數為偶數，也就是最後一個位元 $=0$ ，透過將除數與被除數都一起除以 $2$ 可以得到相同結果。
e.g. div(8,4) = div(4,2)
4. 最核心問題為如何處理除數為奇數的情況?
a) 這邊使用的做法為先將 除數 $+1$ 也就是先得到 $X ÷(Y+1)$ 這個數字，$Y+1$ 為偶數，我們可以用上述提到的方法。
b) 我們希望的目標為 $X ÷ Y$，算得 $X ÷(Y+1)$ 數值之後，我們需要補上其誤差。我們可以在 $Y$ 為大於 $0$ 整數前提下確定 $X ÷(Y+1)$ $<$ $X ÷ Y$ ，也就是我們需要在 $X ÷(Y+1)$ 值加上誤差值。
c) 誤差值計算 :
>$(原)$ ${ - }$ $(修正 (除數 +1)$
$X ÷ Y$ ${ - }$ $X ÷(Y+1)$
$=$ $\dfrac{X}{Y}$ - $\dfrac{X}{Y+1}$
$=$ $\dfrac{X(Y+1)}{Y(Y+1)}$ - $\dfrac{XY}{(Y+1)Y}$
$=$ $\dfrac{X}{Y(Y+1)}$
我們可以發現 $\dfrac{X}{Y+1}$ 這個結果就是我們在 a) 步驟得到的值
也就是說誤差值就是原先到的 $\dfrac{result}{Y}$ ，$Y$ 為原先傳入函式的除數。
### 有了以上概念我們可以非常輕鬆的觀察程式碼 :
```c=
#include <stdio.h>
#include <stdlib.h>
double divop(double orig, int slots) {
if (slots == 1 || orig == 0)
return orig; /*遞迴終止條件，除數為 1 或者 被除數為 0*/
int od = slots & 1;/*除數為奇數或者偶數以決定之後是否加入誤差值*/
double result = divop(orig / 2, od ? (slots + 1) >> 1 : slots >> 1);
/*奇數偶數分別的處理情況，也是 2,3 小點想處理的事情*/
if (od)
result += divop(result, slots); /* c) 小點的誤差值計算*/
return result;
}
```
:::info
除數部分最後一定會有遞迴停止點 : 7 -> 4 -> 2 -> 1
我覺得挺有意思的地方是當處理 "回傳無限循環小數" 的情況，搭配 GDB 可以觀察到整體函數運算過程
:::
### 延伸問題
1. 編譯器最佳化的角度，推測上述程式碼是否可透過 tail call optimization 進行最佳化，搭配對應的實驗
2. 比照 [浮點數運算和定點數操作](https://hackmd.io/@NwaynhhKTK6joWHmoUxc7Q/H1SVhbETQ)，改寫上述程式碼，提出效率更高的實作，同樣避免使用到 FPU 指令 (可用 gcc -S 觀察對應的組合語言輸出)，並與使用到 FPU 的實作進行效能比較
## 測驗 2 $\sqrt{2}$
延續 從 [√2](https://hackmd.io/@sysprog/2020-quiz5) 的運算談浮點數，假設 float 為 32-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式，請補上對應數值，使其得以運作：
```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)
```
避開 dangling else 的巨集 [參考](https://stackoverflow.com/questions/1067226/c-multi-line-macro-do-while0-vs-scope-block)
```c=24
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);
```
sign 為判斷正負之 mask
`EXTRACT_WORDS(ix0, x);` 可以想成 ix0 有與輸入 x 相同的記憶體，但是是以 int 形態存在。
``` c=34
/* 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;
}
```
處理 輸入為`無限大` `NAN` `負數` 的情況
```c=46
/* 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;
```
到目前為止可以確定輸入為 >0 的浮點數值
F : 23 bits
E : 8 bits
第`46`行 : m 會紀錄 E 的值
第`47`~`51`行 : 如果 m = 0 的情況，則為處理 Denormalize floating number 的操作 可以透過 bitwise 移動取得該值
第`53`行 : [IEEE754](https://en.wikipedia.org/wiki/IEEE_754) 所定義該浮點數的指數值
```graphviz
digraph D {
node_B [shape=record label="{0}|{E}|{E}|{E}|{E}|{E}|{E}|{E}|{E}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{...}|{...}|{F}|{F}"];
}
```
第`54`行 : 處理 Denormalize floating number 其格式，透過在 m 存下原本 E 值，ix0 存下 normalize 表示式，透過這個方法讓每個值都可以已 normalize floating number 格式處理。
也就是說原本Denormalize ${ 0.F.... * 2^E}$ 該 $F$ 的部分最高為會被推向 $E$ ，而在 ix0 中的 $E$ 會以 00...01 表示，但是真實的值會放在 m 中
```graphviz
digraph D {
node_B [shape=record label="{0}|{0}|{0}|{0}|{0}|{0}|{0}|{0}|{1}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{F}|{...}|{...}|{F}|{F}"];
}
```
我們可以表示該數字為 ${(-1)^0 * 1.F.... * 2^m}$
第`59`行 :而我們目標為進行數字 N 的方均根，以指數的角度 : ${\sqrt[2]{2^m}}$ = ${2^{m/2}}$
並會在第`88`行位移後寫回 ix0 EE....EE中。
```c=55
if (m & 1) { /* odd m, double x to make it even */
ix0 += ix0;
}
m >>= 1; /* m = [m/2] */
```
之後操作會在 m 格式為偶數情況下操作，所以如果 m 非偶數就和本身的 ix0 借個 2。
為何這麼做?
我們可以很簡單的計算 ${ \sqrt{100} = 10}$ ,${ \sqrt{10000} = 100}$，在這裡 m 可以想成 0的數量。
```c=60
/* generate sqrt(x) bit by bit */
ix0 += ix0;
q = s0 = 0; /* [q] = sqrt(x) */
r = 0x01000000; /* 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;
}
```
這邊得思考剩餘部分如何處理
透過牛頓逼近法
${ \sqrt[n]{N} ≒ \dfrac{(n-1)*a+\dfrac{N}{a^{n-1}}}{n} = a_{1}}$
可以透過多次逼近取得更接近的值
n=2 的情況時 :
${\sqrt{N} ≒ \dfrac{(2-1)*a+\dfrac{N}{a^{2-1}}}{2} = \dfrac{a+\dfrac{N}{a}}{2}}$
:::info
[快速反均方根](http://www.hsfzxjy.site/2016-08-03-uncover-the-secret-of-fast-inverse-square-root-algorithm/) 雷神之槌 中程式碼也是有利用到這點。
```c=
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
```
:::
#### 開均方根直式
這裡用到開均方根直式，以下以 625(1001110001) 方均根後得值 25(11001) 為例子
```
n bits 的數字 N ，開根號後只會為 n/2 bits
從最右邊開始左右兩個兩個一組
---------------
√ 10|01|11|00|01
在 X Y 處放進相同數值讓 X位置的數字 * Y那坨數字 為接近且小於的被除數
X
---------------
√ 10|01|11|00|01 (XY選擇放入 1 )
Y -> 1←
---------------
√ 10|01|11|00|01
1← 1 1
----------------
1
第一輪做完之後繼續第二輪，一直做到最後一位
1 X
---------------
√ 10|01|11|00|01
1 1
----------------
1 01 (XY選擇放入 1)
-> 1 1←
10Y ---------------
√ 10|01|11|00|01
1 1
原本Y * 2 ----------------
↘ 1 01
101← 1 01
-----------------
↙ 0
1 1 X
---------------
√ 10|01|11|00|01
1 1
----------------
1 01
101 1 01
----------------
0 11
110Y
(注意只有最後的 Y *2)
(也就是上輪10"1" 這個 1*2)
↘
1 1 0
---------------
√ 10|01|11|00|01
1 1
----------------
1 01
101 1 01
----------------
11
1100 00
----------------
11
↙
1 1 0 X
---------------
√ 10|01|11|00|01
1 1
----------------
1 01
101 1 01
----------------
11
1100 00
----------------
11 00
1100Y ↘
1 1 0 0
---------------
√ 10|01|11|00|01
1 1
----------------
1 01
101 1 01
----------------
11
1100 00
----------------
11 00
11000 00 00
----------------
↙ 11 00
1 1 0 0 X
---------------
√ 10|01|11|00|01
1 1
----------------
1 01
101 1 01
----------------
11
1100 00
----------------
11 00
11000 00 00
----------------
11 00
11000Y
↘ 1 1 0 0 1
---------------
√ 10|01|11|00|01
1 1
----------------
1 01
101 1 01
----------------
11
1100 00
----------------
11 00
11000 00 00
----------------
11 00 01
110001 11 00 01
----------------
0
最上面的 1 1 0 0 1 便是 √ 1001110001 的結果
也是就是 25
```
熟悉以上流程回到程式碼的部分 :
1. **最上面的 q 為所得之方均根答案**
ix0 代表被方均根數，逐步扣除以逼近 0
程式碼第70行 q += r; 以及
第63行 `r = 0x01000000;`
第73行 `r >>= 1;`
使其位移累加完成方均根之答案
2. **左側 藍色框框如下為程式碼中第 68 行 `s0 = t + r;` 之行為，r會再每次判斷後向右為位移，t 為當輪欲扣除的數。**
3. **第 67 行中的 if 判斷為決定能不能進行對 ix0 能不能扣除 X * Y那坨數字(也就是 t )，而 X 只可能為 0 或者 1 ，所以第 69 行程式碼 `ix0 -= t;` 讓 ix0 扣除逼近。**
4. **程式碼中的第62行 第72行 都寫有到 `ix0 += ix0;` 讓整個被方均根數直接向左位移一位，透過這個方法可以不用移動 t 以及 s0。**
```c=75
/* 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) + 0x3f000000;
ix0 += (m << 23);
INSERT_WORDS(z, ix0);
return z;
}
```
1. 最後一位小數會進行捨入，方法為透過 `z = one - tiny`，
`static const float one = 1.0, tiny = 1.0e-30;`,`z > one`判斷 rounding
2. float 可表達的範圍約3.4e-38~e.4e-38
這裡這樣實作原因我不曉得，不過對於操作 rounding 方法大致就 [這些](https://en.wikipedia.org/wiki/IEEE_754#Rounding_rules)
選擇自己希望的 rounding rule 即可
3. 最後把算得的數字拼回去，並且打包成 float 型態就可以回傳了。
### 延伸問題
1. 解釋上述程式碼何以運作
>搭配研讀 [以牛頓法求整數開二次方根的近似值](http://science.hsjh.chc.edu.tw/upload_works/106/44f6ce7960aee6ef868ae3896ced46eb.pdf) (第 19 頁)
2. 嘗試改寫為更有效率的實作，並與 FPU 版本進行效能和 precision /accuracy 的比較
3. 練習 [LeetCode 69. Sqrt(x)](https://leetcode.com/problems/sqrtx/)，提交不需要 FPU 指令的實作
利用牛頓逼近法計算輸出 nhat 的每一個 bit 是 0 或者 1。
```c=
int mySqrt(int v)
{
unsigned long temp, nHat = 0, b = 0x8000, bshft = 15;
do{
if (v >= (temp = (((nHat << 1) + b) << bshft--)))
{
nHat += b;
v -= temp;
}
} while (b >>= 1);
return nHat;
}
```
參考資料 : [快速整数平方根算法](https://andy1314chen.github.io/2017/06/28/FastSqrt/)
## 測驗 3 Consecutive Numbers Sum
[LeetCode 829. Consecutive Numbers Sum ](https://leetcode.com/problems/consecutive-numbers-sum/)
先參考 [這裡](https://hackmd.io/@sysprog/2020-quiz5#%E6%B8%AC%E9%A9%97-3)
```c=
int consecutiveNumbersSum(int N)
{
if (N < 1)
return 0;
int ret = 1;
int x = N;
for (int i = 2; i < x; i++) {
x += 1-i;
if (x % i == 0)
ret++;
}
return ret;
}
```
::: info
看起來超簡單
上課有提到跟著題目走就可以知道怎麼做，可是看到程式的我覺得
:::
### 程式概念
直接從例子做起，以 N = 9 為例子 :
第一個最容易看到的為數量為 1 ，也就是 9
觀察表格以及其對應公式 (n>0,n位任意數) ，而我們只需要考慮 0~N 的範圍即可。
> e.g. 若 k=2 ，連續兩個數字之和我們可以發現為 1+2,2+3,3+4,4+5,...
> 也就是表格中的 k=2 該列 : 3,5,7,9,....
```graphviz
digraph D {
node_B [shape=record label="{連續數字的數量(k)|1|2|3|4|5|6|...}|{{對於該數量 k 可能 N 的值}|{1|2|3|4|5|6|7|8|9}|{03|05|07|09|11|13|15|17}|{06|09|12|15|18|21|24}|{10|14|18|22|26|30}|{15|20|25|30|35}|{21|27|33|39}|{....}}|{公式: func(n) =|0+1*n|1+2*n|3+3*n|6+4*n|10+5*n|15+6*n|{...}}"];
}
```
我們只要進行 for 迴圈對 "k" 檢查，也就是(連續之數量)
而判斷是否存在存在連續數量之正整數只需要觀察公式是否符合，
程式碼實作 :
```c=
int consecutiveNumbersSum(int N)
{
if (N < 1)
return 0; /*若是輸入0 回傳0*/
int ret = 1; /*每一個 N 都存在 1*N，所以保存 1 個結果*/
int x = N;
for (int i = 2; i < x; i++) {
x += 1-i; /* 這行是為了應對公式的 "XX" + k*n */
/* 這邊比較特別會在下面解釋 */
if (x % i == 0) /* i可以觀察表格公式該列 */
ret++; /* k 個連續整數加總應要符 */
} /* 合 XX + k*n (這裡i=k) */
return ret;
}
```
對於 `x += 1-i;` 我們可以觀察迴圈在對 N 檢查每個 k 。
參考最右行所列舉，檢查 N 是否有連續數值和 : (粗字體為程式碼中的 -(i-1) )
>(N - **1**) % 2
>((N - 1) -**2**) % 3
>(((N - 1) -2) -**3**) % 3
>..... 以此類推
### 延伸問題
1.嘗試改寫為更有效率的實作
Sign in with Wallet
Connect another wallet