--- tags: sysprog2018 --- # 2018q1 Homework2 (assessment) contributed by <`HexRabbit`> ## 第三周測驗題 延續 [從 $\sqrt{2}$ 的運算談浮點數](https://hackmd.io/s/ryOLwrVtz)，假設 double 為 64-bit 寬度，考慮以下符合 IEEE 754 規範的平方根程式，請補上對應數值，使其得以運作： ```clike #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; } ``` ## 想法 & 思考 > 為避免畫面紊亂，部分解說直接記錄於註解中 > [color=orange] **64-bit double presition float data fromat** | Sign | Exponent | Fraction | | --- | --- | --- | | 1 | 11 | 52 | $Exponent > 0: Sign * 2^{Exponent-1023} * 1.Fraction$ $Exponent == 0: Sign * 2^{-1022} * 0.Fraction$ 題目落落長，總之先瞭解各個function與macro在做些什麼吧 - `ieee_double_shape_type`: `union` 用來簡化操作double - `INSERT_WORDS(d, ix0, ix1)`: `macro` 將兩個32bit的int ix0,ix1放入一個double d 中，ix0在高位 ix1在低位 (只適用於little-endian的機器上) - `EXTRACT_WORDS(ix0, ix1, d)`: `macro` 功能與INSERT_WORDS相反 - `ieee754_sqrt(x)`: `function` return x 的平方根 瞭解這些後就可以開始解題了，題目首先把 x 提取至 ix0, ix1中，再來出現了第一個問題 ### Handle INF & NaN ```clike /* take care of INF and NaN */ if ((ix0 & KK1) == KK2) { /* sqrt(NaN) = NaN, sqrt(+INF) = +INF, sqrt(-INF) = sNaN */ return x * x + x; } ``` 首先觀察到了一個奇怪的詞`sNaN`，查詢維基百科後發現NaN透過判斷Fraction的第一位是否為零，可以分為"signal NaN" 與 "quiet NaN"兩種，差異在於前者在FPU運算時會拋出 Exception (Linux下是 [SIGFPE](http://pubs.opengroup.org/onlinepubs/009696699/basedefs/signal.h.html)) 而後者不會。 另外也發現對於浮點數 X 的任何有關 NaN 的數學判斷式，只有 != 會成立 | Comparison | NaN >= X | NaN <= X | NaN > X | NaN < X | NaN == X | NaN != X | --- | --- | --- | --- | --- | --- | --- | Result | Always False | Always False | Always False | Always False | Always False | Always True | 不過只要瞭解這段式子是用來處理 INF 和 NaN， 就可以知道 KK1 的作用為 mask Exponant 部分以判斷 INF 和 NaN， 而 INF 和 NaN 在 Exponent 中的 bit 皆為 1，也就知道 KK2 了。 ::: success KK1 = 0x7FF00000 KK2 = 0x7FF00000 ::: ### Handle zero and negative 程式的下一部份是處理 x 為零(其實還有負數)的情形，稍微將其簡化 ```clike /* take care of zero and negative */ if (ix0 <= 0) { if (x == 0) return x; /* sqrt(+-0) = +-0 */ if (ix0 < 0) return (x - x) / (x - x); /* sqrt(-ve) = sNaN */ } ``` :::warning 不太瞭解 (x - x) / (x - x) 的意思，這跟 0 / 0 差在哪? ::: > 提示： x-x 真的 = 0 嗎？ > [name=ryanpatiency][color=green] > 像是 (-1.0 - -1.0) 這應該是 0 吧 > [name=HexRabbit][color=orange] :::info 新發現，0 / 0 會觸發 FPU 的 exception，而 0.0 / 0.0 不會 ::: 查閱 [IEEE 754 浮點數標準](http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4610935) 7.2 節 Invalid operation 部分， :::info **For operations producing results in floating-point format, the default result of an operation that signals the invalid operation exception shall be a quiet NaN** that should provide some diagnostic information (see 6.2). These operations are: - a) any general-computational or signaling-computational operation on a signaling NaN (see 6.2), except for some conversions (see 5.12) - b) multiplication: multiplication(0, ∞) or multiplication(∞, 0) - c) fusedMultiplyAdd: fusedMultiplyAdd(0, ∞, c) or fusedMultiplyAdd(∞, 0, c) unless c is a quiet NaN; if c is a quiet NaN then it is implementation defined whether the invalid operation exception is signaled - d) addition or subtraction or fusedMultiplyAdd: magnitude subtraction of infinities, such as: addition(+∞, −∞) - e) **division: division(0, 0)** or division(∞, ∞) - f) remainder: remainder(x, y), when y is zero or x is infinite and neither is NaN - g) squareRoot if the operand is less than zero - h) quantize when the result does not fit in the destination format or when one operand is finite and the other is infinite ::: 發現若一個運算的結果是以浮點數格式呈現， 則運算時遇到 invalid operation 時將不會觸發 Exception ( 即為 quiet NaN ) ### Fraction normalization 再來是對 Fraction 部分做 normalize， 將形式全部改為`1.XXXXXXXXXX`，並提取出正確的 Exponent 做 sqrt ```clike /* normalize x */ m = (ix0 >> 20); /* 擷取Exponent部分 */ if (m == 0) { /* subnormal x */ /* 移動Fraction直到指數位調整為1 -> 1.fraction (已經失去指數的意義，只為了後續計算使用)*/ /* 先調整到ix0不為零，使用21是因為若ix1的MSB是1的話就不用繼續操作了 */ while (ix0 == 0) { m -= 21; /* 把乘上的次方數從m中扣除 */ ix0 |= (ix1 >> 11); ix1 <<= 21; } /* 繼續調整到指數位為1 */ 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; /* 去除ix0的 Exponent 其餘部分，只留下1 */ if (m & 1) { /* 將奇數指數位乘入 Fraction 部分 */ ix0 += ix0 + ((ix1 & sign) >> 31); /* 可能進位 */ ix1 += ix1; } m >>= 1; /* 指數位sqrt m = [m/2] */ ``` ::: success KK3 = 1023 ::: :::info 常常出現的程式碼片段，會將 Fraction 部分左移一位 ``` ix0 += ix0 + ((ix1 & sign) >> 31); ix1 += ix1; ``` ::: ### Digit-by-digit square root method ```clike /* 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; /* 實際上要考慮 r*(s0 + r) <= ix0 對於所有 r，但是因為 r 不是 0 就是 1， * 所以先考慮 r=1 的情況，再以 t <= ix0 判斷是否繼續操作 */ if (t <= ix0) { s0 += 2*r; /* 以更易讀的做法取代 s0 = t + r; */ ix0 -= t; q += r; } /* Fraction 部分向左位移，r向右位移，可以讓產生出來的q保持連續 */ 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; } ``` 本段是對小數部分 ( 1.Fraction ) 做 sqrt， 我自己看了幾個小時其實看不太出甚麼端倪，直接代數字進去做也是霧煞煞 後來在網路上找到一個[影片](https://www.youtube.com/watch?v=G_4HE3ek4c4)，才發現原來有這樣一種奇技淫巧可以對應這裡的實作 再翻閱 wiki 上 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 的部分，有提到此作法的原理和證明， 試著套用到二進位上: :::info 假設 $Z$ 是一個有兩個位數的二進位數 $'XY'$，則 $Z^{2} = (2*X+Y)^{2} = 4*X^{2} + 2*2*XY + Y^{2}$ 所以若是可以找到一個 $X$ 使 $4*X^{2}$ 小於等於 $Z^{2}$，也就是$X^{2}<=Z^{2}$的前兩位，即為找到$X$ (因為乘以 4 與左移 2 位相等) 接著從 $Z$ 中扣除 $4*X^{2}$ 後得到 $Z2$， 接著只要想辦法找出一個 $Y$ 滿足 $Y(2*2*X+Y) = Z2$ 即可 這樣便可以繼續推廣到 N 位 ::: 若以圖像表示程式碼中變數和算式的關係的話 ``` >>>>>>> _ ------------------------ _ | 01 .10 11 10 01 ``` ``` >>>>>>> 1 _ ------------------------ 1 | 01 .10 11 10 01 (1*2 = 10) - 1 ------ 10_ | 10 ``` ``` >>>>>>> 1 0 _ ------------------------ 1 | 01 .10 11 10 01 - 1 ------ 100 | 10 (0*2 = 0) - 0 ------ 100_ |10 11 ``` ``` >>>>>>> 1 0 1 X _ <=== q, q1 ------------------------ 1 | 01 .10 11 10 01 - 1 ------ r = any single digit in q 100 | 10 - 0 ------ 1001 |10 11 (1*2 = 10) - 10 01 ------ t ===> 1010X | 10 10 (X*2 = AB) - (X*1010X) (10100+AB = ABCDE) ------ ^ ABCDE_ | ..... | s0 ``` :::warning 本段開始的這個左位移在我看來是不必要的， ``` ix0 += ix0 + ((ix1 & sign) >> 31); ix1 += ix1; r = 0x00200000; ``` 可以將 r 的值修改為 0x00100000 並和後一段中出現的 >> 1 一併移除 ``` ix0 = (q >> 1) + 0x3fe00000; ix1 = q1 >> 1; ``` 可是結果會出現微小誤差($10^{-14}$)，仍待釐清 ::: :::success => 為了多計算一位，做完 rounding 後能夠更加逼近正確值 ::: ### Rounding 再來是做捨入的部分，在上一段的最初將需要計算的部分多左移了一位， 所以會多出一位可以做捨入，讓數值更加精確。 ``` /* use floating add to find out rounding direction */ /* 不為零代表不是完全平方數 */ if ((ix0 | ix1) != 0) { z = one - tiny; /* trigger inexact flag */ if (z >= one) { z = one + tiny; /* 全為1的情況要進位(最後一位四捨五入) */ 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); } } ``` :::warning 實在不是很瞭解這是在做甚麼，還要去查查 ```clike z = one - tiny; /* trigger inexact flag */ if (z >= one) ``` ::: ### Finalize 最後將 ix0 加上 bias 並且把 Exponent 和 Fraction 放回來，再跟 ix1 組裝即完成作業。 ```clike /* 右移丟棄 Fraction 最後一位 */ ix0 = (q >> 1) + 0x3fe00000; /* +1022+1 (q 本身在指數位有 1) */ ix1 = q1 >> 1; if ((q & 1) == 1) ix1 |= sign; ix0 += (m << KK5); INSERT_WORDS(z, ix0, ix1); ``` :::success KK5 = 20 :::