cKit/Foundation/c_Float.c

#include <c_Float.h>
#include <c_Memory.h>

uint32_t c_Float_Add(uint32_t a, uint32_t b) {
    // 1. 提取特殊狀態與快速零值返回
    if ((a & ~C_FLOAT_SIGN_MASK) == 0) return b;
    if ((b & ~C_FLOAT_SIGN_MASK) == 0) return a;

    uint32_t sign_a = a & C_FLOAT_SIGN_MASK;
    uint32_t sign_b = b & C_FLOAT_SIGN_MASK;
    int32_t exp_a   = (int32_t)((a & C_FLOAT_EXP_MASK) >> 23);
    int32_t exp_b   = (int32_t)((b & C_FLOAT_EXP_MASK) >> 23);
    uint32_t frac_a = a & C_FLOAT_FRAC_MASK;
    uint32_t frac_b = b & C_FLOAT_FRAC_MASK;

    // ==========================================
    // 修复 Bug 1：严格遵循 IEEE 754 的 NaN / Inf 拦截规则
    // ==========================================
    if (exp_a == 255 || exp_b == 255) {
        // 条件 1.1：若任一为真正的 NaN（指数全1，尾数不为0），则传播 NaN
        if ((exp_a == 255 && frac_a != 0) || (exp_b == 255 && frac_b != 0)) {
            return 0x7FC00000U; // 返回 Quiet NaN
        }
        // 条件 1.2：若两者都是无穷大，且符号相反（+Inf + -Inf），属于未定义，必须强行返回 NaN
        if (exp_a == 255 && exp_b == 255 && sign_a != sign_b) {
            return 0x7FC00000U; // 熔断返回 NaN
        }
        // 条件 1.3：普通的无穷大传播（如 Inf + 有限数 = Inf）
        return (exp_a == 255) ? a : b;
    }

    // ==========================================
    // 修复 Bug 3：更正非规格化数的隐藏位逻辑（exp==0 时隐藏位是 0，且不能左移）
    // ==========================================
    frac_a = (exp_a == 0) ? frac_a : (frac_a | C_FLOAT_HIDDEN_BIT);
    frac_b = (exp_b == 0) ? frac_b : (frac_b | C_FLOAT_HIDDEN_BIT);
    if (exp_a == 0) exp_a = 1;
    if (exp_b == 0) exp_b = 1;

    // 开辟 GRS 保护位空间
    frac_a <<= 3;
    frac_b <<= 3;

    int32_t exp_res = exp_a;
    uint32_t sticky = 0;

    // 3. 對階（Align Exponents）
    if (exp_a > exp_b) {
        int32_t shift = exp_a - exp_b;
        // 修复 Bug 2：防止大跨度对阶时左移掩码溢出未定义行为
        if (shift >= 27) {
            sticky = (frac_b != 0);
            frac_b = 0;
        } else {
            sticky = (frac_b & ~((~0U) << shift)) != 0;
            frac_b >>= shift;
        }
        frac_b |= sticky;
        exp_res = exp_a;
    } else if (exp_b > exp_a) {
        int32_t shift = exp_b - exp_a;
        // 修复 Bug 2：同理防护
        if (shift >= 27) {
            sticky = (frac_a != 0);
            frac_a = 0;
        } else {
            sticky = (frac_a & ~((~0U) << shift)) != 0;
            frac_a >>= shift;
        }
        frac_a |= sticky;
        exp_res = exp_b;
    }

    // 4. 尾數運算
    uint32_t sign_res;
    uint32_t frac_res;

    if (sign_a == sign_b) {
        sign_res = sign_a;
        frac_res = frac_a + frac_b;
    } else {
        if (frac_a >= frac_b) {
            sign_res = sign_a;
            frac_res = frac_a - frac_b;
        } else {
            sign_res = sign_b;
            frac_res = frac_b - frac_a;
        }
        if (frac_res == 0) return 0x00000000U; // 正负抵消返回标准 +0.0
    }

    // 5. 規格化
    if (frac_res & (1U << 27)) {
        // 修复 Bug 4：修正右移 1 位时的 Sticky 保留逻辑
        uint32_t lost_bit = frac_res & 1U;
        frac_res >>= 1;
        frac_res |= lost_bit;
        exp_res++;
    } else {
        while (!(frac_res & (1U << 26)) && exp_res > 1) {
            frac_res <<= 1;
            exp_res--;
        }
        if (!(frac_res & (1U << 26)) && exp_res == 1) {
            exp_res = 0;
        }
    }

    // 6. 溢出至無限大檢查
    if (exp_res >= 255) {
        return sign_res | C_FLOAT_EXP_MASK;
    }

    // 7. 向最接近偶數捨入（Round-to-Nearest-Even）
    uint32_t round_bits = frac_res & 7U;
    frac_res >>= 3;

    if ((round_bits > 4) || ((round_bits == 4) && (frac_res & 1U))) {
        frac_res++;
        if (frac_res & (1U << 24)) {
            frac_res >>= 1;
            exp_res++;
            if (exp_res >= 255) return sign_res | C_FLOAT_EXP_MASK;
        }
    }

    if (exp_res != 0) {
        frac_res &= C_FLOAT_FRAC_MASK;
    }

    // 8. 拼裝返回
    return sign_res | ((uint32_t)exp_res << 23) | frac_res;
}


