# 2020q3 Homework5 (quiz5)
contributed by < `Tim096` >
###### tags: `進階電腦系統應用2020` `quiz5`
> [浮點數運算 ＆ 數值系統 ＆ Bitwise 練習題](https://hackmd.io/@sysprog/2020-quiz5)
:::spoiler 期中心(ㄈㄟˋ)得(ㄏㄨㄚˋ)
文長慎入
課程也差不多過一半了，打一點小心得，其實在上第一堂課程完畢以後心裡有點被打擊到，想說我剛考上成大來或許也沒有那麼爛了吧，或許未來的路會更順利一些吧，並且事先打聽了許多的消息都說只要上了這一堂課你未來的人生一定可以更順利並且朝者自己夢想更前進一步，以為這一堂和其他課程一樣「輕鬆」只是價值比較高。
對，一開始就是抱持者這種不成熟的態度來進去聽完了第一堂課，第一堂課上完以後也的確和自己打聽到的相同可以朝者自己夢想跨一大步，但是其實更多的是被打擊的心情，想者自己是不是其實也沒那麼喜歡資工，想者自己是不是其實不適合讀資工，甚至是想者如果不適合乾脆明天我就去辦修學不要繼續讀下去了，又或是反正就是去台積電輪班阿，對XDD，一開始自己想了很多，但是其實自己根本還沒去「嘗試」。（這邊感謝某一個後來突然被我叫出來以為我得癌症，聊一聊的朋友XDD)
克服了自己心裡的恐懼感，再來老師每一堂的實做，都會讓我崩潰一次，並且深深的感受的自己到底有多少東西不會，而因為真的太多東西不會，一開始甚至連 .h 可以自己撰寫並且匯入都不知道，更別說 git 的操作又有多少到現在其實都還是不太會的東西，但是神奇的事情老師每一次的填表單我又幾乎都可以答對個幾題XD（存粹的看題目和選項和老師上課的內容猜測出題者的心裡這樣去投機的寫每一次的表單），因此絕大部分的寫出來的作業幾乎都是參(ㄑ一ㄝˋ)考(ㄑㄩˇ)別的共筆裡面根本就沒有自己的想法在裡頭，頂多就只是弄懂和寫下來學會的東西，就已經花費了我幾乎所有的時間，而這時就會體會到的是何謂「好課值得你一修在修」因為修一次根本無法跟上阿XDDDD（我前幾天才知道 for 如果兩層迴圈，即使你只有一行 coment 也一定要有括號包起來不然它就不會執行那一行程式，對，我就爛）
進階題的部份我目前根本就沒有去寫過，目前自己對到期末的期許就是至少能夠把老師教導的所有「武器」都碰過一次，讓自己即使這一學期真的跟不上就算了，並且在未來增強再來一步一步的加強自己並且可以一個一個去精通「武器」，把自己成為一個可以擔負的起自己夢想的人，而不是像是現在只能用嘴談夢想的人。
最後一段讓我小抱怨一下，老師開學的時候雖然說只要有 「GUTS」 就可以來修這堂課，但是我建議還要有一點基礎再來修，這樣你才可以讀到那些精華，那些這堂課最有價值的地方，而不是只能像我現在一樣都在熟悉環境和看基礎題偶爾碰一點進階題，但是進階題的實驗以及思考比較部份會讓你真正體驗到「資工」真的很有趣，而且可以很厲害，真的因為你而改變那麼一點點點點點點這個世界。（快還要更快）
:::
## Q1. 浮點數除法程式: (`fdiv.c`)
### 程式目的
$orig$(浮點數) $\div slot$ (整數)
```cpp=
#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` 型態能表達的數值上界
### 程式碼解析:
#### 大致想法
當進行除法時，以 $orig \div slot$
> 依據題意， `orig` 為 被除數， `slots` 為除數
`orig` 的型態 為 雙精度浮點數 ( `double` )
1. 在二進位中我們可以很簡單的觀察**浮點數 orig**和**整數 slot** 有無 2 的因數，浮點數 orig 可以透過將 exp-1 快速達成，另外整數 slot 可以透過右移
#### 逐行分析
```cpp=4
double divop(double orig, int slots){
```
```cpp=5
if (slots == 1 || orig == 0)
return orig;
```
* 當 除數為 1 或 被除數為 0 時，除法結果還是被除數
```cpp=7
int od = slots & 1;
double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);
```
* 若 被除數 與 除數 都有 2 因子 ，可將兩者都先除以 2 ，除法結果不變
:::danger
其實這邊說除法結果不變，不是很精確，正確來說是在有限制的情況下結果會不變。
e.g. div(8.0 , 4) = div(4.0 , 2) 這種情況其實很好理解結果會相同
但是 e.g. div(3.1415926 , 3) $\neq$ 預期結果 (1.0471975333..3) ， 實際結果 (0.628319)
我測試的臨界值就在此，如果是 3.141592 結果會正確。
在我想要繼續深究的試者把 div(3.1415926 , 3) 過程中的 `result` 和 `slots` 印出來檢查發現了一件神奇的是情他最多就只會印出小數點後6位，因此我猜測只要輸入的 `orig` 超過此範圍就都會有誤差。
詳細數據還需更多實驗
* `od` 判斷 除數 `slots` 是否為 偶數
* 若 `slots` 為偶數，則 `od` 為 0
* 若 `slots` 為奇數，則 `od` 為 1
* 若 除數 `slots` 為 偶數，則 被除數 `orig` 與 除數 `slots` 同除以 2，除法結果不變
* `divop(double orig, int slots)` = `divop(orig / 2, slots >> 1)`
* 因此， **D1 = `2`**
```cpp=8
double result = divop(orig / D1, od ? (slots + D2) >> 1 : slots >> 1);
if (od)
result += divop(result, slots);
```
* 先看一個數學推導：
$A: floating \space number, \space n: odd \space integer$
$原 - \space 修正(除數 +1)$
$A ÷ n$ ${ - }$ $A ÷(n+1)$
$=$ $\dfrac{A}{n}$ - $\dfrac{A}{n+1}$
$=$ $\dfrac{A(n+1)}{n(n+1)}$ - $\dfrac{An}{(n+1)n}$
$=$ $\dfrac{A}{n(n+1)} = 誤差值$
並且我們可以發現 $\dfrac{A}{n+1}$ 這個結果就是我們在前面就得到的值 $result$
也就是說誤差值就是原先到的 $\dfrac{result}{n}$ ，$n$ 為原先傳入函式的除數。
結論 :
$原 = 修正 + 誤差值$
```c=11
result += divop(result, slots);
```
## Q2. 浮點數開二次方根的近似值
* [IEEE 754 Standard: Scope](https://hackmd.io/tbHqxe19SdafIq0XdBnJzQ?both#IEEE-754-Standard-Scope)
* [浮點數 與 定點數](https://hackmd.io/tbHqxe19SdafIq0XdBnJzQ?both#%E6%B5%AE%E9%BB%9E%E6%95%B8%E5%92%8C%E5%AE%9A%E9%BB%9E%E6%95%B8)
* 假設 float 為 32-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式
```cpp=
#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;
}
```
:::info
* [C 語言中的typedef、struct、與union](https://zhung.com.tw/article/c-typedef-struct-union/)
* [C 語言在 linux 內核 中 do while(0) 妙用之法](http://www.aspphp.online/bianchen/cyuyan/C/cyyrm/201701/398.html)
:::
* 單精度浮點數
* 程式碼解析：
```cpp=34
EXTRACT_WORDS(ix0, x);
```
* 將 浮點數 `x` 型態轉換成 32 bit 的整數 `ix0`
* 方便接下來的 bitwise 操作
##### Part1. 根號運算的特殊狀況處理
```cpp=36
/* take care of zero */
if (ix0 <= 0) {
if ((ix0 & (~sign)) == 0) //若 ix0 為零
return x; /* sqrt(+-0) = +-0 */
if (ix0 < 0) //若 ix0 為負數
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 )`
* ==負數開根號後為虛數==，應該在實數運算中排除
```cpp=43
/* take care of +INF and NaN */
if ((ix0 & 0x7ff00000) == 0x7ff00000) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
/* 若 ix0 為 無限 或 NaN 的情況
return x;
}
```
* $\sqrt{\infty} = \infty$
* $\infty$ 以浮點數表示 為 `( 0 11111111 0000...0000 )`
* $NaN 開根號 仍為 NaN$
* $NaN$ 以浮點數表示 為 `( * 11111111 ****...**** )`
##### Part2. 處理 Denormalize 型 並 取出 Exponent 與 Fraction
```cpp=49
m = (ix0 >> 23);
if (m == 0) { /* subnormal x */
for (i = 0; (ix0 & 0x00800000) == 0; i++)
ix0 <<= 1;
m -= i - 1;
}
m -= 127; /* unbias exponent */
```
F : 23 bits
E : 8 bits
**特別注意 : `ix0` 的形式是 IEEE754 的樣子**
* 因此`#49` : `m` 為 浮點數表示 的 Exponent 部份
* 若 `m` 為 0，代表該浮點數為 Denormalized
( i.e. `( * 00000000 ****...**** )` $=0.$**** ... **** $\times 2^{-126}$ )
* 此時，我們還能進行 subnormalize ( `#51` ~ `#53` )
>舉例：
> ( 0.00001***...*** )~2~ = ( 1.*** ... *** )~2~ * $2^{-5}$ = 0 0000000==0== 00001***
> for 迴圈每一次往左位移 1 bit 直到 標示處為 ==1==
> 因此`i` = 6 ( for 迴圈 i++ 特性 導致最後還會多加到一次，因此最後要扣掉那一次 )
> `m` = m - i - 1 = 0 - 6 - 1 = -5
* 0x00800000 = ( 0000 0000 1000 0000 ... 0000 )~2~
= `( 0 00000001 0000...0000 )`
* `#55` : 由於浮點數使用 [Exponent Bias](https://en.wikipedia.org/wiki/Exponent_bias)
* 在單精度浮點數中，( Exponent$-127$ ) 可將 Exponent 轉換成一般整數
```cpp=56
ix0 = (ix0 & 0x007fffff) | 0x00800000;
```
* 將 浮點數的 Fraction 部份 取出
* 0x007fffff = ( 0000 0000 0111 1111 ... 1111 )~2~
= `( 0 00000000 1111...1111 )`
* 走到這裡，`ix0` 只剩 Fraction 的部份
* 從數學上可理解成：`ix0` $= 1.Fraction$
* $1 \leq$ `ix0` $<2$
處理 Denormalize floating number 其格式，透過在 m 存下原本 E 值，ix0 存下 normalize 表示式，透過這個方法讓每個值都可以用 normalize floating number 格式處理。
==這邊不懂，==
##### Part3. 小數開根號
```cpp=57
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{{2x} \times 2^{m-1}}=\sqrt{2x} \times 2^{\left\lfloor \frac{m}{2} \right\rfloor}$
`#57~58` 就是在把奇數的情況下 $ix0*2$
* 若 $m$ 為 偶數，則 $\sqrt{x \times 2^{m}}=\sqrt{x} \times 2^{\frac{m}{2}}$
`#60` 最後統一不管是 奇數 or 偶數將 $m \div 2$
==後面的部分看不懂==
## Q3. LeetCode [829. Consecutive Numbers Sum](https://leetcode.com/problems/consecutive-numbers-sum/)
* Given a positive integer `N`, how many ways can we write it as a sum of consecutive positive integers?
#### 數學推導 I
```cpp=
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;
}
```
* 思路 $k \epsilon 任意自然數$
>Let $N = k + (k+1) + (k+2) + (k+3) + … + (k+i-1) = i k + (1+2+3+… + i- 1)$
Let $sum(i) = (1+2+3+…+i-1)$, then we have
$N = sum(i) + i k => ik = N - sum(i)$
代表 $N = ik +sum(i)$ 有解
* 程式碼解析：
```cpp=5
int ret = 1;
int x = N;
for (int i = 2; i < x; i++) {
x += ZZZ;
if (x % i == 0)
ret++;
}
```
* 假設等差數列 ( 差值為 1 ) 從 $k$ 開始，且個數共有 $i$ 個則此等差數列為
$$
k, k+1,k+2,...,k+(i-1)
$$
且令其 Input 為 $N$ ， $k \epsilon 任意自然數$
則
$$
N=ik+\frac{(i-1)i}{2} \\
\Rightarrow ik =N - \frac{(k-1)k}{2} \\
\Rightarrow ik =N - [1+2+...+(i-1)] \\
\Rightarrow \begin{Bmatrix} N - [1+2+...+(i-1)] \end{Bmatrix}\space mod \space i = 0\\
$$
* (N - 1) % 2 , 此時(i = 2)
((N - 1) -2) % 3 , 此時(i = 3)
(((N - 1) -2) -3) % 4 , 此時(i = 4)
...類推
* 故 **ZZZ = `1 - i`** (e.g. i = 2,3,4..)
* `ret` : 當 `i` 符合的時，則 `ret++`
* 因為 $k \in 任意自然數$，因此程式當中只需要判斷是否存在，因此不需要
* 效能評估
## Q3. 進階題：嘗試更有效率的實作
#### 修改迴圈中止條件
```cpp
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 自然數$
$$
N - \frac{(k-1)k}{2} >0 \\
\Rightarrow k< \sqrt{2N}\\
\Rightarrow k^2 < 2N\\
$$
Sign in with Wallet
Connect another wallet