/***********************************************************************/
/*                                                                     */
/*              Finding mt2 by Bisection                               */
/*                                                                     */
/*              Authors: Hsin-Chia Cheng, Zhenyu Han                   */ 
/*              Dec 11, 2008, v1.01a                                     */
/*                                                                     */  
/***********************************************************************/


/*******************************************************************************
  Usage: 

  1. Define an object of type "mt2":
     
     mt2_bisect::mt2 mt2_event;
 
  2. Set momenta and the mass of the invisible particle, mn:
 
     mt2_event.set_momenta( pa, pb, pmiss );
     mt2_event.set_mass( mn );
 
     where array pa[0..2], pb[0..2], pmiss[0..2] contains (mass,px,py) 
     for the visible particles and the missing momentum. pmiss[0] is not used. 
     All quantities are given in double.    

  3. Use mt2::get_mt2() to obtain the value of mt2:

     double mt2_value = mt2_event.get_mt2();       
          
*******************************************************************************/ 
              
#include <iostream>
#include <math.h>
#include "mt2_bisect.h"

using namespace std;

namespace mt2_bisect
{
mt2::mt2()
{
   solved = false;
   momenta_set = false;
   mt2_b  = 0.;
   scale = 1.;
}

double mt2::get_mt2()
{
   if (!momenta_set)
   {
       cout <<" Please set momenta first!" << endl;
       return 0;
   }
        
   if (!solved) mt2_bisect();
   return mt2_b*scale;
}

void mt2::set_momenta(double* pa0, double* pb0, double* pmiss0)
{
   solved = false;     //reset solved tag when momenta are changed.
   momenta_set = true;

   ma = fabs(pa0[0]);  // mass cannot be negative

   if (ma < ZERO_MASS) ma = ZERO_MASS;

   pax  = pa0[1]; 
   pay  = pa0[2];
   masq = ma*ma;
   Easq = masq+pax*pax+pay*pay;
   Ea   = sqrt(Easq);
   
   mb = fabs(pb0[0]);

   if (mb < ZERO_MASS) mb = ZERO_MASS;

   pbx  = pb0[1]; 
   pby  = pb0[2];
   mbsq = mb*mb;
   Ebsq = mbsq+pbx*pbx+pby*pby;
   Eb   = sqrt(Ebsq);
     
   pmissx   = pmiss0[1]; pmissy = pmiss0[2];
   pmissxsq = pmissx*pmissx;
   pmissysq = pmissy*pmissy;

// set ma>= mb
   if(masq < mbsq)
   {
      double temp;
      temp = pax;  pax  = pbx;  pbx  = temp;
      temp = pay;  pay  = pby;  pby  = temp;
      temp = Ea;   Ea   = Eb;   Eb   = temp;
      temp = Easq; Easq = Ebsq; Ebsq = temp;
      temp = masq; masq = mbsq; mbsq = temp;
      temp = ma;   ma   = mb;   mb   = temp;   
   }
//normalize max{Ea, Eb} to 100
   if (Ea > Eb) scale = Ea/100.;
   else scale = Eb/100.;
   
   if (sqrt(pmissxsq+pmissysq)/100 > scale) scale = sqrt(pmissxsq+pmissysq)/100;
   //scale = 1;
   double scalesq = scale * scale;
   ma  = ma/scale;
   mb  = mb/scale;
   masq = masq/scalesq;
   mbsq = mbsq/scalesq;
   pax = pax/scale; pay = pay/scale;
   pbx = pbx/scale; pby = pby/scale;
   Ea  = Ea/scale;  Eb = Eb/scale;
   
   Easq = Easq/scalesq;
   Ebsq = Ebsq/scalesq;
   pmissx = pmissx/scale;
   pmissy = pmissy/scale;
   pmissxsq = pmissxsq/scalesq;
   pmissysq = pmissysq/scalesq;
   mn   = mn_unscale/scale; 
   mnsq = mn*mn;
  
   if (ABSOLUTE_PRECISION > 100.*RELATIVE_PRECISION) precision = ABSOLUTE_PRECISION;
   else precision = 100.*RELATIVE_PRECISION;
}

void mt2::set_mn(double mn0)
{
   solved = false;    //reset solved tag when mn is changed.
   mn_unscale   = fabs(mn0);  //mass cannot be negative
   mn = mn_unscale/scale;
   mnsq = mn*mn;
}

void mt2::print()
{
   cout << " pax = " << pax*scale << ";   pay = " << pay*scale << ";   ma = " << ma*scale <<";"<< endl;
   cout << " pbx = " << pbx*scale << ";   pby = " << pby*scale << ";   mb = " << mb*scale <<";"<< endl;
   cout << " pmissx = " << pmissx*scale << ";   pmissy = " << pmissy*scale <<";"<< endl;
   cout << " mn = " << mn_unscale<<";" << endl;
}

//special case, the visible particle is massless
void mt2::mt2_massless()
{
   
//rotate so that pay = 0 
   double theta,s,c;
   theta = atan(pay/pax);
   s = sin(theta);
   c = cos(theta);

   double pxtemp,pytemp;
   Easq   = pax*pax+pay*pay;
   Ebsq   = pbx*pbx+pby*pby;
   Ea     = sqrt(Easq);
   Eb     = sqrt(Ebsq);
  
   pxtemp = pax*c+pay*s;
   pax    = pxtemp;
   pay    = 0;
   pxtemp = pbx*c+pby*s;
   pytemp = -s*pbx+c*pby;
   pbx    = pxtemp;
   pby    = pytemp;
   pxtemp = pmissx*c+pmissy*s;
   pytemp = -s*pmissx+c*pmissy;
   pmissx = pxtemp;
   pmissy = pytemp;

   a2  = 1-pbx*pbx/(Ebsq);
   b2  = -pbx*pby/(Ebsq);
   c2  = 1-pby*pby/(Ebsq);

   d21 = (Easq*pbx)/Ebsq;
   d20 = - pmissx +  (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
   e21 = (Easq*pby)/Ebsq;
   e20 = - pmissy +  (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
   f22 = -(Easq*Easq/Ebsq);
   f21 = -2*Easq*(pbx*pmissx + pby*pmissy)/Ebsq;
   f20 = mnsq + pmissxsq + pmissysq - (pbx*pmissx + pby*pmissy)*(pbx*pmissx + pby*pmissy)/Ebsq;

   double Deltasq0    = 0; 
   double Deltasq_low, Deltasq_high;
   int    nsols_high, nsols_low;

   Deltasq_low = Deltasq0 + precision;
   nsols_low = nsols_massless(Deltasq_low);
   
   if(nsols_low > 1) 
   { 
      mt2_b = (double) sqrt(Deltasq0+mnsq);
      return;
   }

/*   
   if( nsols_massless(Deltasq_high) > 0 )
   {
      mt2_b = (double) sqrt(mnsq+Deltasq0);
      return;
   }*/

//look for when both parablos contain origin  
   double Deltasq_high1, Deltasq_high2;
   Deltasq_high1 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy;
   Deltasq_high2 = 2*Ea*mn;
 
   if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high2;
     else Deltasq_high = Deltasq_high1;

   nsols_high=nsols_massless(Deltasq_high);
  
   int foundhigh;
   if (nsols_high == nsols_low)
   {
      
      
      foundhigh=0;
      
      double minmass, maxmass;
      minmass  = mn ;
      maxmass  = sqrt(mnsq + Deltasq_high);
      for(double mass = minmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
      {
	 Deltasq_high = mass*mass - mnsq;
      
         nsols_high = nsols_massless(Deltasq_high);
         if(nsols_high>0)
         {
            foundhigh=1;
            Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
            break;
         }
      }
      if(foundhigh==0) 
      {
       
	cout<<"Deltasq_high not found at event " << nevt <<endl;
        
       
         mt2_b = (double)sqrt(Deltasq_low+mnsq);
         return;
      }
   }

   if(nsols_high == nsols_low)
   { 
      cout << "error: nsols_low=nsols_high=" << nsols_high << endl;
      cout << "Deltasq_high=" << Deltasq_high << endl;
      cout << "Deltasq_low= "<< Deltasq_low << endl;
    
      mt2_b = sqrt(mnsq + Deltasq_low);
      return;
   }
   double minmass, maxmass;
   minmass = sqrt(Deltasq_low+mnsq);
   maxmass = sqrt(Deltasq_high+mnsq);
   while(maxmass - minmass > precision)
   {
      double Delta_mid, midmass, nsols_mid;
      midmass   = (minmass+maxmass)/2.;
      Delta_mid = midmass * midmass - mnsq;
      nsols_mid = nsols_massless(Delta_mid);
      if(nsols_mid != nsols_low) maxmass = midmass;
      if(nsols_mid == nsols_low) minmass = midmass;
   }
   mt2_b = minmass;
   return;
}

int mt2::nsols_massless(double Dsq)
{
  double delta;
  delta = Dsq/(2*Easq);
  d1    = d11*delta;
  e1    = e11*delta;
  f1    = f12*delta*delta+f10;
  d2    = d21*delta+d20;
  e2    = e21*delta+e20;
  f2    = f22*delta*delta+f21*delta+f20;
  
  double a,b;
  if (pax > 0) a = Ea/Dsq;
  else         a = -Ea/Dsq;
  if (pax > 0) b = -Dsq/(4*Ea)+mnsq*Ea/Dsq;
  else         b = Dsq/(4*Ea)-mnsq*Ea/Dsq;
  
  double A4,A3,A2,A1,A0;

  A4 = a*a*a2;
  A3 = 2*a*b2/Ea;
  A2 = (2*a*a2*b+c2+2*a*d2)/(Easq);
  A1 = (2*b*b2+2*e2)/(Easq*Ea);
  A0 = (a2*b*b+2*b*d2+f2)/(Easq*Easq);
  
  long  double A0sq, A1sq, A2sq, A3sq, A4sq;
  A0sq = A0*A0;
  A1sq = A1*A1;
  A2sq = A2*A2;
  A3sq = A3*A3;
  A4sq = A4*A4;

  long double B3, B2, B1, B0;
  B3 = 4*A4;
  B2 = 3*A3;
  B1 = 2*A2;
  B0 = A1;
  long double C2, C1, C0;
  C2 = -(A2/2 - 3*A3sq/(16*A4));
  C1 = -(3*A1/4. -A2*A3/(8*A4));
  C0 = -A0 + A1*A3/(16*A4);
  long double  D1, D0;
  D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
  D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
 
  long double E0;
  E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
 
  long  double t1,t2,t3,t4,t5;
   
//find the coefficients for the leading term in the Sturm sequence  
   t1 = A4;
   t2 = A4;
   t3 = C2;
   t4 = D1;
   t5 = E0;
  
   int nsol;
   nsol = signchange_n(t1,t2,t3,t4,t5)-signchange_p(t1,t2,t3,t4,t5);
   if( nsol < 0 ) nsol=0;

   return nsol;
  
}

void mt2::mt2_bisect()
{
  
   
   solved = true;
   cout.precision(11);

//if masses are very small, use code for massless case.  
   if(masq < MIN_MASS && mbsq < MIN_MASS) 
   { 
      mt2_massless();
      return;
   }
 

   double Deltasq0;     
   Deltasq0 = ma*(ma + 2*mn); //The minimum mass square to have two ellipses 
 
// find the coefficients for the two quadratic equations when Deltasq=Deltasq0.
  
   a1 = 1-pax*pax/(Easq);
   b1 = -pax*pay/(Easq);
   c1 = 1-pay*pay/(Easq);
   d1 = -pax*(Deltasq0-masq)/(2*Easq);
   e1 = -pay*(Deltasq0-masq)/(2*Easq);
   a2 = 1-pbx*pbx/(Ebsq);
   b2 = -pbx*pby/(Ebsq);
   c2 = 1-pby*pby/(Ebsq);
   d2 = -pmissx+pbx*(Deltasq0-mbsq)/(2*Ebsq)+pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
   e2 = -pmissy+pby*(Deltasq0-mbsq)/(2*Ebsq)+pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
   f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq0-mbsq)/(2*Eb)+
        (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq0-mbsq)/(2*Eb)+
        (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
   
// find the center of the smaller ellipse 
   double x0,y0;
   x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
   y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);

   
// Does the larger ellipse contain the smaller one? 
   double dis=a2*x0*x0+2*b2*x0*y0+c2*y0*y0+2*d2*x0+2*e2*y0+f2;

   if(dis<=0.01)
   { 
      mt2_b  = (double) sqrt(mnsq+Deltasq0);
      return;
   }
   

/* find the coefficients for the two quadratic equations           */
/* coefficients for quadratic terms do not change                  */
/* coefficients for linear and constant terms are polynomials of   */
/*       delta=(Deltasq-m7sq)/(2 E7sq)                             */  
   d11 = -pax;
   e11 = -pay;
   f10 = mnsq;
   f12 = -Easq;
   d21 = (Easq*pbx)/Ebsq;
   d20 = ((masq - mbsq)*pbx)/(2.*Ebsq) - pmissx +
         (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
   e21 = (Easq*pby)/Ebsq;
   e20 = ((masq - mbsq)*pby)/(2.*Ebsq) - pmissy +
         (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
   f22 = -Easq*Easq/Ebsq;
   f21 = (-2*Easq*((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb))/Eb;
   f20 = mnsq + pmissx*pmissx + pmissy*pmissy - 
         ((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb)
         *((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb);

//Estimate upper bound of mT2
//when Deltasq > Deltasq_high1, the larger encloses the center of the smaller 
   double p2x0,p2y0;
   double Deltasq_high1;
   p2x0 = pmissx-x0;
   p2y0 = pmissy-y0;
   Deltasq_high1 = 2*Eb*sqrt(p2x0*p2x0+p2y0*p2y0+mnsq)-2*pbx*p2x0-2*pby*p2y0+mbsq;
   
//Another estimate, if both ellipses enclose the origin, Deltasq > mT2

   double Deltasq_high2, Deltasq_high21, Deltasq_high22;
   Deltasq_high21 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy+mbsq;
   Deltasq_high22 = 2*Ea*mn+masq;
  
   if ( Deltasq_high21 < Deltasq_high22 ) Deltasq_high2 = Deltasq_high22;
   else Deltasq_high2 = Deltasq_high21;

//pick the smaller upper bound   
   double Deltasq_high;
   if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high1;
   else Deltasq_high = Deltasq_high2;
   
  
   double Deltasq_low; //lower bound
   Deltasq_low = Deltasq0;

//number of solutions at Deltasq_low should not be larger than zero
   if( nsols(Deltasq_low) > 0 )
   {
     //cout << "nsolutions(Deltasq_low) > 0"<<endl;
     mt2_b = (double) sqrt(mnsq+Deltasq0);
     return;
   }
  
   int nsols_high, nsols_low;

   nsols_low  = nsols(Deltasq_low);
   int foundhigh;
  

//if nsols_high=nsols_low, we missed the region where the two ellipse overlap 
//if nsols_high=4, also need a scan because we may find the wrong tangent point.

   nsols_high = nsols(Deltasq_high);
  
   if(nsols_high == nsols_low || nsols_high == 4)
   {
      //foundhigh = scan_high(Deltasq_high);
      foundhigh = find_high(Deltasq_high);
      if(foundhigh == 0) 
      {
 	 cout << "Deltasq_high not found at event " << nevt << endl;
         mt2_b = sqrt( Deltasq_low + mnsq );
         return;
      }
      
   }

   while(sqrt(Deltasq_high+mnsq) - sqrt(Deltasq_low+mnsq) > precision)
   {
      double Deltasq_mid,nsols_mid;
      //bisect
      Deltasq_mid = (Deltasq_high+Deltasq_low)/2.;
      nsols_mid = nsols(Deltasq_mid);
      // if nsols_mid = 4, rescan for Deltasq_high
      if ( nsols_mid == 4 ) 
      {
         Deltasq_high = Deltasq_mid;
         //scan_high(Deltasq_high);
         find_high(Deltasq_high);
         continue;
      } 
         
      
      if(nsols_mid != nsols_low) Deltasq_high = Deltasq_mid;
      if(nsols_mid == nsols_low) Deltasq_low  = Deltasq_mid;
   }
   mt2_b = (double) sqrt( mnsq + Deltasq_high);
   return;
}

int mt2::find_high(double & Deltasq_high)
{
   double x0,y0;
   x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
   y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
   double Deltasq_low = (mn + ma)*(mn + ma) - mnsq;
   do 
   {
      double Deltasq_mid = (Deltasq_high + Deltasq_low)/2.;
      int nsols_mid = nsols(Deltasq_mid);
      if ( nsols_mid == 2 )
      {
         Deltasq_high = Deltasq_mid;
         return 1;
      }
      else if (nsols_mid == 4)
      {
         Deltasq_high = Deltasq_mid;
         continue;
      }
      else if (nsols_mid ==0)
      {
         d1 = -pax*(Deltasq_mid-masq)/(2*Easq);
         e1 = -pay*(Deltasq_mid-masq)/(2*Easq);
         d2 = -pmissx + pbx*(Deltasq_mid - mbsq)/(2*Ebsq)
              + pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
         e2 = -pmissy + pby*(Deltasq_mid - mbsq)/(2*Ebsq)
              + pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
         f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq_mid-mbsq)/(2*Eb)+
              (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq_mid-mbsq)/(2*Eb)+
              (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
// Does the larger ellipse contain the smaller one? 
         double dis = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;
         if (dis < 0) Deltasq_high = Deltasq_mid;
           else Deltasq_low = Deltasq_mid;
      }
      
   } while ( Deltasq_high - Deltasq_low > 0.001);
   return 0;
}  
int mt2::scan_high(double & Deltasq_high)
{
   int foundhigh = 0 ;
   int nsols_high;

   
   double Deltasq_low;
   double tempmass, maxmass;
   tempmass = mn + ma;
   maxmass  = sqrt(mnsq + Deltasq_high);
   if (nevt == 32334) cout << "Deltasq_high = " << Deltasq_high << endl;
   for(double mass = tempmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
   {
      Deltasq_high = mass*mass - mnsq;
      nsols_high   = nsols(Deltasq_high);
      
      if( nsols_high > 0)
      {
	 Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
         foundhigh   = 1;
         break;
      }
    }
    return foundhigh;
}
int mt2::nsols(  double Dsq)
{
   double delta = (Dsq-masq)/(2*Easq);
  
//calculate coefficients for the two quadratic equations
   d1 = d11*delta;
   e1 = e11*delta;
   f1 = f12*delta*delta+f10;
   d2 = d21*delta+d20;
   e2 = e21*delta+e20;
   f2 = f22*delta*delta+f21*delta+f20;

//obtain the coefficients for the 4th order equation 
//devided by Ea^n to make the variable dimensionless
   long double A4, A3, A2, A1, A0;

   A4 = 
   -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 + 
   4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 + 
   a1*a1*c2*c2;  

   A3 =
     (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 + 
   8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 + 
   4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 - 
     4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Ea;

   
   A2 =
     (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 - 
   8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 + 
   4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 - 
   8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 + 
   4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 + 
     4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Easq;
  
   A1 =
     (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 - 
   4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 - 
   4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 - 
     4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Easq*Ea);
  
   A0 =
     (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 + 
   4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 + 
     a1*a1*f2*f2)/(Easq*Easq);
   
   long  double A0sq, A1sq, A2sq, A3sq, A4sq;
   A0sq = A0*A0;
   A1sq = A1*A1;
   A2sq = A2*A2;
   A3sq = A3*A3;
   A4sq = A4*A4;
   
   long double B3, B2, B1, B0;
   B3 = 4*A4;
   B2 = 3*A3;
   B1 = 2*A2;
   B0 = A1;
   
   long double C2, C1, C0;
   C2 = -(A2/2 - 3*A3sq/(16*A4));
   C1 = -(3*A1/4. -A2*A3/(8*A4));
   C0 = -A0 + A1*A3/(16*A4);
   
   long double D1, D0;
   D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
   D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
   
   long double E0;
   E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
   
   long  double t1,t2,t3,t4,t5;
//find the coefficients for the leading term in the Sturm sequence  
   t1 = A4;
   t2 = A4;
   t3 = C2;
   t4 = D1;
   t5 = E0;
 

//The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
   int nsol;
   nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5);

//Cannot have negative number of solutions, must be roundoff effect
   if (nsol < 0) nsol = 0;

   return nsol;
  
}  

inline int mt2::signchange_n( long double t1, long double t2, long double t3, long double t4, long double t5)
{
   int nsc;
   nsc=0;
   if(t1*t2>0) nsc++;
   if(t2*t3>0) nsc++;
   if(t3*t4>0) nsc++;
   if(t4*t5>0) nsc++;
   return nsc;
}
inline int mt2::signchange_p( long double t1, long double t2, long double t3, long double t4, long double t5)
{
   int nsc;
   nsc=0;
   if(t1*t2<0) nsc++;
   if(t2*t3<0) nsc++;
   if(t3*t4<0) nsc++;
   if(t4*t5<0) nsc++;
   return nsc;
}

}//end namespace mt2_bisect