x87 可以对 UNsigned QUADword 整数执行精确除法吗？

Question

8086/8087/8088 宏汇编语言参考手册 (c) 1980 Intel Corporation 提到（重点是我的）：

... 8087 提供了非常好的实数系统近似值。然而，重要的是要记住，它不是精确的表示，并且实数的算术本质上是近似的。
相反，同样重要的是，8087 确实对其实数的整数子集执行精确算术。也就是说，对两个整数进行运算会返回精确的积分结果，前提是真实结果是整数并且在 range 内。

最近的手册更加简洁（强调他们的）：

IA 处理器...它们可以处理最多 18 位的十进制数，而不会出现舍入错误，对大至 2^64（或 10^18）的整数执行精确算术。

FPU 支持的整数数据类型包括有符号字（16 位）、有符号双字（32 位）和有符号 qword（64 位）。从来没有提到过 UNsigned。事实上，FPU 的一切都带有符号性，甚至支持带符号零（+0 和 -0）。
那么，是否可以使用 FPU 将几个无符号64 位数字相除并得到精确的商和余数？

对于几个有符号64 位数字的除法，我编写了下面的代码。商看起来不错，但余数总是返回零。为什么是这样？

; IN (edx:eax,ecx:ebx) OUT (edx:eax,ecx:ebx,CF)
FiDiv:  push    edi ecx ebx edx eax
        mov     edi, esp
        fninit
        fild    qword [edi]     ; Dividend
        fild    qword [edi+8]   ; Divisor
        fld
        fnstcw  [edi]
        or      word [edi], 0C00h ; Truncate Towards Zero
        fldcw   [edi]
        fdivr   st2             ; st0 = st2 / st0
        fld                     ; Duplicate because `fistp` does pop
        fistp   qword [edi]     ; Quotient
        fmulp                   ; st1 *= st0, pop st0
        fsubp                   ; st1 -= st0, pop st0
        fistp   qword [edi+8]   ; Remainder
        fnstsw  ax
        test    ax, 101b        ; #Z (Divide-By-Zero), #I (Invalid Arithmetic)
        pop     eax edx ebx ecx edi
        jnz     .NOK
        ret                     ; CF=0
.NOK:   stc                     ; Overflow on 8000000000000000h / -1
        ret                     ; or Division by zero
; ------------------------------

FPU 舍入模式设置为“向零截断”又名“截断”，以模仿 ALUidiv指令的行为。

Answer 1

余数为零

fdivr   st2             ; st0 = st2 / st0
fld                     ; Duplicate because `fistp` does pop
fistp   qword [edi]     ; Quotient
fmulp                   ; st1 *= st0, pop st0
fsubp                   ; st1 -= st0, pop st0
fistp   qword [edi+8]   ; Remainder

此代码计算余数：

Remainder = Dividend - (Quotient * Divisor)

由于舍入模式设置为“向零截断”，因此该fistp qword [edi]指令将在存储到内存之前将 ST0 中保存的商（商的副本）转换为整数。然而，保留在 (fpu) 堆栈上的商值仍然是带分数的实数。一旦与除数相乘，它将再次产生被除数，导致余数为零。
缺少的是将商舍入为整数并已在 (fpu) 堆栈上执行此操作：

fdivr   st2             ; st0 = st2 / st0
frndint
fld                     ; Duplicate because `fistp` does pop
fistp   qword [edi]     ; Quotient
fmulp                   ; st1 *= st0, pop st0
fsubp                   ; st1 -= st0, pop st0
fistp   qword [edi+8]   ; Remainder

但更快的方法是简单地从内存中重新加载整数商：

fdivr   st2             ; st0 = st2 / st0
fistp   qword [edi]     ; Quotient
fild    qword [edi]
fmulp                   ; st1 *= st0, pop st0
fsubp                   ; st1 -= st0, pop st0
fistp   qword [edi+8]   ; Remainder

无符号除法

在内部，FPU 将 64 位专用于数字的有效数，加上一个单独的位用于数字的符号。FPU 可以表示从 -18'446744073'709551616 到 18'446744073'709551616 范围内的每个整数。64 位有效数允许我们处理范围从 0 到 18'446744073'709551615 的无符号 64 位整数。唯一的麻烦是如何加载和存储这些fild和fistp无法处理的值（因为它们被限制在 -9'223372036'854775808 到 9'223372036'854775807 的范围内操作）。
可以在无符号四字和扩展实数格式之间来回转换，因此我们可以使用fldandfstp来代替。另一种方法是从上半部分和下半部分加载/存储无符号四字。但转换需要时间，所以我发现通过消除足够多的特殊情况，剩下的唯一麻烦的操作就是股息的加载。其他一切都可以照常fild使用。fistp

特殊情况包括：

• 除以零：CF=1
• 除以一：商=被除数，余数=0
• 被除数等于除数的除法：商=1 余数=0
• 被除数小于除数的除法：商=0，余数=被除数

在需要实际值的情况下fdiv，代码首先加载一半的被除数，将其加倍返回 (fpu) 堆栈，如果真正的被除数是奇数，则有条件地加 1。

; IN (edx:eax,ecx:ebx) OUT (edx:eax,ecx:ebx,CF)
FuDiv:  cmp     ebx, 1
        jbe     .TST            ; Divisor could be 0 or 1
.a:     cmp     edx, ecx
        jb      .LT             ; Dividend < Divisor
        ja      .b              ; Dividend > Divisor
        cmp     eax, ebx
        jb      .LT             ; Dividend < Divisor
        je      .GE             ; Dividend = Divisor
