From 84ebd52c2493a38db67d8167fc307ba810cf0b5d Mon Sep 17 00:00:00 2001 From: Chen Peng Date: Sat, 29 Aug 2026 11:40:02 +0800 Subject: [PATCH] =?UTF-8?q?=E6=B5=8B=E8=AF=95=E7=94=A8=E4=BE=8B?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Foundation/c_Float.c | 915 ++++++++++++++++++++++++ Foundation/c_Float.h | 211 ++++++ Foundation/{c_float.t.c => c_Float.t.c} | 0 3 files changed, 1126 insertions(+) create mode 100644 Foundation/c_Float.c create mode 100644 Foundation/c_Float.h rename Foundation/{c_float.t.c => c_Float.t.c} (100%) diff --git a/Foundation/c_Float.c b/Foundation/c_Float.c new file mode 100644 index 0000000..281b872 --- /dev/null +++ b/Foundation/c_Float.c @@ -0,0 +1,915 @@ +#include +#include + +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; +} + + + diff --git a/Foundation/c_Float.h b/Foundation/c_Float.h new file mode 100644 index 0000000..93b6ff0 --- /dev/null +++ b/Foundation/c_Float.h @@ -0,0 +1,211 @@ +#ifndef INCLUDED_C_FLOAT_H +#define INCLUDED_C_FLOAT_H + +#ifndef INCLUDED_C_TYPES_H +#include +#endif /*INCLUDED_C_TYPES_H*/ + +#ifndef INCLUDED_MATH_H +#define INCLUDED_MATH_H +#include +#endif /*INCLUDED_MATH_H*/ + +#ifndef INCLUDED_FLOAT_H +#define INCLUDED_FLOAT_H +#include +#endif /*INCLUDED_FLOAT_H*/ + + +/* ------------------------------------------------------------------------------------------------------------------ */ +/* */ + +typedef union { + float f; + uint32_t raw; + struct { + uint32_t fraction : 23; // 尾数 (M) + uint32_t exponent : 8; // 指数 (E) + uint32_t sign : 1; // 符号位 (S) + } parts; +} c_Float_t; + +typedef union { + double d; + uint64_t raw; + struct { + uint64_t fraction : 52; // 尾数 + uint64_t exponent : 11; // 指数 + uint64_t sign : 1; // 符号 + } parts; +} c_Double_t; + +/* ------------------------------------------------------------------------------------------------------------------ */ +/* */ + +// 32位单精度常量定义 +#define C_FLOAT_SIGN_MASK 0x80000000U +#define C_FLOAT_EXP_MASK 0x7F800000U +#define C_FLOAT_FRAC_MASK 0x007FFFFFU +#define C_FLOAT_HIDDEN_BIT 0x00800000U // 隐藏的最高位1 +#define C_FLOAT_EXP_BIAS 127 + +#define C_FLOAT_NEG_INF 0xFF800000U +#define C_FLOAT_POS_INF 0x7F800000U + +// 快捷提取宏 +#define C_FLOAT_GET_SIGN(u) (((u) & C_FLOAT_SIGN_MASK) >> 31) +#define C_FLOAT_GET_EXP(u) (((u) & C_FLOAT_EXP_MASK) >> 23) +#define C_FLOAT_GET_FRAC(u) ((u) & C_FLOAT_FRAC_MASK) + +#define C_DOUBLE_SIGN_MASK 0x8000000000000000ULL +#define C_DOUBLE_EXP_MASK 0x7FF0000000000000ULL +#define C_DOUBLE_FRAC_MASK 0x000FFFFFFFFFFFFFULL +#define C_DOUBLE_HIDDEN_BIT 0x0010000000000000ULL // 第52位（从0开始算） + +#define C_DOUBLE_POS_INF 0x7FF0000000000000ULL +#define C_DOUBLE_NEG_INF 0xFFF0000000000000ULL + +#define C_DOUBLE_GET_SIGN(u) (((u) & C_DOUBLE_SIGN_MASK) >> 63) +#define C_DOUBLE_GET_EXP(u) (((u) & C_DOUBLE_EXP_MASK) >> 52) +#define C_DOUBLE_GET_FRAC(u) ((u) & C_DOUBLE_FRAC_MASK) + +/* ------------------------------------------------------------------------------------------------------------------ */ +/* */ + +C_STATIC_FORCE_INLINE +uint32_t c_Float_Pack(const uint32_t sign, const uint32_t exp, const uint32_t frac) { + return ((sign << 31) & C_FLOAT_SIGN_MASK) | + ((exp << 23) & C_FLOAT_EXP_MASK) | + (frac & C_FLOAT_FRAC_MASK); +} + +C_STATIC_FORCE_INLINE +int c_Float_IsNAN(uint32_t raw) { + return ((raw & C_FLOAT_EXP_MASK) == C_FLOAT_EXP_MASK) && ((raw & C_FLOAT_FRAC_MASK) != 0); +} + +C_STATIC_FORCE_INLINE +int c_Double_IsNAN(uint64_t raw) { + return ((raw & C_DOUBLE_EXP_MASK) == C_DOUBLE_EXP_MASK) && ((raw & C_DOUBLE_FRAC_MASK) != 0); +} + +/* ------------------------------------------------------------------------------------------------------------------ */ +/* */ + +uint32_t c_Float_Add(uint32_t a, uint32_t b); + +uint32_t c_Float_Mul(uint32_t a, uint32_t b); + +uint32_t c_Float_Div(uint32_t a, uint32_t b); + +int c_Float_Cmp(uint32_t a, uint32_t b); + +C_STATIC_FORCE_INLINE +uint32_t c_Float_Sub(uint32_t a, uint32_t b) { + // 透過與 0x80000000 進行 XOR，直接將 b 的符號位元取反 (0->1, 1->0) + // 隨後將 A - B 轉換為 A + (-B) 傳入加法器 + return c_Float_Add(a, b ^ C_FLOAT_SIGN_MASK); +} + +C_STATIC_FORCE_INLINE +bool c_Float_IsZero(uint32_t raw) { + // Strip away the sign bit; check if the remaining 31 bits are 0 + return (raw & ~C_FLOAT_SIGN_MASK) == 0U; +} + +C_STATIC_FORCE_INLINE +bool c_Float_IsInf(uint32_t raw) { + // Strip the sign bit and check if it exactly matches the exponent mask. + // If any fraction bits were set, it would be a NaN instead of Infinity. + return (raw & ~0x80000000U) == C_FLOAT_EXP_MASK; +} + +/** + * @brief Determines if the float is specifically Negative Infinity (-Inf). + */ +C_STATIC_FORCE_INLINE +bool c_Float_IsNegInf(uint32_t raw) { + return raw == C_FLOAT_NEG_INF; +} + +/** + * @brief Determines if the float is specifically Positive Infinity (+Inf). + */ +C_STATIC_FORCE_INLINE +bool c_Float_IsPosInf(uint32_t raw) { + return raw == C_FLOAT_POS_INF; +} + +/* ------------------------------------------------------------------------------------------------------------------ */ +/* */ + +uint64_t c_Double_Pack(uint32_t sign, int32_t exp, uint64_t frac64); + +int c_Double_Cmp(uint64_t a, uint64_t b); + +uint64_t c_Double_Add(uint64_t a, uint64_t b); + +C_STATIC_FORCE_INLINE +uint64_t c_Double_Sub(uint64_t a, uint64_t b) { + // A - B == A + (-B) + return c_Double_Add(a, b ^ C_DOUBLE_SIGN_MASK); +} + +uint64_t c_Double_Mul(uint64_t a, uint64_t b); + +uint64_t c_Double_Div(uint64_t a, uint64_t b); + + +C_STATIC_FORCE_INLINE +bool c_Double_IsZero(const uint64_t raw) { + // Strip away the sign bit; check if the remaining 63 bits are 0 + return (raw & ~C_DOUBLE_SIGN_MASK) == 0ULL; +} + +C_STATIC_FORCE_INLINE +bool c_Double_IsInf(uint64_t raw) { + // Strip the sign bit and check if it exactly matches the exponent mask. + // If any fraction bits were set, it would be a NaN instead of Infinity. + return (raw & ~C_DOUBLE_SIGN_MASK) == C_DOUBLE_EXP_MASK; +} + +/** + * @brief Determines if the double is specifically Negative Infinity (-Inf). + */ +C_STATIC_FORCE_INLINE +bool c_Double_IsNegInf(uint64_t raw) { + return raw == C_DOUBLE_NEG_INF; +} + +/** + * @brief Determines if the double is specifically Positive Infinity (+Inf). + */ +C_STATIC_FORCE_INLINE +bool c_Double_IsPosInf(uint64_t raw) { + return raw == C_DOUBLE_POS_INF; +} + +/* ------------------------------------------------------------------------------------------------------------------ */ +/* */ + +C_STATIC_FORCE_INLINE +int c_float_cmp(const float a, const float b) { + c_Float_t va; + c_Float_t vb; + va.f = a; + vb.f = b; + return c_Float_Cmp(va.raw, vb.raw); +} + +C_STATIC_FORCE_INLINE +int c_double_cmp(const double a, const double b) { + c_Double_t va; + c_Double_t vb; + va.d = a; + vb.d = b; + return c_Double_Cmp(va.raw, vb.raw); +} + + + +#endif /*INCLUDED_C_FLOAT_H*/ diff --git a/Foundation/c_float.t.c b/Foundation/c_Float.t.c similarity index 100% rename from Foundation/c_float.t.c rename to Foundation/c_Float.t.c