contributed by < sammer1107
>
進階電腦系統理論與實作
#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);
if (od)
result += divop(result, slots);
return result;
}
od
判斷除數是否為奇數
slot
加一再除以二。如此的話,除出來的結果會比較小,之後需要加回來。divop(orig / 2, od ? (slots + 1) >> 1 : slots >> 1)
divop(result, slot)
原題目:
#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 = 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;
}
/* 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;
}
37~42: 這裡一並處理數字為 0 或數字為負的情況。
<= 0
來判斷。44~47: 這裡要處理 NaN 與 Inf,兩者在 exponent 的部份皆是 255,所以要透過 & 0x7F800000
來看 exponent 是不是全 1 。這部份正確的寫法應該是
if ((ix0 & 0x7f800000) == 0x7f800000) {
/* sqrt(NaN) = NaN, sqrt(+INF) = +INF,*/
return x;
}
48~60: 這裡把指數的部份單獨取出來做處理,因為我們已經去掉負數的情況,所以可以直接 >> 23
來取得指數。
/* 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] */
ix0
向左移,直到 mantissa 中有 1 overflow 到 exponent 中,並以負數紀錄下相對的指數在 m 中。m
現在表示真實的指數。& 0x007fffff
只留下 mantissa 的部份,然後再 | 0x00800000
來補回隱藏的 1 ,如此現在 ix0
就反映出實際的值了,像 fixed point 的感覺。>> 1
),但當 m
為奇數時,無法被 2 整除,因此我們要把 ix0
乘以二,如此指數減 1 ,現在就可以直接用 shift 做除以 2 了。62~76: 這裡用的開根號方式類似手算十進位開根號的方式,只不過改成是在 2 進位做。
ix0
這裡所代表的值,在 56 之後會佔用 24 bit,再經過 58,63 行之後最多佔用 26 bit。x.xx...
,則平方後小於 ix0
是 25 bit 或 26 bit,所以第一個處理的位置是 bit 25 (從 1 開始數),對應到 r = 0x01000000
這個答案。之後一次要跳兩格,但是為了計算方便,74~75 把 ix0 << 1
然後 r >> 1
所以和跳兩格是等效的。r
代表目前正在處理的 bit 位置,68,69 行做的事情就是測試下一個開出來的值是不是 1 ,如果是就:
2r
加到 s0
(68,70 總共加了兩個 r),這對應到 10 進位中把最右邊的數字再乘以 2 的動作。t
)與新的 bit (ix0
) 扣掉q += r
)
t = s0 + r;
if (t <= ix0) {
s0 = t + r;
ix0 -= t;
q += r;
}
q
應該會是一個 25 bit 的數字,其中 24 個 bit 皆在小數點底下,這多出來的一個 bit,讓我們之後可以做 rounding。rounding: 這裡我們透過把 1 加減很小的數來判斷系統使用的 rounding mode。
+ 如果減了之後小於 1 ,代表系統是 round down,而因為 ix0 > 0
,所以我們目前的答案已經是偏小的,所以不需再做動作。多餘的 1 等一下就會被 shift 掉。
+ 如果減了之後等於 1 ,且加了之後大於 1,說明系統是 round up,因此直接進位到最低位 (+2
)
+ 如果加減都回到 1 ,說明系統是 round to nearest? 此時如果最低位為 1 ,就進位。
90~91: 做完 rounding 之後,要把 mantissa 與 exponent 再次合併,所以要把 m
再 << 23
回去他的位置。要注意的是這裡的 m
並沒有 bias。所以 90 行要把 m
從 2 補數轉到 bias representation。
127=0b011111111
加回去。但q >> 1
事實上還有一個 1 在 bit 24,所以答案是 0x3f000000
。+ m << 23
的動作可能會影響到 float 的 sign bit,但是因為 m
不可能是 -128 所以一定會 overflow 把 float sign bit 變回 0。利用一樣的原理,我們一樣可以做整數開根號
int mySqrt(int x){
uint32_t pos, t=0, ans=0, s0=0;
if (x <= 1)
return x;
pos = 1 << ((31 - __builtin_clz(x)) & (~1));
while (pos) {
ans <<= 1;
t = s0 | pos;
if (t <= x) {
s0 = t + pos;
x -= t;
ans |= 1;
}
pos >>= 2;
s0 >>= 1;
}
return ans;
}
x = 10011
則 pos = 10000
,x = 110010
則 pos = 10000
。x
,而是 shift s0
與 pos
(原本的 r
)。因為當數字很大的時候,左移 x
有可能會 overflow。pos
一次移動兩格,也不方便透過 pos
來設定 bit 了。
int consecutiveNumbersSum(int N)
{
if (N < 1)
return 0;
int ret = 1;
int x = N;
for (int i = 2; i < x; i++) {
x += ZZZ;
if (x % i == 0)
ret++;
}
return ret;
}
從數學式
看,i
其實就對應到 x
裡面存放的就是 x=N
開始,每次減 x += 1-i
。檢查有沒有解就是等式右邊是否能被 x % i == 0
。