.b:     test    ecx, ecx
        js      .GE             ; Dividend > Divisor > 7FFFFFFFFFFFFFFFh

        shr     edx, 1          ; Halving the unsigned 64-bit Dividend
        rcr     eax, 1          ; (CF)      to get in range for `fild`
        push    edi ecx ebx edx eax
        mov     edi, esp
        fninit
        fild    qword [edi]     ; st0 = int(Dividend / 2)
        fadd    st0             ; st0 = {Dividend - 1, Dividend}
        jnc     .c              ; (CF)
        fld1
        faddp                   ; st0 = Dividend [0, FFFFFFFFFFFFFFFFh]
.c:     fild    qword [edi+8]   ; Divisor is [2, 7FFFFFFFFFFFFFFFh]
        fld
        fnstcw  [edi]
        or      word [edi], 0C00h ; Truncate Towards Zero
        fldcw   [edi]
        fdivr   st2             ; st0 = st2 / st0
        fistp   qword [edi]     ; Quotient
        fild    qword [edi]
        fmulp                   ; st1 *= st0, pop st0
        fsubp                   ; st1 -= st0, pop st0
        fistp   qword [edi+8]   ; Remainder
        pop     eax edx ebx ecx edi
        ret                     ; CF=0

.TST:   test    ecx, ecx
        jnz     .a
        cmp     ebx, 1          ; Divisor is 0 or 1
        jb      .NOK
.ONE:   dec     ebx             ; Remainder ECX:EBX is 0
.NOK:   ret

.GE:    sub     eax, ebx        ; Remainder ECX:EBX is Dividend - Divisor
        sbb     edx, ecx
        mov     ebx, eax
        mov     ecx, edx
        mov     eax, 1          ; Quotient EDX:EAX is 1
        cdq
        ret                     ; CF=0

.LT:    mov     ebx, eax        ; Remainder ECX:EBX is Dividend
        mov     ecx, edx
        xor     eax, eax        ; Quotient EDX:EAX is 0
        cdq
        ret                     ; CF=0
; ------------------------------

更快获得结果

引用一个答案，该答案使用名为 Chunking 的技术来划分几个 64 位整数：

即使您使用 64 位数据类型，实践表明（通用）程序中的大多数除法仍然可以仅使用内置指令来完成div。这就是为什么我在代码中添加了一个检测机制作为前缀，该机制检查 ECX:EBX 中的除数是否小于 4GB（因此适合 EBX），以及 EDX 中的除数扩展是否小于 EBX 中的除数。如果满足这些条件，正常的div指令就会完成工作，而且速度也会更快。如果由于某种原因（例如学校）div不允许使用，则只需删除前缀代码即可。

今天的代码可以受益于相同的前缀，但事实证明，首先继续检测特殊情况会更有利：

; IN (edx:eax,ecx:ebx) OUT (edx:eax,ecx:ebx,CF)
FuDiv:  cmp     ebx, 1
        jbe     .TST            ; Divisor could be 0 or 1
.a:     cmp     edx, ecx
        jb      .LT             ; Dividend < Divisor
        ja      .b              ; Dividend > Divisor
        cmp     eax, ebx
        jb      .LT             ; Dividend < Divisor
        je      .GE             ; Dividend = Divisor
.b:     test    ecx, ecx
        js      .GE             ; Dividend > Divisor > 7FFFFFFFFFFFFFFFh
; - - - - - - - - - - - - - - -
        jnz     .fdiv
        cmp     edx, ebx
        jnb     .fdiv
.div:   div     ebx             ; EDX:EAX / EBX --> EAX Quotient, EDX Remainder
        mov     ebx, edx        ; Remainder to ECX:EBX
        xor     edx, edx        ; Quotient to EDX:EAX
        ret                     ; CF=0
; - - - - - - - - - - - - - - -
.fdiv:  shr     edx, 1          ; Halving the unsigned 64-bit Dividend
        rcr     eax, 1          ; (CF)      to get in range for `fild`
        push    edi ecx ebx edx eax
        mov     edi, esp
        fninit
        fild    qword [edi]     ; st0 = int(Dividend / 2)
        fadd    st0             ; st0 = {Dividend - 1, Dividend}
        jnc     .c              ; (CF)
        fld1
        faddp                   ; st0 = Dividend [0, FFFFFFFFFFFFFFFFh]
.c:     fild    qword [edi+8]   ; Divisor is [2, 7FFFFFFFFFFFFFFFh]
        fld
        fnstcw  [edi]
        or      word [edi], 0C00h ; Truncate Towards Zero
        fldcw   [edi]
        fdivr   st2             ; st0 = st2 / st0
        fistp   qword [edi]     ; Quotient
        fild    qword [edi]
        fmulp                   ; st1 *= st0, pop st0
        fsubp                   ; st1 -= st0, pop st0
        fistp   qword [edi+8]   ; Remainder
        pop     eax edx ebx ecx edi
        ret                     ; CF=0

.TST:   test    ecx, ecx
        jnz     .a
        cmp     ebx, 1          ; Divisor is 0 or 1
        jb      .NOK
.ONE:   dec     ebx             ; Remainder ECX:EBX is 0
.NOK:   ret

.GE:    sub     eax, ebx        ; Remainder ECX:EBX is Dividend - Divisor
        sbb     edx, ecx
        mov     ebx, eax
        mov     ecx, edx
        mov     eax, 1          ; Quotient EDX:EAX is 1
        cdq
        ret                     ; CF=0

.LT:    mov     ebx, eax        ; Remainder ECX:EBX is Dividend
        mov     ecx, edx
        xor     eax, eax        ; Quotient EDX:EAX is 0
        cdq
        ret                     ; CF=0
; ------------------------------