; quad.asm: integer power function callable from Borland C++.
; Copyright (C) 1991 by Nicholas Wilt.  All rights reserved.

.MODEL  LARGE,C

.CODE

; int solve_quadratic(double a, double b, double c, double *x1, double *x2);
; solve_quadratic takes the coefficients of a quadratic polynomial
; 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
; and solve_quadratic returns 0.
; The quadratic formula is as follows:
;   (-b +/- sqrt(b*b-4*a*c)) / 2*a

    PUBLIC  solve_quadratic

solve_quadratic PROC    A:QWORD,B:QWORD,C:QWORD,X1:DWORD,X2:DWORD
                ; Comments show stack contents
                ; 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.
    jmp short LeaveQuadratic
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
LeaveQuadratic:
    ret         ; Return
solve_quadratic ENDP

    END