在 nasm 中证明三角恒等式

Question

用于证明三角恒等式：

cos(fi+theta) = cos(theta)cos(fi)-sin(theta).sin(fi)  

以下用 NASM 编写的程序应该在验证身份时打印 1，否则打印 0。我总是得到 0 作为输出。在我看来，我认为我的逻辑是正确的。内存加载或堆栈溢出存在问题。

[bits 32]

extern printf
extern _exit


section .data

hello: db "Value =  %d ", 10, 0
theta: dd 1.0
fi: dd 2.0
threshold: dd 1e-1
SECTION .text                   

global  main                ; "C" main program 

main:
start:


; compute cos(theta)cos(fi)-sin(theta).sin(fi)
fld dword [theta]
fcos
fld dword [fi]
fcos
fmul st0,st1
fstp dword [esp]
fld dword [theta]
fsin
fld dword [fi]
fsin
fmul st0,st1
fld dword [esp]
fmul st0,st1
fstp dword [esp] 

;compute cos(fi+theta)
fld dword [theta]
fld dword [fi]
fadd st0,st1
fcos

; compare
fld dword [esp]  
fsub st0, st1
fld dword [threshold]
fcomi st0, st1
ja .equal
mov eax, 0
jmp .exit

.equal:
    mov eax, 1
.exit:

    mov     dword[ esp ],       hello
    mov     dword [ esp + 4 ],   eax
    call   printf

    ; Call 'exit': exit( 0 );
    mov     dword[ esp ],       0
    call   _exit

Answer 1

测试一个值并不能证明一个身份！如果我为x=2 评估x*x 和x+x，这是否意味着我已经证明了x^2 = x+x？

测试所有可能的floats（其中 2^32 个包括 NaN）是检查您编写的函数是否适用于所有可能的极端情况的好方法，但浮点数仍然不是与数学实数相同。特别是。当结果有重大错误时（见下文关于fsin）。

尝试这样的事情并发现它适用于所有floats 是确认值得寻找数学证明，而不是你已经找到它。

---

另请注意，在 Intel CPU（显然除了AMD k5 之外的所有其他 x86 CPU）上，fsin insn 不对于 Pi 附近的输入是准确的。根据 Bruce Dawson 的优秀博客文章，在最坏的情况下，错误是 1.37 quintillion units in the last place, leaving fewer than four bits correct。这是由于只有 66 位 Pi 常数的范围缩小，并且由于向后兼容的原因无法修复。 （即fsin 的确切行为本质上是 ISA 的一部分，是缺点和全部）。

标准的 add/sub/mul/div/sqrt 操作都产生最多 0.5 ulp 误差的结果（即，即使尾数的最后一位被正确舍入）。

Bruce 拥有完整系列的 FP 文章，索引在 this one about FP comparisons。我已将其中一些文章的链接添加到 x86 标签 wiki，因为它们有很好的信息。它们不全是关于 x87 的。他们中的大多数人提到了 x87 和 SSE 之间的不同之处，这对于当前版本的 MSVC 和/或 gcc 很重要。

---

有趣的事实：浮点三角恒等式：

sin(double(pi)) != 0.0。你不应该期望它是，因为double 不能准确地表示无理数 pi 的值。然而，

 pi ~= double(pi) + sin(double(pi)) 

（如果 sin(double(pi)) 被准确评估，而不是 x87 fsin，并且您执行的加法比 double 更精确）

再次引用布鲁斯·道森的话：

这是有效的，因为一个非常接近 pi 的数字的正弦几乎等于 pi 的估计误差。这只是微积分 101，是牛顿方法的一种变体，但我仍然觉得它很迷人。

更多信息，see Bruce's FP comparison article，并在其中搜索sin(pi)。他有一整节关于这个浮点三角恒等式。

Answer 2

编辑：找到！
在回答您的问题时，您的代码有一条指令错误：我添加了 cmets...

fld dword [theta]  ; Load theta
fcos               ; Get its cos()
fld dword [fi]     ; Load fi
fcos               ; Get its cos()
fmul st0,st1       ; Multiply them
fstp dword [esp]   ; Store away on stack [Why?]

fld dword [theta]  ; Load theta
fsin               ; Get its sin()
fld dword [fi]     ; Load fi
fsin               ; Get its sin()
fmul st0,st1       ; Multiply them

fld dword [esp]    ; Get what was stored away
fmul st0,st1       ; Multiply them??? [ You meant fsub!]
fstp dword [esp]   ; Store it away

倒数第二行应该是fsub，而不是fmul

---

一般 cmets：
你的代码中有一个符号pi，它只设置为3.14——幸运的是你没有在你的代码中使用它。不过，不要这样做：x87 对 pi 了如指掌：FLDPI

您正确做的另一件事是与阈值进行比较 - 不过，1e-1 是一个相当大的增量...

你有一个更多更根本的问题，当你谈论“堆栈溢出”时你提到了它。您的代码将值存储在（程序）堆栈中，而没有先为它们腾出空间：

fstp dword [esp]
...
fld dword [esp]
...
mov     dword[ esp ],       hello
mov     dword [ esp + 4 ],   eax
call   printf

最后一个可能会使您的程序崩溃：esp+4 可能超出堆栈顶部！如果你想使用mov 而不是push，你应该以sub esp,8 开始你的程序，以便首先为两个dwords 腾出空间。

您知道 x87 已经有一个由 8 个浮点寄存器组成的内部堆栈吗？您可以在计算其他中间结果时“保存”该堆栈中的值，然后再次取回这些值。将它们存储到内存中（例如上面的第一行代码）实际上会降低中间结果的精度。