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

main()
{
        int n,j,ier;
        double a,c,b,x1,x2,discr,err1,err2;

        cprini("stdout","fort.13");
/* 
*       create the coefficients of the equation 
*/

        a=-2.1;
        a=2.1;

        b=1;
        c=-1;

//        a=1.0e-12;

        cprin2("a as created is",&a,1);
        cprin2("b as created is",&b,1);
        cprin2("c as created is",&c,1);

        x1=7;
        x2=8;

        quad_solve(&ier,a,b,c,&x1,&x2,&err1,&err2);

        cprinf("after quad_solve, ier=",&ier,1);




        cprin2("after quad_solve, x1=",&x1,1);
        cprin2("after quad_solve, x2=",&x2,1);

        cprin2("after quad_solve err1=",&err1,1);
        cprin2("after quad_solve err2=",&err2,1);

}


int quad_solve(ierp,a,b,c,x1p,x2p,err1p,err2p)
double a,b,c,*x1p,*x2p,*err1p,*err2p;
int *ierp;
{
double x1,x2,discr,err1,err2;
/* 
*       This function solves the quadratic equation with 
*             the coefficients a,b,c
*
*       ATTENTION: FOR SMALL |a|, this function is 
*       liable to produce garbage!
*
*                Input parameters:
*
*  a,b,c - coefficients of the equation
*  
*                Output parameters:
*
*  ierp - error return code: 
*            ierp=0 means successful conclusion
*            ierp=16 means that the equation has no real
*                   roots
*  x1p, x2p - solutions of the equation
*  err1p, err2p discrepancies corresponding to the 
*          solutions x1p, x2p 
*
*/
        *ierp=0;

        discr=b*b-4*a*c;

        if(discr < 0) {
            *ierp=16;
            return;
        }


        x1=(-b-sqrt(discr)) /(2*a);

        x2=(-b+sqrt(discr)) /(2*a);

//        cprin2("x1=",&x1,1);
//        cprin2("x2=",&x2,1);

        err1=a*x1*x1+b*x1+c;
        err2=a*x2*x2+b*x2+c;

//        cprin2("and err1=",&err1,1);
//        cprin2("and err2=",&err2,1);

        *x1p=x1;
        *x2p=x2;
        *err1p=err1;
        *err2p=err2;

}