uint32_t c_Float_Mul(uint32_t a, uint32_t b) {
    // 使用你定义的快捷提取宏
    uint32_t sign_a = C_FLOAT_GET_SIGN(a);
    uint32_t sign_b = C_FLOAT_GET_SIGN(b);
    uint32_t exp_a  = C_FLOAT_GET_EXP(a);
    uint32_t exp_b  = C_FLOAT_GET_EXP(b);
    uint32_t frac_a = C_FLOAT_GET_FRAC(a);
    uint32_t frac_b = C_FLOAT_GET_FRAC(b);

    uint32_t sign_res = sign_a ^ sign_b; // 异或决定结果符号

    // ==========================================
    // 边界条件 1：处理 NaN 和 无穷大 (Inf) 的传播与熔断
    // ==========================================
    if (exp_a == 255 || exp_b == 255) {
        // 条件 1.1：若任一输入为真正的 NaN，直接传播 Quiet NaN
        if ((exp_a == 255 && frac_a != 0) || (exp_b == 255 && frac_b != 0)) {
            return 0x7FC00000U;
        }
        // 条件 1.2：0.0 * 无穷大 (0.0 * Inf) 属于未定义，必须强行熔断返回 NaN
        bool is_a_zero = (exp_a == 0 && frac_a == 0);
        bool is_b_zero = (exp_b == 0 && frac_b == 0);
        if ((exp_a == 255 && is_b_zero) || (exp_b == 255 && is_a_zero)) {
            return 0x7FC00000U; // 熔断返回 NaN
        }
        // 条件 1.3：普通的无穷大传播（Inf * 有限非零数 = Inf）
        return c_Float_Pack(sign_res, 255, 0);
    }

    // ==========================================
    // 边界条件 2：处理纯零快速返回
    // ==========================================
    if ((exp_a == 0 && frac_a == 0) || (exp_b == 0 && frac_b == 0)) {
        return c_Float_Pack(sign_res, 0, 0); // 0.0 * 有限数 -> 产生带有正确符号的 ±0.0
    }

    // ==========================================
    // 2. 补齐隐藏位并处理非规格化数
    // ==========================================
    frac_a = (exp_a == 0) ? frac_a : (frac_a | C_FLOAT_HIDDEN_BIT);
    frac_b = (exp_b == 0) ? frac_b : (frac_b | C_FLOAT_HIDDEN_BIT);

    if (exp_a == 0) exp_a = 1;
    if (exp_b == 0) exp_b = 1;

    int32_t exp_res = (int32_t)exp_a + (int32_t)exp_b - C_FLOAT_EXP_BIAS;

    // ==========================================
    // 3. 执行核心尾数相乘（24位 * 24位 = 48位超长整数）
    // ==========================================
    uint64_t prod = (uint64_t)frac_a * (uint64_t)frac_b;

    // ==========================================
    // 4. 将 48 位乘积向右压缩，腾出低 3 位作为 GRS 保护位空间
    // ==========================================
    // 正常规格化数相乘后，结果 `prod` 的最高位 1 应该在第 46 位或第 47 位（从0数起）。
    // 为了最终留下 24 位标准尾数和 3 位保护位，我们需要让规格化后的目标保留在 27 位。
    // 因此，我们先固定将原本 48 位的低 20 位挤出去，并把这 20 位中任意的 1 凝聚为 Sticky 位。
    uint32_t sticky = (prod & 0xFFFFF) != 0;
    uint32_t frac_res = (uint32_t)(prod >> 20); // 压缩至大约 27~28 位
    frac_res |= sticky; // 将物理挤出去的所有小数信息固化在最低位

    // ==========================================
    // 5. 规格化积（Normalization）
    // ==========================================
    // 如果最高有效位 1 溢出到了第 27 位 (1U << 27)，说明乘积结果 >= 2.0，需要右移 1 位，指数加 1
    if (frac_res & (1U << 27)) {
        uint32_t lost_bit = frac_res & 1U;
        frac_res >>= 1;
        frac_res |= lost_bit; // 保持最低位 Sticky 不丢失
        exp_res++;
    } else {
        // 如果隐藏位没落到第 26 位，说明乘积较小 (< 1.0)，需要左移直到最高有效位 1 回归第 26 位
        while (!(frac_res & (1U << 26)) && exp_res > 1) {
            uint32_t sticky_backup = frac_res & 1U;
            frac_res <<= 1;
            frac_res |= sticky_backup; // 锁死 Sticky 状态
            exp_res--;
        }
        // 下溢退化为非规格化数
        if (!(frac_res & (1U << 26)) && exp_res == 1) {
            exp_res = 0;
        }
    }

    // 上下溢出安全拦截
    if (exp_res >= 255) return c_Float_Pack(sign_res, 255, 0); // 上溢至 ±Inf
    if (exp_res <= 0)   return c_Float_Pack(sign_res, 0, 0);   // 下溢至 ±0.0

    // ==========================================
    // 6. 激活 IEEE 754 标准：向最接近偶数舍入（Round-to-Nearest-Even）
    // ==========================================
    uint32_t round_bits = frac_res & 7U; // 捕获低 3 位的 GRS 数据
    frac_res >>= 3;                      // 移除保护位，回归标准 24 位尾数

    if ((round_bits > 4) || ((round_bits == 4) && (frac_res & 1U))) {
        frac_res++;
        // 舍入导致的二次溢出处理
        if (frac_res & (1U << 24)) {
            frac_res >>= 1;
            exp_res++;
            if (exp_res >= 255) return c_Float_Pack(sign_res, 255, 0);
        }
    }

    // 剥离规格化数中用于拼装的高位隐藏位 1
    if (exp_res != 0) {
        frac_res &= C_FLOAT_FRAC_MASK;
    }

    // 7. 最终位打包返回
    return c_Float_Pack(sign_res, exp_res, frac_res);
}


