# 2025q1 Homework5 (assessment)
contributed by < `gnkuan0712` >
## 一對一檢討
### 2025/05/09: bitwise, IEEE 754
#### 題目
實作`geomean`，計算幾何平均數，不能用 **FPU**
```c
int geomean(float *in, int n) { } // n 表示 in 的總數
```
#### 補足背景知識
首先根據 `IEEE 754` 單精度浮點數會有 32 bits，其編碼方式為：
$$value = (-1)^s \times M \times 2^E$$
* fraction(M)部分存放小數點後資料
* 舉例: `8.5 = 2³ + 1/2¹` 轉成二進制為 `1000.1` 調整成 `1.0001 x 2³`
* E的部分就是3，但exp編碼以`01111111 (127)`表示指數為0，所以指數為3時應先+3，表示成`10000010 (130)`
* frac(M) 的部分就是`0001`因此放到mantissa裡就是`00010000000000000000000` (共23bits)
* 把`uint32_t a`的三個部分取出方法
```c
uint32_t sign_a = a >> 31;
uint32_t exp_a = (a >> 23) & 0xFF; //0xFF = 1111 1111₂ 把以外的位元都清0
uint32_t frac_a = a & 0x7FFFFF;
```
接著兩個單精度浮點數相乘a*b
$$value = (-1)^{s_a+s_b} \times (M_a \times M_b) \times 2^{E_a+E_b}$$
* sign的部分其實就是把${s_a和s_b}$做XOR
* 接下來M的部分因為一開始在存放時，已省略1.0001最前面的1了，所以在做乘法時應當要補回 mantissa 最前面的1: `M_a = (1 << 23) | frac_a; ` 再做相乘
* exp的部分就是直接相加
除了以上大概念外，還會遇到像是 mantissa 相乘時會有 overflow 和正規化的問題的問題
* IEEE‑754 要求 mantissa 必須是 `1.xxxx` 的形式，也就是在 1.0~2.0 之間。如果乘出來 ≥ 2.0，就要把尾數右移一位，並把 exponent 加 1。
* 參考[你所不知道的 C 語言: 浮點數運算](https://hackmd.io/@sysprog/c-floating-point)中有提到浮點數的加法，若超過1也是把其中一個 mantissa 往右移，然後 exp 加一。
* 所以乘法應當也會有兩種情況
| 乘積大小 | Bit47 | 代表值範圍 | 處理方式 |
| ---------- | ----- | ----------- | ----------------- |
| ≥ 2.0 | 1 | \[2.0, 4.0) | `P >>= 24; exp++` |
| \[1.0,2.0) | 0 | \[1.0, 2.0) | `P >>= 23;` |
```c
uint64_t P = M_a * M_b;
// 正規化：若 ≥2.0，右移 24, exp+1；否則右移 23
if (P & (1ULL<<47)) { //假設Bit47是1則代表乘積超過2
P >>= 24;
exp_r += 1;
} else {
P >>= 23;
}
```
最後我認為若要擴展到多個浮點數相乘，應該要一邊做正規化，保持 mantissa 在 23 bits，但這樣會一直遇到**精準度遺失**的問題。
==TODO: 如何最大保留精準度==
原本上面的作法是無條件捨取，所以誤差會很顯著，接著參考了 [CMU: computer system](https://www.cs.cmu.edu/~213/lectures/m22/04-float.pdf) 中 page 30 有提到可以採用 round 的方式，也就是對 `mantissa` 做四捨五入。
四捨五入的實作用 bitwise 操作可以是 P 加上 P 的一半再往右移 n 個 bits `P = (P + (1ULL << (n - 1))) >> n;`
## 我的學習狀況
### 上課狀況
> 在第11周的隨堂測驗有提到「不管大小寫字母比較」的實做，這是我上課就聽的懂的程式碼
實作流程：
* 因為ASCII碼大小寫差32，所以只要往左移5，強迫它變成小寫
* ((ca - 'A') <= ('Z' - 'A'))) 看是不是大寫
* 是:回傳1
* 否:回傳0
* 接著就把1往左移5變成10000
* 然後再跟原本做or所以就會強制轉小寫
心得：
我認為這個寫法真的超酷的，尤其判斷大小寫不管是哪一個都是做一樣的操作（往左移五個bit），不用再做額外不同的操作
>
```c
static inline bool ci_ascii_eq(const char *a, const char *b, size_t n)
{
while (n--) {
unsigned char ca = *(const unsigned char *) a++;
unsigned char cb = *(const unsigned char *) b++;
ca |= ((unsigned char) ((ca - 'A') <= ('Z' - 'A'))) << 5;
cb |= ((unsigned char) ((cb - 'A') <= ('Z' - 'A'))) << 5;
if (ca != cb)
return false;
}
return true;
}
```
### 作業狀況
* [hw1](https://hackmd.io/ukV-iqKuS4e6-ysqDtpQLA?view)
* 實作 queue.[ch]
* 實作 shuffle 命令
* [hw3](https://hackmd.io/@kirua/linux2025-homework3)
* 把畫面呈現這部份從 kernal space 移到 kernal space
## [因為自動飲料機而延畢的那一年](https://www.opasschang.com/docs/the-story-of-auto-beverage-machine-1)：閱讀啟發
我看到了不斷在試錯，然後做實驗，透過數據來調整下一步方向的開發過程。
## 研讀第 1 到第 6 週「[課程教材](https://wiki.csie.ncku.edu.tw/sysprog/schedule)」和 [CS:APP 3/e](https://csapp.cs.cmu.edu/) (至少到第二章)
### 第一章
1. DMA
## 專題的想法
> 參考[2023](https://hackmd.io/@sysprog/linux2023-projects)、[2024](https://hackmd.io/@sysprog/linux2024-showcase)
* [回顧〈並行和多執行緒程式設計〉](https://hackmd.io/@sysprog/linux2023-projects#%E4%B8%A6%E8%A1%8C%E7%A8%8B%E5%BC%8F%E8%A8%AD%E8%A8%88)
* 這個我有興趣回顧
* [高效網頁伺服器](https://hackmd.io/@sysprog/linux2023-projects#%E9%AB%98%E6%95%88%E7%B6%B2%E9%A0%81%E4%BC%BA%E6%9C%8D%E5%99%A8)
* 我想要研讀、撰寫document: [參考此同學](https://hackmd.io/@sysprog/HJgX4_MH3)
* [透過 Netfilter 自動過濾廣告](https://hackmd.io/@sysprog/linux2023-projects#%E9%80%8F%E9%81%8E-Netfilter-%E8%87%AA%E5%8B%95%E9%81%8E%E6%BF%BE%E5%BB%A3%E5%91%8A)
* 我覺得這個很酷
### Backup
[1 on 1](https://hackmd.io/@sysprog/HkNNMsjaJl)
### 2的幾次方
```c
#include <stdio.h>
#include <stdint.h>
#include <string.h>
float my_ldexpf(float x, int n){
union{
float f;
uint32_t u32;
}val;
val.f = x;
uint32_t sign = val.u32 & 0x80000000;
uint32_t exp = (val.u32 >> 23) & 0xff;
uint32_t frac = val.u32 & 0x7fffff;
// x = +-inf or 0 or n=0
if (exp >= 0xff || exp == 0 || n ==0)
return x;
uint32_t new_exp = exp + n;
val.u32 = sign | new_exp << 23 | frac;
return val.f;
}
int main()
{
printf("%f", my_ldexpf(2.0, 2));
return 0;
}
```
### 兩數相乘
```c
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
// 把 float 轉成 uint32_t bits
static inline uint32_t float_to_bits(float f) {
uint32_t u;
memcpy(&u, &f, sizeof(u));
return u;
}
// 把 uint32_t bits 轉回 float
static inline float bits_to_float(uint32_t u) {
float f;
memcpy(&f, &u, sizeof(f));
return f;
}
// 純 bitwise 模擬單精度浮點乘法
float multiply(float a, float b) {
uint32_t A = float_to_bits(a);
uint32_t B = float_to_bits(b);
// 拆出 sign / exp / frac
uint32_t sign_a = A >> 31;
uint32_t exp_a = (A >> 23) & 0xFF;
uint32_t frc_a = A & 0x7FFFFF;
uint32_t sign_b = B >> 31;
uint32_t exp_b = (B >> 23) & 0xFF;
uint32_t frc_b = B & 0x7FFFFF;
// 結果 sign = xor
uint32_t sign_r = sign_a ^ sign_b;
// 重建 24‑bit mantissa（加回隱藏的 1）
uint64_t M_a = (1ULL<<23) | frc_a;
uint64_t M_b = (1ULL<<23) | frc_b;
// 整數乘法 → 48-bit
uint64_t P = M_a * M_b;
// 新 exponent = exp_a + exp_b - bias(127)
int32_t exp_r = (int32_t)exp_a + (int32_t)exp_b - 127;
// 正規化：若 ≥2.0，右移 24, exp+1；否則右移 23
if (P & (1ULL<<47)) {
P >>= 24;
exp_r += 1;
} else {
P >>= 23;
}
// 取出 23-bit fraction
uint32_t frc_r = (uint32_t)(P & 0x7FFFFF);
// 處理 under/overflow
if (exp_r <= 0) return sign_r ? -0.0f : 0.0f;
if (exp_r >= 255) return sign_r ? -INFINITY : INFINITY;
// 重組 bits
uint32_t R = (sign_r << 31)
| ((exp_r & 0xFF) << 23)
| frc_r;
return bits_to_float(R);
}
int main(void) {
float x, y, z;
// 測試 1
x = 2.0f; y = 3.5f;
z = multiply(x, y);
printf("%f * %f = %f\n", x, y, z);
return 0;
}
```
### GeoMean
```c
#include <stdint.h>
#include <string.h>
// Q16.16 固定小數常數
#define Q16_ONE (1<<16) // 1.0 in Q16.16
#define LOG2_K 94548 // ≈ (1/ln2)*2^16
#define EXP2_K 45426 // ≈ ln2*2^16
// 將 float bit 轉成 uint32
static inline uint32_t float_to_bits(float f) {
uint32_t u;
memcpy(&u, &f, sizeof(u));
return u;
}
// 將 uint32 bit 轉回 float
static inline float bits_to_float(uint32_t u) {
float f;
memcpy(&f, &u, sizeof(f));
return f;
}
int geomean(float *in, int n) {
int64_t sum_log2 = 0; // Q16.16 下的累加
for (int i = 0; i < n; i++) {
uint32_t v = float_to_bits(in[i]);
int32_t exp = ((v >> 23) & 0xFF) - 127; // 真實指數 E
uint32_t frac = v & 0x7FFFFF; // fraction bits
// 將 mantissa 轉成 Q16.16： f = mant/2^23 → mant>>7
int32_t f_fp = (int32_t)(frac >> 7);
// log2(x) ≈ E + f*(1/ln2)
int32_t log2i = (exp << 16)
+ (int32_t)(((int64_t)f_fp * LOG2_K) >> 16);
sum_log2 += log2i;
}
// 平均後的 log2 (Q16.16)
int32_t avg_log2 = (int32_t)(sum_log2 / n);
// 分離整數、小數部分
int32_t I = avg_log2 >> 16;
int32_t f_fp = avg_log2 & 0xFFFF;
// exp2(f) ≈ 1 + f*ln2 → Q16.16
int32_t mant_fp = Q16_ONE + (int32_t)(((int64_t)f_fp * EXP2_K) >> 16);
// 重組 IEEE‑754 bits
int32_t exp_out = I + 127;
if (exp_out <= 0) return 0; // underflow → 0
if (exp_out >= 255) return 0x7F800000; // overflow → +∞
uint32_t frac_out = (uint32_t)(((int64_t)(mant_fp - Q16_ONE) * (1<<23)) >> 16) & 0x7FFFFF;
uint32_t result = ((uint32_t)exp_out << 23) | frac_out;
// 轉回 float 並四捨五入成 int
float gf = bits_to_float(result);
return (int)(gf + 0.5f);
}
int main(void) {
float data1[] = {2.0f, 8.0f, 4.0f};
int gm1 = geomean(data1, 3);
printf("Test 1: geomean({2,8,4}) = %d (預期約 4)\n", gm1);
}
```
Sign in with Wallet
Connect another wallet