2025q1 Homework5 (assessment)

contributed by < gnkuan0712 >

一對一檢討

2025/05/09: bitwise, IEEE 754

題目

實作geomean，計算幾何平均數，不能用 FPU

int geomean(float *in, int n) { } // n 表示 in 的總數

補足背景知識

首先根據 IEEE 754 單精度浮點數會有 32 bits，其編碼方式為：

v a l u e = (- 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的三個部分取出方法

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

v a l u e = (- 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 語言: 浮點數運算中有提到浮點數的加法，若超過1也是把其中一個 mantissa 往右移，然後 exp 加一。
Image Not Showing Possible Reasons
- The image was uploaded to a note which you don't have access to
- The note which the image was originally uploaded to has been deleted
Learn More →

所以乘法應當也會有兩種情況

乘積大小	Bit47	代表值範圍	處理方式
≥ 2.0	1	[2.0, 4.0)	`P >>= 24; exp++`
[1.0,2.0)	0	[1.0, 2.0)	`P >>= 23;`

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 中 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），不用再做額外不同的操作

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
- 實作 queue.[ch]
- 實作 shuffle 命令
hw3
- 把畫面呈現這部份從 kernal space 移到 kernal space

因為自動飲料機而延畢的那一年：閱讀啟發

我看到了不斷在試錯，然後做實驗，透過數據來調整下一步方向的開發過程。

研讀第 1 到第 6 週「課程教材」和 CS:APP 3/e (至少到第二章)

第一章

專題的想法

參考2023、2024

回顧〈並行和多執行緒程式設計〉
- 這個我有興趣回顧
高效網頁伺服器
- 我想要研讀、撰寫document: 參考此同學
透過 Netfilter 自動過濾廣告
- 我覺得這個很酷

Backup

1 on 1

2的幾次方

#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;
}

兩數相乘

#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

#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);
}