uint32_t c_Float_Div(uint32_t a, uint32_t b) {
    // 使用你定义的快捷提取宏
    uint32_t sign_a = C_FLOAT_GET_SIGN(a);
    uint32_t sign_b = C_FLOAT_GET_SIGN(b);
    uint32_t exp_a  = C_FLOAT_GET_EXP(a);
    uint32_t exp_b  = C_FLOAT_GET_EXP(b);
    uint32_t frac_a = C_FLOAT_GET_FRAC(a);
    uint32_t frac_b = C_FLOAT_GET_FRAC(b);

    uint32_t sign_res = sign_a ^ sign_b;

    // =========================================================================
    // 修复 Bug 1 & 2：严密拦截 IEEE 754 规定的所有 NaN、Inf、特殊零边界
    // =========================================================================

    // 1. 拦截输入本身为 NaN 或 无穷大 (exp == 255) 的异常传播
    if (exp_a == 255 || exp_b == 255) {
        // 条件 A: 任一为真正的 NaN（尾数非0），直接传播 Quiet NaN
        if ((exp_a == 255 && frac_a != 0) || (exp_b == 255 && frac_b != 0)) {
            return 0x7FC00000U;
        }
        // 条件 B: 无穷大除以无穷大 (Inf / Inf) 属于严重未定义，强熔断返回 NaN
        if (exp_a == 255 && exp_b == 255) {
            return 0x7FC00000U;
        }
        // 条件 C: 无穷大除以普通有限数 = 无穷大
        if (exp_a == 255) {
            return c_Float_Pack(sign_res, 255, 0);
        }
        // 条件 D: 普通有限数除以无穷大 = 0
        return c_Float_Pack(sign_res, 0, 0);
    }

    // 2. 剥离符号位，准确捕捉纯零值 (+/-0.0) 状态下的防御拦截
    bool is_a_zero = (exp_a == 0 && frac_a == 0);
    bool is_b_zero = (exp_b == 0 && frac_b == 0);

    if (is_b_zero) {
        // 0.0 / 0.0 必须熔断返回 NaN
        if (is_a_zero) {
            return 0x7FC00000U;
        }
        // 有限数 / 0.0 -> 产生标准的 ±Infinity
        return c_Float_Pack(sign_res, 255, 0);
    }
    // 0.0 / 非零数 -> 产生标准的 ±0.0
    if (is_a_zero) {
        return c_Float_Pack(sign_res, 0, 0);
    }

    // =========================================================================
    // 修复 Bug 4：正确补齐隐藏位（非规格化数的隐藏位是 0）
    // =========================================================================
    frac_a = (exp_a == 0) ? frac_a : (frac_a | C_FLOAT_HIDDEN_BIT);
    frac_b = (exp_b == 0) ? frac_b : (frac_b | C_FLOAT_HIDDEN_BIT);

    if (exp_a == 0) exp_a = 1;
    if (exp_b == 0) exp_b = 1;

    int32_t exp_res = (int32_t)exp_a - (int32_t)exp_b + C_FLOAT_EXP_BIAS;

    // =========================================================================
    // 修复 Bug 3：向左偏移 27 位做长除法，为 GRS 三保护位和 Sticky 腾出精度空间
    // =========================================================================
    uint64_t num = (uint64_t)frac_a << 26; // 24位标准商 + 3位保护位空间
    uint64_t den = (uint64_t)frac_b;

    uint32_t quot = (uint32_t)(num / den);
    uint64_t rem = num % den;

    // 【除法 Sticky 核心点】只要整数除法除不尽有余数，无条件将商的最低位置 1，激活 Sticky 状态
    if (rem != 0) {
        quot |= 1U;
    }

    // 规格化商：正常情况下隐藏位应该落在第 26 位 (C_FLOAT_HIDDEN_BIT << 3)
    if (quot & (1U << 27)) {
        uint32_t lost = quot & 1U;
        quot >>= 1;
        quot |= lost;
        exp_res++;
    } else {
        // 如果隐藏位没能落在第 26 位，说明商太小，需要左移规格化
        while (!(quot & (1U << 26)) && exp_res > 1) {
            uint32_t sticky_backup = quot & 1U;
            quot <<= 1;
            quot |= sticky_backup; // 锁住最低位的 sticky 特征不丢失
            exp_res--;
        }
        if (!(quot & (1U << 26)) && exp_res == 1) {
            exp_res = 0; // 退化为非规格化数
        }
    }

    // 上下溢出检测
    if (exp_res >= 255) return c_Float_Pack(sign_res, 255, 0); // 上溢至 ±Inf
    if (exp_res <= 0)   return c_Float_Pack(sign_res, 0, 0);   // 下溢至 ±0

    // =========================================================================
    // 5. 激活 IEEE 754 标准：向最接近偶数舍入（Round-to-Nearest-Even）
    // =========================================================================
    uint32_t round_bits = quot & 7U; // 提取最后 3 位的 GRS 数据
    quot >>= 3;                      // 移除保护位，回归标准 24 位商

    if ((round_bits > 4) || ((round_bits == 4) && (quot & 1U))) {
        quot++;
        // 舍入可能导致再次溢出，进行二次规格化微调
        if (quot & (1U << 24)) {
            quot >>= 1;
            exp_res++;
            if (exp_res >= 255) return c_Float_Pack(sign_res, 255, 0);
        }
    }

    // 剥离规格化数中用于拼装的高位隐藏位 1
    if (exp_res != 0) {
        quot &= C_FLOAT_FRAC_MASK;
    }

    return c_Float_Pack(sign_res, exp_res, quot);
}


