; quad.asm: integer power function callable from Borland C++.

.MODEL LARGE,C

.CODE

; int solve_quadratic(double a, double b, double c, double *x1, double *x2);
; and finds the roots of that polynomial. If there are two real
; roots, it writes them back to x1 and x2 and returns 1. If the
; discriminant b^2-4*a*c is less than 0, the roots are not real
; The quadratic formula is as follows:
; (-b +/- sqrt(b*b-4*a*c)) / 2*a

; separated by | signs
; Stack top is at left
fld A ; a Here, ST(0) contains a.
fld B ; b | a Now ST(1) has a. Etc.

; Negate b -- squaring it is sign-independent, and we need b
; negated in all the other instances in the formula.
fchs

fld C ; c | b | a
fld st(1) ; b | c | b | a
fmul st,st(0) ; b*b | c | b | a

; Just plain fxch has an implicit operand of ST(1).

fxch ; c | b*b | b | a
fadd st,st(0) ; 2*c | b*b | b | a
fadd st,st(0) ; 4*c | b*b | b | a
fmul st,st(3) ; 4*a*c | b*b | b | a
fsubp st(1), st ; b*b-4*a*c | b | a
ftst ; Compare against 0

; We've computed b*b-4*a*c. If negative, we return 0.
; To do the comparison, we have to store the 80x87 status
; word and use sahf to store it into the flags. Once it's
; in the flags, we can jump on above or equal (jae) to jump if the
; number tested is greater than 0 after the FTST instruction above.
sub sp,2 ; Quick, allocate a local
mov bx,sp ; Point BX at it
fstsw ss:[bx] ; Store the 80x87's status word there
mov ax,ss:[bx] ; AX <- status word
add sp,2 ; Deallocate the local
sahf ; Get SW into flags
jae ComputeResults ; Jump if positive
fstp st ; This instruction clears the stack
fstp st ; we've got three values on the stack
fstp st ; to clear
xor ax,ax ; Return 0 - no roots found.
ComputeResults:
fsqrt ; Find square root of discriminant
fxch st(2) ; a | b | sqrt(b*b-4*a*c)
fadd st,st(0) ; 2*a | b | sqrt(b*b-4*a*c)
fld st(1) ; b | 2*a | b | sqrt(b*b-4*a*c)
fadd st,st(3) ; b+sqrt(b*b-4*a*c) | 2*a | b | sq..
fdiv st,st(1) ; x1 | 2*a | b | sqrt(b*b-4*a*c)
les bx,X1 ;
fstp qword ptr es:[bx] ; 2*a | b | sqrt(b*b-4*a*c)
fxch ; b | 2*a | sqrt(b*b-4*a*c)
fsubrp st(2),st ; 2*a | -b - sqrt(b*b-4*a*c)
fdiv ; x2
les bx,X2 ; Store x2
fstp qword ptr es:[bx]
mov ax,1 ; Return 1