Programming

1 ; function quadratic
2 ; finds solutions to the quadratic equation:
3 ; a*x^2 + b*x + c = 0
4 ; C prototype:
5 ; int quadratic( double a, double b, double c,
6 ; double * root1, double *root2 )
7 ; Parameters:
8 ; a, b, c - coefficients of powers of quadratic equation (see above)
9 ; root1 - pointer to double to store first root in
10 ; root2 - pointer to double to store second root in
11 ; Return value:
12 ; returns 1 if real roots found, else 0
13
14 %define a qword [ebp+8]
15 %define b qword [ebp+16]
16 %define c qword [ebp+24]
17 %define root1 dword [ebp+32]
18 %define root2 dword [ebp+36]
19 %define disc qword [ebp-8]
132 CHAPTER 6. FLOATING POINT
20 %define one_over_2a qword [ebp-16]
21
22 segment .data
23 MinusFour dw -4
24
25 segment .text
26 global _quadratic
27 _quadratic:
28 push ebp
29 mov ebp, esp
30 sub esp, 16 ; allocate 2 doubles (disc & one_over_2a)
31 push ebx ; must save original ebx
32
33 fild word [MinusFour]; stack -4
34 fld a ; stack: a, -4
35 fld c ; stack: c, a, -4
36 fmulp st1 ; stack: a*c, -4
37 fmulp st1 ; stack: -4*a*c
38 fld b
39 fld b ; stack: b, b, -4*a*c
40 fmulp st1 ; stack: b*b, -4*a*c
41 faddp st1 ; stack: b*b - 4*a*c
42 ftst ; test with 0
43 fstsw ax
44 sahf
45 jb no_real_solutions ; if disc < 0, no real solutions
46 fsqrt ; stack: sqrt(b*b - 4*a*c)
47 fstp disc ; store and pop stack
48 fld1 ; stack: 1.0
49 fld a ; stack: a, 1.0
50 fscale ; stack: a * 2^(1.0) = 2*a, 1
51 fdivp st1 ; stack: 1/(2*a)
52 fst one_over_2a ; stack: 1/(2*a)
53 fld b ; stack: b, 1/(2*a)
54 fld disc ; stack: disc, b, 1/(2*a)
55 fsubrp st1 ; stack: disc - b, 1/(2*a)
56 fmulp st1 ; stack: (-b + disc)/(2*a)
57 mov ebx, root1
58 fstp qword [ebx] ; store in *root1
59 fld b ; stack: b
60 fld disc ; stack: disc, b
61 fchs ; stack: -disc, b
6.3. THE NUMERIC COPROCESSOR 133
62 fsubrp st1 ; stack: -disc - b
63 fmul one_over_2a ; stack: (-b - disc)/(2*a)
64 mov ebx, root2
65 fstp qword [ebx] ; store in *root2
66 mov eax, 1 ; return value is 1
67 jmp short quit
68
69 no_real_solutions:
70 mov eax, 0 ; return value is 0
71
72 quit:
73 pop ebx
74 mov esp, ebp
75 pop ebp
76 ret

gcc -O3 -S code.c