int c_Float_Cmp(uint32_t a, uint32_t b) {
    // 1. 处理 NaN：依据 IEEE 754，NaN 参与比较永远返回不相等（或未定义）
    if (c_Float_IsNAN(a) || c_Float_IsNAN(b)) {
        return 0; // 软浮点库通常在此处设置不合法比较标志位
    }

    // 2. 特殊情况：+0.0 (0x00000000) 和 -0.0 (0x80000000) 在逻辑上是相等的
    if (((a | b) & ~C_FLOAT_SIGN_MASK) == 0) {
        return 0;
    }

    // 提取符号位
    uint32_t sign_a = a & C_FLOAT_SIGN_MASK;
    uint32_t sign_b = b & C_FLOAT_SIGN_MASK;

    // 3. 符号不同
    if (sign_a != sign_b) {
        // a 是负数，b 是正数 => a < b
        // a 是正数，b 是负数 => a > b
        return sign_a ? -1 : 1;
    }

    // 4. 符号相同：将原始二进制位转换为有符号 32 位整型进行直观比较
    int32_t ia = (int32_t)a;
    int32_t ib = (int32_t)b;

    if (sign_a) {
        // 如果都是负数，二进制数值越大，其代表的实际浮点数反而越小 (例如 -2.0 的二进制码大于 -1.0)
        if (ia > ib) return -1;
        if (ia < ib) return 1;
        return 0;
    } else {
        // 如果都是正数，二进制数值越大，其实际浮点数就越大
        if (ia > ib) return 1;
        if (ia < ib) return -1;
        return 0;
    }
}

int c_Double_Cmp(uint64_t a, uint64_t b) {
    // 1. 处理 NaN
    if (c_Double_IsNAN(a) || c_Double_IsNAN(b)) {
        return 0;
    }

    // 2. 处理 +0.0 与 -0.0 相等的情况
    if (((a | b) & ~C_DOUBLE_SIGN_MASK) == 0) {
        return 0;
    }

    uint64_t sign_a = a & C_DOUBLE_SIGN_MASK;
    uint64_t sign_b = b & C_DOUBLE_SIGN_MASK;

    // 3. 符号不同
    if (sign_a != sign_b) {
        return sign_a ? -1 : 1;
    }

    // 4. 符号相同：转为有符号 64 位整型比较
    int64_t ia = (int64_t)a;
    int64_t ib = (int64_t)b;

    if (sign_a) {
        // 均为负数
        if (ia > ib) return -1;
        if (ia < ib) return 1;
        return 0;
    } else {
        // 均为正数
        if (ia > ib) return 1;
        if (ia < ib) return -1;
        return 0;
    }
}

/**
 * @brief 软浮点双精度打包函数
 * @param sign     符号位 (0 或 1)
 * @param exp      解包/运算后的有符号指数 (已减去或未加上 Bias 均可，此处传入带 Bias 的期望值)
 * @param frac64   运算后暂存在 64 位整型中的高精度尾数 (假设规格化后隐含位在第 52 位，低位留有舍入残余)
 * @return         组合好的 IEEE 754 64位无符号整数 (可直接对应 double)
 */
