Prove trigonometric identities in nasm

Question

Given to prove trigonometric identity:

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

The following program written in NASM should print 1 when the identity is verified else 0. I always get 0 as the output. According to me, I think my logic is right. There exists some problem in the loading of memory or stack overflow.

[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

Testing one value doesn't prove an identity! If I evaluate x*x and x+x for x=2, does that mean I've proved that x^2 = x+x?

Testing with all possible floats (2^32 of them including NaNs) is a good way to check that a function you've written works properly for all possible corner cases, but floating point numbers still aren't the same thing as mathematical real numbers. Esp. when the results have significant errors (see below about fsin).

Trying something like this and finding it works for all floats is confirmation that it's worth looking for mathematical proof, not that you've found it.

---

Also note that on Intel CPUs (and apparently all other x86 CPUs except AMD k5) the fsin insn is not accurate for inputs near Pi. In the worst case, the error is 1.37 quintillion units in the last place, leaving fewer than four bits correct, according to Bruce Dawson's excellent blog post. This is due to range-reduction with only a 66 bit Pi constant, and can't be fixed for backwards-compatibility reasons. (i.e. the exact behaviour of fsin is essentially part of the ISA, warts and all).

The standard add/sub/mul/div/sqrt operations all produce results with at most 0.5 ulp of error (i.e. even the last bit of the mantissa is rounded correctly).

Bruce has a whole series of FP articles, with an index in this one about FP comparisons. I've added links to some of those articles to the x86 tag wiki, because they have such good info. They're not all about x87; most of them mention how things are different between x87 and SSE when it matters for current versions of MSVC and/or gcc.

---

Fun fact: A floating point trig identity:

sin(double(pi)) != 0.0. You shouldn't expect it to be, because double can't exactly represent the value of the irrational number pi. However,

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

(If sin(double(pi)) is evaluated accurately, not with x87 fsin, and you do the addition with more precision than double)

Quoting Bruce Dawson again:

This works because the sine of a number very close to pi is almost equal to the error in the estimate of pi. This is just calculus 101, and a variant of Newton’s method, but I still find it charming.

For more, see Bruce's FP comparison article, and search in it for sin(pi). He has a whole section about this floating-point trig identity.

Answer 2

EDIT: Found!
In answer to your question, your code has one instruction wrong: I've added comments...

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

The second-last line should be fsub, not fmul

---

General comments:
You've got a symbol pi in your code, and it is only set to 3.14 - luckily you don't use it in your code. Still, don't do this: the x87 knows all about pi: FLDPI

Another thing you're doing correctly is comparing against a threshold - still, 1e-1 is a pretty large delta...

You have a much more fundamental problem though, and you mention it when you talk about "stack overflow". Your code is storing values in the (program) stack without making room for them first:

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

This last one could crash your program: esp+4 could be past the top of the stack! If you want to use mov rather than push, you should start your program with sub esp,8 to make room for two dwords first.

And you do know that the x87 already has an internal stack of 8 floating point registers? You can 'save away' values inside that stack while working out other interim results, then get those values back again. Storing them to memory (such as the first line of code above) actually reduces the precision of your interim results.