#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// function for checking that assembly code is computing the correct result
double quadraticRootC(double a, double b, double c)
{
    return (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
}

double quadraticRoot(double a, double b, double c)
{
    double root;

    asm(
        "fldl       %1       # a                   \n"
        "fadd       %%ST     # 2a                  \n"
        "fldl       %1       # a 2a                \n"
        "fldl       %3       # c a 2a              \n"
        "fmulp      %%ST(1)  # ac 2a               \n"
        "fadd       %%ST     # 2ac 2a              \n"
        "fadd       %%ST     # 4ac 2a              \n"
        "fldl       %2       # b 4ac 2a            \n"
        "fldl       %2       # b b 4ac 2a          \n"
        "fmulp      %%ST(1)  # b^2 4ac 2a          \n"
        "fsubp      %%ST(1)  # b^2-4ac 2a          \n"
        "fsqrt               # sqrt(b^2-4ac) 2a    \n"
        "fldl       %2       # b sqrt(b^2-4ac) 2a  \n"
        "fchs                # -b sqrt(b^2-4ac) 2a \n"
        "faddp      %%ST(1)  # -b+sqrt(b^2-4ac) 2a \n"
        "fdivp      %%ST(1)  # -b+sqrt(b^2-4ac)/2a \n"
        :"=t"(root)
        :"m"(a), "m"(b), "m"(c)
        );
    return(root);
}

int main(int argc, char **argv)
{
    double  a, b, c;
    double  root, rootC;

    a = 3.0;
    b = -5.0;
    c = 1.0;
    root = quadraticRoot(a, b, c);
    rootC = quadraticRootC(a, b, c);

    printf("quadraticRoot(%.3f, %.3f, %.3f) = %.3f, %.3f\n", a, b, c, root, rootC);

    return 0;
}