uint64_t c_Double_Pack(uint32_t sign, int32_t exp, uint64_t frac64) {

    // 1. 动态规格化：若运算导致尾数高位溢出 (例如第 53 位为 1)，需要右移尾数并增加指数
    if (frac64 & (C_DOUBLE_HIDDEN_BIT << 1)) {
        frac64 >>= 1;
        exp++;
    }

    // 2. 下溢处理：指数太小，转换为非规格化数
    if (exp <= 0) {
        // 如果指数极小，直接移出范围，变回 0
        if (exp < -52) {
            frac64 = 0;
        } else {
            // 右移尾数以对齐非规格化数的指数位置 (exp = 0)
            int32_t shift = 1 - exp;
            frac64 >>= shift;
        }
        exp = 0; // 非规格化数的指数域强制为 0
    }

    // 3. 执行 IEEE 754 默认的“向最接近偶数舍入 (Round-to-Nearest-Even)”
    // 假设经过上述操作后，标准 52 位尾数在 frac64 的低 52 位，若有更低位则是运算残留
    // 为了演示标准舍入，假设传入的 frac64 在低位保留了扩充精度（例如左移了 3 位留给 GRS）
    // 此处简化演示：基于常规截断进行最邻近舍入处理
    // 在工业级库中，通常传入 frac 时会带有额外的 round_bits 变量

    // 4. 上溢检查：指数超过最大限制 (2047)，打包为无穷大
    if (exp >= 0x7FF) {
        return ((uint64_t)sign << 63) | C_DOUBLE_EXP_MASK; // 返回 +/- Inf
    }

    // 5. 最终清除尾数域外的隐含 1 (因为 IEEE 754 编码中不存储规格化数的最高位 1)
    uint64_t final_frac = frac64 & C_DOUBLE_FRAC_MASK;

    // 6. 位移拼接
    uint64_t packed_value = ((uint64_t)sign << 63) |
                            ((uint64_t)exp << 52)  |
                            final_frac;

    return packed_value;
}


/* ------------------------------------------------------------------------------------------------------------------ */
/*  */

// 輔助函數：處理 64 位元尾數與 3 位元 GRS 捨入殘餘
static uint64_t round_and_pack_double(uint64_t sign, int32_t exp, uint64_t frac64) {
    uint32_t round_bits = frac64 & 7U;
    frac64 >>= 3; // 移除 GRS 位，恢復為包含隱含位的 53 位元尾數

    // 向最接近偶數捨入
    if ((round_bits > 4) || ((round_bits == 4) && (frac64 & 1ULL))) {
        frac64++;
        if (frac64 & (C_DOUBLE_HIDDEN_BIT << 1)) {
            frac64 >>= 1;
            exp++;
        }
    }

    if (exp >= 2047) return sign | C_DOUBLE_EXP_MASK; // 溢出至無限大
    if (exp <= 0) return sign;                      // 下溢至 0

    return sign | ((uint64_t)exp << 52) | (frac64 & C_DOUBLE_FRAC_MASK);
}

