如何在不转为double/long double的情况下实现高精度汇编fmodf函数？

你的问题核心在于当|x/y| > 2^23时，单精度浮点数无法表示x/y的小数部分，导致trunc(x/y)等于x/y本身，最终计算出的余数完全错误。由于无法使用更高精度的浮点数，我们可以利用x86 FPU的扩展精度特性（即使输入是float，FPU内部会用80位计算），结合frem指令来实现，再将结果调整为符合fmodf定义的形式。

步骤1：处理特殊情况 #

首先按照IEEE-754标准处理边界输入：

• 如果y为0，返回NaN（单精度NaN的十六进制是0x7FC00000）
• 如果x是无穷大，返回NaN
• 如果x为0，返回0

步骤2：用FPU计算余数（frem指令） #

frem指令会计算x - y * n，其中n是最接近x/y的整数（四舍五入），FPU内部用80位扩展精度计算，完全不会有|x/y|>2^23时的精度丢失问题。

步骤3：将frem结果转换为fmodf结果 #

fmodf要求余数的符号与x相同，且绝对值小于|y|，而frem的余数绝对值小于等于|y|/2，符号可能与x相反（当x/y的小数部分≥0.5且为正，或≤-0.5且为负时）。转换逻辑如下：

1. 比较余数和x的符号：
   • 如果符号相同，余数就是fmodf的结果
   • 如果符号相反，给余数加上y（当x为正时）或减去y（当x为负时），本质等价于余数 += y * sign(x)

汇编示例（x86 32位） #

; 输入：x在eax（单精度float的内存地址或寄存器），y在ebx
; 输出：结果在eax（单精度float）

fmodf:
    ; 加载x和y到FPU栈
    fld     dword [eax]
    fld     dword [ebx]

    ; 处理y=0的情况
    ftst
    fstsw   ax
    sahf
    jz      return_nan

    ; 处理x是无穷大/NaN的情况
    fld     st1                  ; 复制x到ST0
    fldz
    fcompp
    fstsw   ax
    sahf
    jp      return_nan           ; x是NaN或无穷大
    ja      return_nan           ; x是正无穷
    jb      return_nan           ; x是负无穷

    ; 处理x=0的情况
    fld     st1                  ; 加载x到ST0
    ftst
    fstsw   ax
    sahf
    jz      return_zero

    ; 计算frem
    fxch    st1                  ; ST0=y，ST1=x
    frem                         ; ST0 = x - y*round(x/y)，ST1=y
    fxch    st1                  ; ST0=y，ST1=余数

    ; 比较余数和x的符号
    fld     dword [eax]          ; ST0=x，ST1=余数，ST2=y
    fmul    st1, st0             ; ST1=余数*x，判断符号
    ftst
    fstsw   ax
    sahf
    jns     sign_ok              ; 符号相同，无需调整

    ; 符号相反，调整余数：余数 += y * sign(x)
    fld     st2                  ; ST0=y，ST1=x，ST2=余数，ST3=y
    fabs
    fmul    st1, st0             ; ST1=x的绝对值
    fdiv    st1, st1             ; ST1=1.0（x正）或-1.0（x负）
    fmul    st0, st1             ; ST0=y * sign(x)
    faddp   st2, st0             ; ST2=余数 + y*sign(x)
    fxch    st2                  ; ST0=调整后的余数

sign_ok:
    ; 保存结果到eax指向的内存
    fstp    dword [eax]
    fstp    st0                  ; 清理栈
    fstp    st0
    ret

return_nan:
    mov     dword [eax], 0x7FC00000
    fstp    st0
    fstp    st0
    ret

return_zero:
    mov     dword [eax], 0x00000000
    fstp    st0
    fstp    st0
    ret

替代方案（无FPU扩展精度支持） #

如果你的汇编环境连FPU扩展精度都不支持，只能手动操作IEEE-754单精度的指数和尾数：

1. 提取x和y的符号位、指数、尾数（加上隐含的1，变成24位整数）
2. 计算指数差diff = ex - ey：
   • 若diff <=23，直接用原公式计算
   • 若diff >23，将x的尾数左移diff位（得到一个大整数），与y的尾数做整数取余，再将余数转换回单精度浮点数（符号位与x相同，指数设为ey + 余数的最高位偏移）

这种方法需要处理大整数运算（比如用64位寄存器存储左移后的尾数），但完全不依赖更高精度的浮点数类型。

内容的提问来源于stack exchange，提问作者Megumin