c*********************************************************************
c Filename:      quadratic.f
c Version:       2.0
c Description:   Simple program to demonstate use of quadratic solver
c Author:        Eddie Baron <baron@nhn.ou.edu>
c Created at:    Mon Aug 25 13:14:00 2003
c Modified at:   Wed Aug 27 15:33:26 2003
c Modified by:   Eddie Baron <baron@nhn.ou.edu>
c*********************************************************************
c
      program quadratic
c     -----------------
      implicit none
      integer,parameter :: sp=kind(0.0)
      real (kind=sp) :: x1,x2,a,b,c,discrim
c
      write(*,*)"y=a*x^2 + b*x +c = 0"
      write(*,*)"give a,b,c: "
      read(*,*)a,b,c
c
c
c      put equation in standard form
c      x^2 + 2bx + c
c
      if(a .eq. 0.d0) then
       write(*,*)"this is not a quadratic equation"
       call stop_exit(1,"this is not a quadratic equation")
      endif
      
      b = b/(2.d0*a)
      c = c/a
c--      
c-- compute the roots from the standard formula
c--
      discrim = sqrt(b**2 - c)
      if(discrim .ne. discrim) then
       write(*,*)"roots are imaginary"
       call stop_exit(1,"roots are imaginary")
      endif
      x1 = -b + discrim
      x2 = -b - discrim
      print *,"standard method gives x1 = ",x1," x2 = ",x2
c--      
c-- compute the roots using the larger root
c--
      x1 = -b - sign(discrim,b)
      x2 = c/x1
      print *,"accurate  method gives x1 = ",x1," x2 = ",x2
      call stop_exit(0,"sucessful completion")
      end
c
c
c
      subroutine stop_exit(status,message)
c     ------------------------------------
      integer status
      character(len=*) :: message
************************************************************************
*  this subroutine is used to print a message to stderr and 
*  stop the code with a return value, if possible.
*        version 1.0 of 20/may/94 by eab
*-- notes:
*  this routine is machine dependent (call to exit-routine only)
*  and this version is by far not complete ...
* exit_() is for IBM RS/6000
************************************************************************
      write(0,'(a)')message
c      call exit_(status) ! use on IBMs
      call exit(status) ! use on SUNs
c     stop
      return
      end