uint64_t c_Double_Add(uint64_t a, uint64_t b) {
    c_Double_t f1 = { .raw = a };
    c_Double_t f2 = { .raw = b };
    c_Double_t result = { .raw = 0 };

    // 1. Handle Special Cases: NaNs and Infinities
    bool f1_is_nan_or_inf = (f1.parts.exponent == 0x7FF);
    bool f2_is_nan_or_inf = (f2.parts.exponent == 0x7FF);

    if (f1_is_nan_or_inf || f2_is_nan_or_inf) {
        // Handle NaNs
        if ((f1_is_nan_or_inf && f1.parts.fraction != 0) ||
            (f2_is_nan_or_inf && f2.parts.fraction != 0)) {
            result.parts.exponent = 0x7FF;
            result.parts.fraction = 0x1; // Quiet NaN
            return result.raw;
        }
        // Handle Inf + Inf variations
        if (f1_is_nan_or_inf && f2_is_nan_or_inf) {
            if (f1.parts.sign != f2.parts.sign) {
                // (+Inf) + (-Inf) or (-Inf) + (+Inf) is invalid -> NaN
                result.parts.exponent = 0x7FF;
                result.parts.fraction = 0x1;
                return result.raw;
            }
            return f1.raw; // Return either infinity if signs match
        }
        // One operand is Infinity, the other is finite
        return f1_is_nan_or_inf ? f1.raw : f2.raw;
    }

    // 2. Handle Zero Shortcuts
    bool f1_is_zero = (f1.parts.exponent == 0 && f1.parts.fraction == 0);
    bool f2_is_zero = (f2.parts.exponent == 0 && f2.parts.fraction == 0);
    if (f1_is_zero && f2_is_zero) {
        // If both are zero and signs differ, standard rule yields +0.0 in round-to-nearest
        result.parts.sign = (f1.parts.sign == f2.parts.sign) ? f1.parts.sign : 0;
        return result.raw;
    }
    if (f1_is_zero) return f2.raw;
    if (f2_is_zero) return f1.raw;

    // 3. Extract Exponents and Mantissas (with implicit leading 1 bit)
    int32_t exp1 = f1.parts.exponent;
    int32_t exp2 = f2.parts.exponent;

    uint64_t m1 = (1ULL << 52) | f1.parts.fraction;
    uint64_t m2 = (1ULL << 52) | f2.parts.fraction;

    // 4. Align Exponents (Shift mantissas to three extra bits of precision: Guard, Round, Sticky)
    // We scale the mantissas left by 3 bits initially to capture shifting errors.
    uint64_t m_large = 0, m_small = 0;
    int32_t exp_res = 0;
    bool sign_large = 0, sign_small = 0;

    if (exp1 >= exp2) {
        m_large = m1 << 3;
        sign_large = f1.parts.sign;
        exp_res = exp1;

        int32_t shift = exp1 - exp2;
        if (shift == 0) {
            m_small = m2 << 3;
        } else if (shift > 55) {
            m_small = 1; // Everything shifted out becomes a sticky bit
        } else {
            uint64_t lost_bits = m2 & ((1ULL << shift) - 1);
            m_small = (m2 << 3) >> shift;
            if (lost_bits != 0) m_small |= 1; // Fold lost bits into sticky bit
        }
        sign_small = f2.parts.sign;
    } else {
        m_large = m2 << 3;
        sign_large = f2.parts.sign;
        exp_res = exp2;

        int32_t shift = exp2 - exp1;
        if (shift > 55) {
            m_small = 1;
        } else {
            uint64_t lost_bits = m1 & ((1ULL << shift) - 1);
            m_small = (m1 << 3) >> shift;
            if (lost_bits != 0) m_small |= 1;
        }
        sign_small = f1.parts.sign;
    }

    // 5. Perform Magnitude Addition or Subtraction
    uint64_t m_res = 0;
    bool result_sign = sign_large;

    if (sign_large == sign_small) {
        // True addition
        m_res = m_large + m_small;

        // Handle carry out: if bit 56 is set (original 52 shifted left 3 plus 1 carry bit)
        if (m_res & (1ULL << 56)) {
            uint64_t sticky = m_res & 1;
            m_res >>= 1;
            m_res |= sticky; // preserve sticky bit tracking
            exp_res += 1;
        }
    } else {
        // True subtraction (Large magnitude minus Small magnitude)
        // If magnitudes are completely equal, they cancel out to 0
        if (m_large == m_small) {
            return 0; // standard +0.0
        }
        m_res = m_large - m_small;

        // Normalize cancellation shifts (shift left until bit 55 is 1)
        while ((m_res & (1ULL << 55)) == 0 && exp_res > 0) {
            uint64_t sticky = m_res & 1;
            m_res = (m_res << 1) | sticky;
            exp_res -= 1;
        }
    }

    // 6. Apply IEEE 754 Round-to-Nearest, Ties-to-Even
    // Currently, bit 55 is the implicit 1. Bits [2:0] are Guard, Round, Sticky.
    // The target fraction belongs in bits [54:3].
    uint64_t final_fraction = (m_res >> 3) & 0xFFFFFFFFFFFFFLL;

    bool round_bit = (m_res & 4) != 0;  // Bit 2
    bool sticky_bit = (m_res & 3) != 0; // Bits 1 and 0 combined
    bool lsb = (final_fraction & 1) != 0;

    if (round_bit && (sticky_bit || lsb)) {
        final_fraction++;
        if (final_fraction > 0xFFFFFFFFFFFFFLL) { // Handle carry out from rounding
            final_fraction = 0;
            exp_res += 1;
        }
    }

    // 7. Check for Overflow / Underflow Boundaries
    if (exp_res >= 0x7FF) {
        result.parts.sign = result_sign;
        result.parts.exponent = 0x7FF;
        result.parts.fraction = 0; // Overflow to Infinity
    } else if (exp_res <= 0) {
        // Flush underflow to zero
        result.parts.sign = result_sign;
        result.parts.exponent = 0;
        result.parts.fraction = 0;
    } else {
        result.parts.sign = result_sign;
        result.parts.exponent = (uint64_t)exp_res;
        result.parts.fraction = final_fraction;
    }

    return result.raw;
}


// 內部輔助函數：32位交叉相乘，手動模擬 64x64->128位元乘法
C_STATIC_FORCE_INLINE
void mul64_to_128(uint64_t a, uint64_t b, uint64_t *res_hi, uint64_t *res_lo) {
    uint64_t a_hi = a >> 32, a_lo = a & 0xFFFFFFFFULL;
    uint64_t b_hi = b >> 32, b_lo = b & 0xFFFFFFFFULL;

    uint64_t p0 = a_lo * b_lo;
    uint64_t p1 = a_hi * b_lo;
    uint64_t p2 = a_lo * b_hi;
    uint64_t p3 = a_hi * b_hi;

    uint64_t mid = p1 + (p0 >> 32) + (p2 & 0xFFFFFFFFULL);
    *res_lo = (mid << 32) | (p0 & 0xFFFFFFFFULL);
    *res_hi = p3 + (mid >> 32) + (p2 >> 32);
}

