2018q1 Homework2 (assessment)
===
contributed by <`blake11235`>
---
# 重做第一週測驗題
>好東西分享，很完整的[HackMD功能介紹](https://hackmd.io/TKNuhom7S62OV6bDyBglXA?both)
## 測驗1
考慮以下程式碼:
```c=
#include <stdlib.h>
int isMultN(unsigned int n) {
int odd_c = 0, even_c = 0; /* variables to count odd and even SET bits */
if (n == 0) // return true if difference is 0.
return 1;
if (n == 1) // return false if the difference is not 0.
return 0;
while (n) {
if (n & 1) // odd bit is SET, increment odd_C
odd_c++;
n >>= 1;
if (n & 1) // even bit is SET, increment even_c
even_c++;
n = n >> 1;
}
/* Recursive call till you get 0/1 */
return(isMultN(abs(odd_c - even_c)));
}
```
其作用為檢查輸入整數是否為 N 的倍數
### 題目研讀
首先先觀察程式碼
```
while (n) {
if (n & 1) // odd bit is SET, increment odd_C
odd_c++;
n >>= 1;
if (n & 1) // even bit is SET, increment even_c
even_c++;
n = n >> 1;
}
```
從這段可以得知他是在計算 n 用 binary 表示的話，奇數位與偶數位的 1 的數目。先是判斷最右邊的那位數，看是否為 1 ，接著向右位移再判斷是否為 1 ，一直操作到 n 為 0 ，代表全部計算完畢。
```
return(isMultN(abs(odd_c - even_c)));
```
最後將 odd bit count 和 even bit count 的差值再次輸入函式做遞迴，直到數值為 0 或 1 。也就是最終的判斷依據，是看**奇數位 1 的數與偶數位 1 的數兩者要一樣**。
既然這段程式碼是依據各 bit 在判斷，那麼就把數值列出來啦~
|數值|binary|奇偶位 1 的數的差值|
|:----:|:----:|:----:|
|1|0001|1|
|2|0010|1|
|3|0011|0|
|5|0101|2|
|7|0111|1|
|11|1011|1|
|13|1101|1|
搭啦~ 答案很明顯的就是 3 啦~
因為只有 3 他的奇位 1 的數目和偶位 1 的數目一樣，所以只能選他了。
若差值是 1 以上的，會再次進入函式遞迴，直到最後的數值是 1 或者 0 。
### 推導
但事情有那麼順利嗎？會不會有除了 3 以外的質數，他的奇位 1 的數目和偶位 1 的數目也一樣，是不是就爆炸了？那麼我們就來證證看囉！
首先，這個題目讓我想到了一個經典個國中數學題目：**如何判斷一個數 3 的倍數？**
只需要把每個位數的值都加起來，如果為 3 的倍數的話，那麼此數就是 3 的倍數。
推導方法很簡單，這裡就有點了的講...
那麼就從那個觀念延伸，來這提證明。
假設一個二進位的數值 ABCDEFGH ，代表的是 $A\cdot2^7+B\cdot2^6+C\cdot2^5+D\cdot2^4+E\cdot2^3+F\cdot2^2+G\cdot2^1+H\cdot2^0
\\=A\cdot128+B\cdot64+C\cdot32+D\cdot16+E\cdot8+F\cdot4+G\cdot2+H\cdot1
\\=(A\cdot129-A)+(B\cdot63+B)+(C\cdot33-C)+(D\cdot15+D)+(E\cdot9-E)+(F\cdot3+F)+(G\cdot3-G)+(H\cdot0+H)
\\=(A\cdot129+B\cdot63+C\cdot33+D\cdot15+E\cdot9+F\cdot3+G\cdot3+H\cdot0)+(-A+B-C+D-E+F-G+H)$
$(A\cdot129+B\cdot63+C\cdot33+D\cdot15+E\cdot9+F\cdot3+G\cdot3+H\cdot0)$ 是 3 的倍數
所以只要判斷$(-A+B-C+D-E+F-G+H)$ 也要是 3 的倍數就好了
可能是 3、6、9 ...等等，再把數值代進迴圈，直到最後奇偶位數相差為 0 ，便知道此數是 3 的倍數
### 其他質數
針對其他質數判斷是否為其倍數，2 和 5 較為簡單，就不多說了，那麼就想針對 7、11、13 做研究
再次回想起國中數學~~
* 7的倍數：由個數起每三位數字一節，各奇數節的和與偶數節的和相減，其差是7的倍數。
* 11的倍數：奇數位數字和與偶數位數字和相差為11的倍數。
* 13的倍數：由個數起每三位數字一節，各奇數節的和與偶數節的和相減，其差是13的倍數。
**目標，只使用加減法位移與迴圈**
#### 7的倍數
針對一個數於乘於 7 直式乘法，可以觀察到一個有趣現象
$ABCDEFG\cdot7 = ABCDEFG\cdot0111$
```
A B C D E F G
A B C D E F G
+ A B C D E F G
-----------------------
```
$(A), (A+B), (A+B+C), (B+C+D), (C+D+E), (D+E+F), (E+F+G), (F+G), (G)$
也就是說，把一個數除以 7 ，就是把最後一位數拿去減倒數第三位、倒數第二位和倒數第一位，減完之後數值向右位移，然後回傳給函式繼續進行，直到除了後三位的其他位數都是 0 ，最後再判斷數值是否為 7 。
```=c
int isMultN(unsigned int n)
{
if(n>>3 == 0) {
return (n==7);//判斷是否餘七
}
else{
int temp;
temp = (n&1);//取最後一位數值
n = n - temp - (temp<<1) - (temp<<2);//除於七
return(isMultN(n>>1));
}
}
```
既然都打出程式了，那麼就來比賽一下，看看是否有比 mod 運算還快，不然而外想到的這種作法沒有比較快，內心會很挫折。
所以阿我用了 clock() 來算一下時間，跟 mod 來比賽看看
```=c
#include <stdlib.h>
#include <time.h>
static double diff_in_second(struct timespec t1, struct timespec t2)
{
struct timespec diff;
if (t2.tv_nsec-t1.tv_nsec < 0) {
diff.tv_sec = t2.tv_sec - t1.tv_sec - 1;
diff.tv_nsec = t2.tv_nsec - t1.tv_nsec + 1000000000;
} else {
diff.tv_sec = t2.tv_sec - t1.tv_sec;
diff.tv_nsec = t2.tv_nsec - t1.tv_nsec;
}
return (diff.tv_sec + diff.tv_nsec / 1000000000.0);
}
int isMultN(unsigned int n)
{
if(n>>3 == 0) {
return (n==7);
}
else{
int temp;
temp = (n&1);
n = n - temp - (temp<<1) - (temp<<2);
return(isMultN(n>>1));
}
}
int main()
{
unsigned int x;
scanf("%u", &x);
struct timespec start1, end1, start2, end2;
clock_gettime(CLOCK_REALTIME, &start1);
printf("%s\n", isMultN(x)?"yes":"no");
clock_gettime(CLOCK_REALTIME, &end1);
printf("my time: %.10f sec\n\n", diff_in_second(start1, end1));
clock_gettime(CLOCK_REALTIME, &start2);
printf("%s\n", (x%7)?"no":"yes");
clock_gettime(CLOCK_REALTIME, &end2);
printf("mod time: %.10f sec\n\n", diff_in_second(start2, end2));
return 0;
}
```
結論是......
我輸了
```
1354684745
no
my time: 0.0000273730 sec
no
mod time: 0.0000087500 sec
```
好吧...先做 11 的倍數，不然做不完啦...
#### 11 的倍數
11 的倍數也可以用類似於 7 的倍數的方法，只是改變了一點數值。
```=c
int isMultN(unsigned int n)
{
if(n>>4 == 0) {
return (n==11);
}
else{
int temp;
temp = (n&1);
n = n - temp - (temp<<1) - (temp<<3);
return(isMultN(n>>1));
}
}
```
但是，這個方法他的運算時間，比 mod 慢更多，差到了100倍左右...
```
6541213
no
my time: 0.0000271370 sec
no
mod time: 0.0000034600 sec
```
~~想一個比 mod 慢的程式有屁用~~<by`e94046165`>
#### 13 的倍數
也是可以用一樣的手法但是速度一樣慢...
```
65
yes
my time: 0.0000257650 sec
yes
mod time: 0.0000036420 sec
```
所以，強者我室友`e94046165`看了我的 code 叫我把函式遞迴改成 while() 迴圈，順便嘴我的資料結構...
結果，速度快的幾乎跟 mod 一樣了！！
---
# 重做第二週測驗題
## 測驗1
以類似 CS:APP 3/e 第 80 頁針對 8 位元浮點格式的示例，擴充至單精度浮點數
|位表示(b x x b)|指數|小數|值|描述|
|:---:|:---:|:---:|:---:|:---:|
|bin hex hex bin|e, E, 2^E^|f, M|V, 十進位|
|0 00 00000 000|0, -126, 2^-126^|0, 0|0, 0.0|零|
|0 00 00000 001|0, -126, 2^-126^|2^-23^, 2^-23^|2^-149^, 0.00..|最小 denormalized|
|...|...|...|...|denormalized|
|0 00 11111 111|0, -126, 2^-126^|2^-1^+2^-2^...+2^-23^, 2^-1^+2^-2^...+2^-23^|...|最大 denormalized
|0 01 00000 000|1, -126, 2^-126^|0, 1|2^-126^, 0.00...|最小 normalized
|0 01 00000 001|1, -126, 2^-126^|2^-23^, 1+2^-23^|2^-126^+2^149^, 0.00...|normalized|
|...|...|...|...|normalized|
|0 fe fffff 111|254, 127, 2^127^|2^-1^+2^-2^...+2^-23^, 1+2^-1^+2^-2^...+2^-23^|...|最大 normalized|
|0 ff 00000 000|...|...|...|正無限大|
在來看題目：
```=c
#include <math.h>
static inline float u2f(unsigned int x) { return *(float *) &x; }
float exp2_fp(int x) {
unsigned int exp /* exponent */, frac /* fraction */;
/* too small */
if (x < 2 - pow(2, Y0) - 23) {
exp = 0;
frac = 0;
/* denormalize */
} else if (x < Y1 - pow(2, Y2)) {
exp = 0;
frac = 1 << (unsigned) (x - (2 - pow(2, Y3) - Y4));
/* normalized */
} else if (x < pow(2, Y5)) {
exp = pow(2, Y6) - 1 + x;
frac = 0;
/* too large */
} else {
exp = Y7;
frac = 0;
}
/* pack exp and frac into 32 bits */
return u2f((unsigned) exp << 23 | frac);
}
```
* too small
* 在程式碼的這部份當中，代表的是此數值要比最小的 denormalized 還小，也就是 < 2^-149^
2^-149^ 是由指數部份 2^-126^ 和小數部份 2^-23^ 相乘得來，所以`(x < 2 - pow(2, 7) - 23)`
+ denormalized
- 此部份代表的是非規格化的數，要符合的話必須小於最大的 normalized 2^-126^ ，也就是 exp - *Bias* = 1 - (2^8-1^ - 1)，所以`(x < 2 - pow(2, 7))`
- frac 的部份，若要表示 2^x^ 必定只會有一個 1 ，範圍在 100...000 = 2^-1^ * 2^-126^ = 2^-127^ 最大，到 000...001 = 2^-23^ * 2^-126^ = 2^-149^ 最小
右移 0 代表的數是 2^-149^ ，左移 22 代表的數是 2^-127^ ，所以`1 << (x - (2 - pow(2, 7) - 23))`
* normalized
+ 要符合規格化的數，必須比 normalized 能表示的最大數再大，也就是比 2^127^ * ( 1 + 2^-1^ + 2^-2^ + ... +2^-23^ ) 大，那麼就是 2^128^ 還小，所以`(x < pow(2, 7))`
+ exp 的部份，E = e - *Bias* = e - (2^7^ - 1)，所以`exp = pow(2, 7) - 1 + x`
* too large
+ 根據定義，無窮大的表示方式為 s 11111111 000...000，所以是`0xff`
### 延伸題
給定一個單精度的浮點數值，輸出較大且最接近的 2^x^ 值，需要充分考慮到邊界
原來是針對原本的題目做反向操作阿，研究了一下，就用 bits 來操作吧
首先針對所得到的浮點數做分類：
* 0
很簡單，輸入是 0 的時候當然就是 0 啦
* 無限大
也很簡單，當數值為 0xffffffff 就是無限大
* denormalized
比較麻煩一點，得先判斷 23~3 位是 0 ，再來判斷後面 frac 的部份，由於題目有要求到輸出較大，所以當考慮 frac 的最左邊的 1 時若其右邊還有 1 ，他的數值就必須在向左移，也就是多 2 倍
* normalized
反而更簡單了，先處理 exp 的部份，將數值向左移 23 位出現 exp 的部份，再減去 *Bias* = 127
frac 的部份，幾乎不用考慮，因為是要輸出較大，所以只要不是全為 0 ，其他都要多乘 2^1^
```=c
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
unsigned int flo2bin(float f);
int main(){
float x;
int ans = 0;
printf("input a float number：");
scanf("%f", &x);
if(!x){ //0
printf("0\n");
}
else if(x == 0xffffffff){ //infinity
printf("INFINITY\n");
}
else{
if(x >> 23 == 0){ //denormalized
x = x << 9;
while(!(x&80000000)){
x = x << 1;
ans++;
}
if(x << 1)
ans = -ans;
else
ans = -ans-1;
}
else{ //normalized
if(x << 9 == 0) //M=(1+0)=2^0
ans = (x >> 23) - 127;
else //M=(1+f)=2^1
ans = (x >> 23) - 127 + 1;
}
printf("2 e%d", ans);
}
return 0;
}
```
完成了...嗎？
```
float.c: In function ‘main’:
float.c:23:8: error: invalid operands to binary >> (have ‘float’ and ‘int’)
if(x >> 23 == 0){
```
原來， float 不能做 bitwise 的運算阿！~我真的現在才知道~
于是乎，~~我~~想到了一招，用 [**pointer**](http://edisonx.pixnet.net/blog/post/83095843-%5B%E6%B5%AE%E9%BB%9E%E6%95%B8%5D-c-%E8%AA%9E%E8%A8%80%E9%A1%AF%E7%A4%BA%E6%B5%AE%E9%BB%9E%E6%95%B816-2%E9%80%B2%E4%BD%8D)
1. 先將 float 取位址
2. 將位址內容轉型成 unsigned int* 指標
3. 將unsigned int* 指標取值
意思就是透過轉換位址的型態，將位址傳給 type size 一樣的 unsigned int ，在對它取值，這樣就可以對 unsigned int 做 bitwise 運算了
```=c
unsigned int flo2bin(float f){
unsigned int i = *(unsigned int*)&f;
return i;
}
```
最後我再多計算輸出 2^x^ 與原本浮點數的差值
>再啦幹
>[name=e94046165]
>>...謝謝...
>>[name=blake11235]
```
input a float number：0.25
2 e-2
diff: 0.000000
input a float number：0.0035
2 e-8
diff: 0.000406
input a float number：5646513.154
2 e23
diff: 2742095.000000
```
---
# 重做第三週測驗題
## 測驗一
延續 [從 $\sqrt{2}$ 的運算談浮點數](https://hackmd.io/s/ryOLwrVtz)，假設 double 為 32-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式，請補上對應數值，使其得以運作：
```=c
#include <stdint.h>
/* A union allowing us to convert between a double and two 32-bit integers.
* Little-endian representation
*/
typedef union {
double value;
struct {
uint32_t lsw;
uint32_t msw;
} parts;
} ieee_double_shape_type;
/* Set a double from two 32 bit ints. */
#define INSERT_WORDS(d, ix0, ix1) \
do { \
ieee_double_shape_type iw_u = { \
.parts.msw = ix0, \
.parts.lsw = ix1, \
}; \
(d) = iw_u.value; \
} while (0)
/* Get two 32 bit ints from a double. */
#define EXTRACT_WORDS(ix0, ix1, d) \
do { \
ieee_double_shape_type ew_u; \
ew_u.value = (d); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
} while (0)
static const double one = 1.0, tiny = 1.0e-300;
double ieee754_sqrt(double x)
{
double z;
int32_t sign = 0x80000000;
uint32_t r, t1, s1, ix1, q1;
int32_t ix0, s0, q, m, t, i;
EXTRACT_WORDS(ix0, ix1, x);
/* take care of INF and NaN */
if ((ix0 & KK1) == KK2) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN */
return x * x + x;
}
/* take care of zero */
if (ix0 <= 0) {
if (((ix0 & (~sign)) | ix1) == 0)
return x; /* sqrt(+-0) = +-0 */
if (ix0 < 0)
return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
}
/* normalize x */
m = (ix0 >> 20);
if (m == 0) { /* subnormal x */
while (ix0 == 0) {
m -= 21;
ix0 |= (ix1 >> 11);
ix1 <<= 21;
}
for (i = 0; (ix0 & 0x00100000) == 0; i++)
ix0 <<= 1;
m -= i - 1;
ix0 |= (ix1 >> (32 - i));
ix1 <<= i;
}
m -= KK3; /* unbias exponent */
ix0 = (ix0 & 0x000fffff) | 0x00100000;
if (m & 1) { /* odd m, double x to make it even */
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
}
m >>= 1; /* m = [m/2] */
/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x00200000; /* 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 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;
}
r = sign;
while (r != 0) {
t1 = s1 + r;
t = s0;
if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
s1 = t1 + r;
if (((t1 & sign) == sign) && (s1 & sign) == 0)
s0 += 1;
ix0 -= t;
if (ix1 < t1)
ix0 -= 1;
ix1 -= t1;
q1 += r;
}
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;
}
/* use floating add to find out rounding direction */
if ((ix0 | ix1) != 0) {
z = one - tiny; /* trigger inexact flag */
if (z >= one) {
z = one + tiny;
if (q1 == (uint32_t) 0xffffffff) {
q1 = 0;
q += 1;
} else if (z > one) {
if (q1 == (uint32_t) KK4)
q += 1;
q1 += 2;
} else
q1 += (q1 & 1);
}
}
ix0 = (q >> 1) + 0x3fe00000;
ix1 = q1 >> 1;
if ((q & 1) == 1)
ix1 |= sign;
ix0 += (m << KK5);
INSERT_WORDS(z, ix0, ix1);
return z;
}
```
乾這題好難，本來想偷偷看有誰有做這題的，結果沒有...
首先先解決那五格填空該填入什麼
```
/* take care of INF and NaN */
if ((ix0 & KK1) == KK2) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN */
return x * x + x;
}
```
可以觀察出此段程式碼要判斷是否為 INF 或 NAN ，而判斷標準就是數字 exp 的部份皆為 1 ，因此 KK2 為 0x7ff00000 ， mask 的部份就是取 ix0 的 62~52 位數，所以 KK1 也是 0x7ff00000
```
m -= KK3; /* unbias exponent */
```
程式碼後面的註解都說了要 unbias 也就是將偏置移動回來，對 double 雙精度而言 *Bias* = 2^k-1^ - 1 = 1023
```
/* use floating add to find out rounding direction */
if ((ix0 | ix1) != 0) {
z = one - tiny; /* trigger inexact flag */
if (z >= one) {
z = one + tiny;
if (q1 == (uint32_t) 0xffffffff) {
q1 = 0;
q += 1;
} else if (z > one) {
if (q1 == (uint32_t) KK4)
q += 1;
q1 += 2;
} else
q1 += (q1 & 1);
}
}
```
:::info
其實這題我不太懂...
:::
```
ix0 += (m << KK5);
```
最後這段，是要把 m exponent 放入 ix0 裡面，而針對 double 他的 exp 是在 62~52 位元，所以需要向左位移 51-32+1=20 位
### 分析程式碼
```=c
typedef union {
double value;
struct {
uint32_t lsw;
uint32_t msw;
} parts;
} ieee_double_shape_type;
/* Set a double from two 32 bit ints. */
#define INSERT_WORDS(d, ix0, ix1) \
do { \
ieee_double_shape_type iw_u = { \
.parts.msw = ix0, \
.parts.lsw = ix1, \
}; \
(d) = iw_u.value; \
} while (0)
/* Get two 32 bit ints from a double. */
#define EXTRACT_WORDS(ix0, ix1, d) \
do { \
ieee_double_shape_type ew_u; \
ew_u.value = (d); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
} while (0)
```
* 這部份的程式碼，是用 struct 的方式解決了我之前所遇到的問題， float 無法做 bitwise 的運算。宣告一個 union ，裏面包含了一個 double 的值與真正可以操作的 uint32_t 部份。之後再 define EXTRACT_WORDS 和 INSERT_WORDS，讓一個 double 的數值，第 63 位道 32 位能夠輸入到 uint_32 ix0 ， 31 位到 0 位輸入到 uint_32 ix1。
```=c
/* take care of INF and NaN */
if ((ix0 & 0x7ff00000) == 0x7ff00000) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN */
return x * x + x;
}
/* take care of zero */
if (ix0 <= 0) {
if (((ix0 & (~sign)) | ix1) == 0)
return x; /* sqrt(+-0) = +-0 */
if (ix0 < 0)
return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
}
```
* 用條件判斷來分別處理 INF 、 NAN 、 ZERO 這幾種狀況
* (ix0 & 0x7ff00000) 是代表 double 的 exp 部份，若皆為 1 則代表 INF or NAN ，為了要使`sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN` 因此使用 x * x + x 來表示，[參考](http://www.cnblogs.com/dosrun/p/3908617.html)
* 要處理輸入為 0 的部份，則判斷 ix0 的後 31 位及 ix1 的所有位數都要是 0 ， ix0 的第一位是 sign 用來判斷 +0 還是 -0
* 若輸入是小於 0 ，則可以用 0/0 來使 return 為 NAN
:::info
這裡不太清楚為何能讓 透過 (x - x) / (x - x); 使得 sqrt(-ve) = -NaN ，是否因為產生了 -0/-0 ？
:::
```=c
/* normalize x */
m = (ix0 >> 20);
if (m == 0) { // subnormal x
while (ix0 == 0) {
m -= 21;
ix0 |= (ix1 >> 11);
ix1 <<= 21;
}
for (i = 0; (ix0 & 0x00100000) == 0; i++)
ix0 <<= 1;
m -= i - 1;
ix0 |= (ix1 >> (32 - i));
ix1 <<= i;
}
```
* 將非規格化的數規格化。這裡的 m 代表的是 sign exp 的部份，若為 0 則需要先規格化。
```=c
m -= 1023; /* unbias exponent E = e - Bias */
ix0 = (ix0 & 0x000fffff) | 0x00100000;
if (m & 1) {
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
}
m >>= 1;
```
* 針對一個 normalized 的數值可以知道 $V = (-1)^s\cdot M\cdot 2^E = (1+f)\cdot(e-Bias)$
* 將 e 減去 1023 之後就可得到真正二的幕次
* 把 ix0 取出 frac 的部份，也就是第 51 到 32 位，之後再加上 2^0^ ，得到 (1+f) = M
* $V = M\cdot 2^E$ ，若 E 為偶數，則 $\sqrt{M\cdot 2^E} = 2^{E/2}\cdot\sqrt{M}$ 所以要將 E 化為偶數，所以當 E 為奇數 (m & 1)，將 ix0 和 ix1 的部份乘於 2 ， m 可以直接右移 1 代表除於 2 ，指數部份開根號就完成了。
```=c
/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x00200000; /* 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 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;//r one time / 2
}
r = sign;
while (r != 0) {
t1 = s1 + r;
t = s0;
if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
s1 = t1 + r;
if (((t1 & sign) == sign) && (s1 & sign) == 0)
s0 += 1;
ix0 -= t;
if (ix1 < t1)
ix0 -= 1;
ix1 -= t1;
q1 += r;
}
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;
}
```
* 原本以為這是用二分法去實現，因為看到了 ix0 和 t 比較大小，所以一直往二分法想，但結果不對，又不可能是十分法，牛頓法也不像，最後找了好久...原來是[這個](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_(base_2))阿。
* $(a\cdot2^a+b\cdot2^b+c\cdot2^c+d\cdot2^d...)^2=\\2^{2a}\cdot a^2 \\+(2^{2a-1}\cdot2a+2^{2a-2}\cdot b)\cdot b\\+(2^{2a-2}\cdot2a+2^{2a-3}\cdot 2b +2^{2a-4}\cdot c)\cdot c\\+(2^{2a-3}\cdot2a+2^{2a-4}\cdot 2b +2^{2a-5}\cdot 2c+2^{2a-6}\cdot d)\cdot d\\+ ...$
* 分成 ix0 和 ix1 的部份執行，執行結束後可以得到 q 和 q1 兩個 root
```
/* use floating add to find out rounding direction */
if ((ix0 | ix1) != 0) {
z = one - tiny; /* trigger inexact flag */
if (z >= one) {
z = one + tiny;
if (q1 == (uint32_t) 0xffffffff) {
q1 = 0;
q += 1;
} else if (z > one) {
if (q1 == (uint32_t) 0xfffffffe)
q += 1;
q1 += 2;
} else
q1 += (q1 & 1);
}
}
```
* 這部份是用來捨入的，只是詳細的原因我不太懂...
```
ix0 = (q >> 1) + 0x3fe00000;
ix1 = q1 >> 1;
if ((q & 1) == 1)
ix1 |= sign;
ix0 += (m << KK5);
INSERT_WORDS(z, ix0, ix1);
return z;
```
* 終於到最後了，此時把 q 的部份放進 ix0 並且把之前扣掉的 Bias 加回去，在把 q1 也放進去 ix1，特別注意由於在前面有把 frac 的部份向右移，這裡要全部移回去。最後把 exp m 的部份加到 ix0 裏面就完成了。
### float版本
```=c
#include <stdint.h>
/* A union allowing us to convert between a double and two 32-bit integers.
* Little-endian representation
*/
typedef union {
float value;
struct {
uint16_t lsw;
uint16_t msw;
} parts;
} ieee_float_shape_type;
/* Set a float from two 16 bit ints. */
#define INSERT_WORDS(f, ix0, ix1) \
do { \
ieee_float_shape_type iw_u = { \
.parts.msw = ix0, \
.parts.lsw = ix1, \
}; \
(f) = iw_u.value; \
} while (0)
/* Get two 16 bit ints from a float. */
#define EXTRACT_WORDS(ix0, ix1, f) \
do { \
ieee_float_shape_type ew_u; \
ew_u.value = (f); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
} while (0)
static const float one = 1.0, tiny = 1.0e-300;//TODO
double ieee754_sqrt(float x)
{
float z;
int16_t sign = 0x8000;
uint16_t r, t1, s1, ix1, q1;
int16_t ix0, s0, q, m, t, i;
EXTRACT_WORDS(ix0, ix1, x);
/* take care of INF and NaN */
if ((ix0 & 0x7f80) == 0x7f80) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN */
return x * x + x;
}
/* take care of zero */
if (ix0 <= 0) {
if (((ix0 & (~sign)) | ix1) == 0)
return x; /* sqrt(+-0) = +-0 */
if (ix0 < 0)
return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
}
// normalize x
m = (ix0 >> 7);
/* if (m == 0) { // subnormal x
while (ix0 == 0) {
m -= 21;
ix0 |= (ix1 >> 11);
ix1 <<= 21;
}
for (i = 0; (ix0 & 0x00100000) == 0; i++)
ix0 <<= 1;
m -= i - 1;
ix0 |= (ix1 >> (32 - i));
ix1 <<= i;
}
*/
m -= 127; /* unbias exponent E = e - Bias */
ix0 = (ix0 & 0x007f) | 0x0080; /* M = 1 + f = (2^0 == |00100000) + f */
if (m & 1) { /*if m is odd , double x to make it even , in order to make 2^E--->2^(E/2) , E need to be even 2^(E+1)*M = 2^(E)*M*2 */
ix0 += ix0 + ((ix1 & sign) >> 15);/* if need to carry in, ix1 & sign will be 1, ix1*2 overflow is OK */
ix1 += ix1;
}
m >>= 1; /* m = [m/2] 2^(E+1)--->2^(E) */
/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1 & sign) >> 15);
ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x0100; /* 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 + ((ix1 & sign) >> 15);
ix1 += ix1;
r >>= 1;//r one time / 2
}
r = sign;
while (r != 0) {
t1 = s1 + r;
t = s0;
if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
s1 = t1 + r;
if (((t1 & sign) == sign) && (s1 & sign) == 0)
s0 += 1;
ix0 -= t;
if (ix1 < t1)
ix0 -= 1;
ix1 -= t1;
q1 += r;
}
ix0 += ix0 + ((ix1 & sign) >> 15);
ix1 += ix1;
r >>= 1;
}
/* use floating add to find out rounding direction beacuse we have a bit that need to be rounding*/
if ((ix0 | ix1) != 0) {
z = one - tiny; /* trigger inexact flag */
if (z >= one) {
z = one + tiny;
if (q1 == (uint16_t) 0xffff) {
q1 = 0;
q += 1;
} else if (z > one) {
if (q1 == (uint16_t) 0xfffe)
q += 1;
q1 += 2;//make q1 == 0
} else
q1 += (q1 & 1);
}
}
ix0 = (q >> 1) + 0x3f00;//Bias exponent E = e - 1023
ix1 = q1 >> 1;
if ((q & 1) == 1)
ix1 |= sign;
ix0 += (m << 7);
INSERT_WORDS(z, ix0, ix1);
return z;
}
int main(){
float f;
scanf("%f", &f);
printf("%f\n", ieee754_sqrt(f));
}
```
感覺不難，只需要把數值改成對應 float 格式就好。
* 綜觀這一題，我想這就是老師所說的**好的程式碼**，能夠適應任何情況，不需要其他函式庫，透過簡單的 bitwise 操作，沒有複雜的乘除法，還考慮正負零、無限、非數等情況，捨入的部份也考慮到不同的捨入可能方式，這就是一段有價值的程式碼。
Sign in with Wallet
Connect another wallet