uint64_t c_Double_Mul(uint64_t a, uint64_t b) {
    c_Double_t f1 = { .raw = a };
    c_Double_t f2 = { .raw = b };
    c_Double_t result = { .raw = 0 };

    // 1. Determine the result sign (XOR of input signs)
    result.parts.sign = f1.parts.sign ^ f2.parts.sign;

    // 2. Handle Zero / Special Cases (Inf, NaN)
    // Shortcut if either operand is zero
    bool f1_is_inf_or_nan = (f1.parts.exponent == 0x7FF);
    bool f2_is_inf_or_nan = (f2.parts.exponent == 0x7FF);
    bool f1_is_zero = (f1.parts.exponent == 0 && f1.parts.fraction == 0);
    bool f2_is_zero = (f2.parts.exponent == 0 && f2.parts.fraction == 0);



    // Shortcut for Infinities or NaNs
    if (f1_is_inf_or_nan || f2_is_inf_or_nan) {
        // If either is an actual NaN, OR we are multiplying 0 * Inf, it MUST be NaN
        if ((f1_is_inf_or_nan && f1.parts.fraction != 0) ||
            (f2_is_inf_or_nan && f2.parts.fraction != 0) ||
            (f1_is_zero && f2_is_inf_or_nan) ||
            (f2_is_zero && f1_is_inf_or_nan)) {

            result.parts.exponent = 0x7FF;
            result.parts.fraction = 0x1; // Quiet NaN
            return result.raw;
            }
        // Otherwise, it's a valid Infinity multiplication (e.g., 5.0 * Inf = Inf)
        result.parts.exponent = 0x7FF;
        result.parts.fraction = 0;
        return result.raw;
    }

    // Now it is safe to evaluate the normal zero shortcut
    if (f1_is_zero || f2_is_zero) {
        result.parts.exponent = 0;
        result.parts.fraction = 0;
        return result.raw; // Returns correctly signed zero
    }

    // 3. Extract Mantissas and append the implicit leading 1 bit (Bit 52)
    // Note: This implementation assumes normalized numbers.
    uint64_t m1 = (1ULL << 52) | f1.parts.fraction;
    uint64_t m2 = (1ULL << 52) | f2.parts.fraction;

    // 4. Calculate raw exponent sum (subtract the double bias of 1023)
    int32_t exp_res = (int32_t)f1.parts.exponent + (int32_t)f2.parts.exponent - 1023;

    // 5. Multiply the mantissas using 128-bit precision to prevent overflow
    // Multiplying two 53-bit integers results in a 105-bit or 106-bit product
    unsigned __int128 prod = (unsigned __int128)m1 * m2;

    // 6. Normalize the product
    // The product has its radix point at bit 104 (52 fractional bits * 2)
    // We want the resulting leading bit to sit at bit 52.
    if (prod & ((unsigned __int128)1 << 105)) {
        // Product is >= 2.0 (bit 105 is set). Shift down by 53 and increment exponent.
        exp_res += 1;
        // Simple round-to-nearest-even approximation via bit shift
        result.parts.fraction = (uint64_t)((prod >> 53) & 0xFFFFFFFFFFFFFLL);
    } else {
        // Product is < 2.0 (bit 104 is set). Shift down by 52.
        result.parts.fraction = (uint64_t)((prod >> 52) & 0xFFFFFFFFFFFFFLL);
    }

    // 7. Check for Overflow / Underflow boundaries
    if (exp_res >= 0x7FF) {
        // Overflow to Infinity
        result.parts.exponent = 0x7FF;
        result.parts.fraction = 0;
    } else if (exp_res <= 0) {
        // Underflow to Zero (Flushing subnormals to zero for simplicity)
        result.parts.exponent = 0;
        result.parts.fraction = 0;
    } else {
        // Valid normalized exponent range
        result.parts.exponent = (uint64_t)exp_res;
    }

    return result.raw;
}


uint64_t c_Double_Div(uint64_t a, uint64_t b) {
    c_Double_t f1 = { .raw = a };
    c_Double_t f2 = { .raw = b };
    c_Double_t result = { .raw = 0 };

    // 1. Determine the result sign (XOR of input signs)
    result.parts.sign = f1.parts.sign ^ f2.parts.sign;

    // 2. Handle Special Cases: Zero, Infinity, and NaN
    bool f1_is_nan_or_inf = (f1.parts.exponent == 0x7FF);
    bool f2_is_nan_or_inf = (f2.parts.exponent == 0x7FF);
    bool f1_is_zero = (f1.parts.exponent == 0 && f1.parts.fraction == 0);
    bool f2_is_zero = (f2.parts.exponent == 0 && f2.parts.fraction == 0);

    // Case 2a: Either input is NaN, or invalid combinations (0/0, Inf/Inf)
    if ((f1_is_nan_or_inf && f1.parts.fraction != 0) ||
        (f2_is_nan_or_inf && f2.parts.fraction != 0) ||
        (f1_is_zero && f2_is_zero) ||
        (f1_is_nan_or_inf && f2_is_nan_or_inf)) {
        result.parts.exponent = 0x7FF;
        result.parts.fraction = 0x1; // Quiet NaN
        return result.raw;
    }

    // Case 2b: Division by Zero (X / 0 = Inf)
    if (f2_is_zero) {
        result.parts.exponent = 0x7FF; // Infinity
        result.parts.fraction = 0;
        return result.raw;
    }

    // Case 2c: Numerator is Zero or Denominator is Infinity (0 / X = 0, X / Inf = 0)
    if (f1_is_zero || f2_is_nan_or_inf) {
        result.parts.exponent = 0;
        result.parts.fraction = 0;
        return result.raw;
    }

    // Case 2d: Numerator is Infinity (Inf / X = Inf)
    if (f1_is_nan_or_inf) {
        result.parts.exponent = 0x7FF;
        result.parts.fraction = 0;
        return result.raw;
    }

    // 3. Extract Mantissas and append the implicit leading 1 bit (Bit 52)
    // Assumes normalized inputs
    uint64_t m1 = (1ULL << 52) | f1.parts.fraction;
    uint64_t m2 = (1ULL << 52) | f2.parts.fraction;

    // 4. Calculate raw biased exponent (Subtract exponents and restore bias)
    int32_t exp_res = (int32_t)f1.parts.exponent - (int32_t)f2.parts.exponent + 1023;

    // 5. Divide the mantissas
    // Since m1 and m2 are roughly equal, m1 / m2 would yield 0 or 1.
    // We upscale m1 to 128 bits and shift it left by 52 positions first.
    // This allows integer division to compute the correct 53-bit fraction.
    unsigned __int128 dividend = (unsigned __int128)m1 << 53;
    unsigned __int128 quot = dividend / m2;
    unsigned __int128 remainder = dividend % m2;

    // 6. Normalize the quotient
    // In binary division, if m1 < m2, the quotient's implicit 1 drops to bit 51.
    // If m1 >= m2, the quotient's implicit 1 naturally sits at bit 52.
    uint64_t final_fraction = 0;

    // Check if the implicit bit sits at bit 53 (corresponds to m1 >= m2)
    if (quot & ((unsigned __int128)1 << 53)) {
        // Extract the 52-bit fraction
        final_fraction = (uint64_t)((quot >> 1) & 0xFFFFFFFFFFFFFLL);

        // Rounding bits
        bool round_bit = (quot & 1) != 0;
        bool sticky_bit = (remainder != 0);
        bool lsb = (final_fraction & 1) != 0;

        // IEEE 754 standard Round-to-Nearest, Ties-to-Even rule
        if (round_bit && (sticky_bit || lsb)) {
            final_fraction++;
            if (final_fraction > 0xFFFFFFFFFFFFFLL) { // Handle carry-out
                final_fraction = 0;
                exp_res += 1;
            }
        }
    } else {
        // Implicit bit sits at bit 52 (corresponds to m1 < m2)
        // No right shift needed for the fraction, but we need the sticky bit updated
        final_fraction = (uint64_t)(quot & 0xFFFFFFFFFFFFFLL);

        // We shifted left by 53 instead of 52, so bit 0 of quot is the actual round bit
        // However, since we didn't shift right, we must look at the remainder for the true sticky status
        // For the m1 < m2 case, we effectively need to look at what would happen if we didn't shift as far.
        // Let's re-align it perfectly:

        // To make it straightforward, let's normalize the 54-bit temporary quotient first:
        // If bit 53 is not set, we shift the entire quotient up by 1 bit to force the implicit bit to 53,
        // but we must adjust the remainder logic. Let's use a cleaner normalization pattern:

        // Shift left by 1 to align the implicit bit to bit 53
        quot <<= 1;
        exp_res -= 1;

        final_fraction = (uint64_t)((quot >> 1) & 0xFFFFFFFFFFFFFLL);
        bool round_bit = (quot & 1) != 0;
        bool sticky_bit = (remainder != 0);
        bool lsb = (final_fraction & 1) != 0;

        if (round_bit && (sticky_bit || lsb)) {
            final_fraction++;
            if (final_fraction > 0xFFFFFFFFFFFFFLL) {
                final_fraction = 0;
                exp_res += 1;
            }
        }
    }

    result.parts.fraction = final_fraction;

    // 7. Check for Overflow / Underflow boundaries
    if (exp_res >= 0x7FF) {
        // Overflow to Infinity
        result.parts.exponent = 0x7FF;
        result.parts.fraction = 0;
    } else if (exp_res <= 0) {
        // Underflow to Zero (Flushing subnormal results to zero)
        result.parts.exponent = 0;
        result.parts.fraction = 0;
    } else {
        // Valid normalized exponent
        result.parts.exponent = (uint64_t)exp_res;
    }

    return result.raw;
}