Namespaces

namespace rdslv
namespace sbdrl

Classes

struct method_base
struct select_mth
struct select_mth< system, 0 >
struct select_mth< system, 1 >
struct select_mth< system, 3 >
struct select_mth< system, 4 >
struct select_mth< system, 2 >
struct select_mth< system, 5 >
struct select_mth< system, 6 >
struct select_mth< system, 7 >
class solver_mv_monomial
class box_rep
Box representation. More...
struct method_debug
struct method
struct system_ctrl

Functions

template<class T , class eenv_t , class MPOL > void fill_data (T *data, eenv_t *env, const MPOL &p, binomials< T > &benv)
template<class monom > void scan_monom_env (std::set< int > &env, const monom &m)
template<class POL , class FT > int topological_degree_2d (box_rep< POL > *b)
template<class POL > polynomial< double, with
< MonomialTensor > > * jacobian (Seq< POL * > &S0)
template<class T > void precondition (Seq< T > val, polynomial< double, with< MonomialTensor > > *mat, int n)
template<class T > T * eval_poly_matrix (Seq< T > val, polynomial< double, with< MonomialTensor > > *mat, int n)
template<class T > int signof_det (T *mat, const int &n)
template<class T > T * submatrix (T *matrix, int order, int i, int j)
template<class T > T det (T *mat, int order)
template<class FT > int signof (polynomial< double, with< MonomialTensor > > *p, Seq< FT > t)
template<class POL > bool miranda_test (Seq< POL > &S, const int i, const int j)
template<class POL > bool include1 (box_rep< POL > *b, polynomial< double, with< MonomialTensor > > *J)
Inclusion criteria (Miranda Test)
template<class POL > bool include2 (box_rep< POL > *b, polynomial< double, with< MonomialTensor > > *J)
Inclusion criteria (Jacobian+Topological Degree)
template<class POL > bool include3 (box_rep< POL > *b, polynomial< double, with< MonomialTensor > > *J, Seq< POL * > &S0)
Inclusion criteria based on Rump's test.
template<class POL > bool exclude1 (box_rep< POL > *b, Seq< POL * > &S0)
Exclusion criteria (inteval arithmetic)
template<class POL > bool exclude2 (box_rep< POL > *b, polynomial< double, with< MonomialTensor > > *J)
Exclusion criteria (topological degree)
template<class POL > bool exclude3 (box_rep< POL > *b)
template<class POL > bool no_variation (POL &p)
True iff all coefficients of p have the same sign.
template<class FT > int sgn (FT a)
Sign of a.
template<class POL , class FT > int sgn (POL &f, FT a)
Sign of f at point a>=0 (-1=inf)
template<class real_t > unsigned clean_result (real_t *isols, int nvars, int nsols, const real_t &prec)

Variables

const int E_CTRL = 0
const int E_STRGY = 1
const int E_RDSLV = 2
const int E_SBDRL = 3
const int E_INIT = 4
const int C_ACCEPT = 5
const int R_REJECT = 6
const int R_ISOK = 7
const int R_WEAK = 8
const int R_FAIL = 9
const int R_ERROR = 10
const int D_REJECT = 11

Function Documentation

unsigned mmx::realroot::clean_result	(	real_t *	isols,
		int	nvars,
		int	nsols,
		const real_t &	prec
	)

Definition at line 34 of file system_support.hpp.

References mmx::abs(), mmx::sparse::copy(), and mmx::max().

Referenced by solver< C, ProjRd< MTH > >::run(), and solver< C, ProjRd< MTH > >::solve_monomial().

  {
    int a,c,n,i;
    real_t cd;
    if ( nsols == 0 ) return 0;
    for ( a = 0; a < nsols*nvars; isols[a] = (isols[2*a]+isols[2*a+1])/2.0, a ++ ) ;

    for ( n = 1, c = 0, a = 0; a < nsols*nvars; a += nvars )
      {
        for ( cd = 0, i = 0; i < nvars; cd = std::max(cd,std::abs((isols+c)[i]-(isols+a)[i])), i ++ ) ;
        if ( cd < prec )
        {
          double fa = ((double)n)/(n+1);
          double fb = ((double)1.0)/(n+1);
          for ( int i = 0; i < nvars; i ++ ) (isols+c)[i] = fa*(isols+c)[i]+fb*(isols+a)[i];
          n ++;
        }
        else
          {
            c += nvars;
            std::copy(isols+a,isols+a+nvars,isols+c);
            n  = 1;
          };
      };
    return (c+nvars)/nvars;
  };

T mmx::realroot::det	(	T *	mat,
		int	order
	)

Definition at line 136 of file solver_mv_monomial_tests.hpp.

References submatrix().

Referenced by exclude2(), include1(), include2(), and signof_det().

    {
      if(order == 1)
        return mat[0]; // Matrix is of order one
      int i;
      T ev(0);

      for(i = 0; i < order; i++)
      { 
        //std::cout<<"entering"<<order<<",i="<<i <<std::endl;
        ev +=  T(i%2==0?1:-1)* mat[i*order] *  det(submatrix(mat, order, i, 0), order - 1);

        //std::cout<<"OK, ("<<order<<" "<<i<<")"<<std::endl;
      }
      return ev;
    }

T* mmx::realroot::eval_poly_matrix	(	Seq< T >	val,
		polynomial< double, with< MonomialTensor > > *	mat,
		int	n
	)		`[inline]`

Definition at line 84 of file solver_mv_monomial_tests.hpp.

References mmx::eval(), mmx::rep(), and Seq< C, R >::size().

Referenced by exclude2(), include1(), include2(), and include3().

    {
      int m(val.size());
      T *ev;//[m*n];
      ev= new T [m*n];

      for( int i=0;i<m;i++)
      {
        for( int k=0;k<n;k++)
        {
          tensor::eval( ev[i*m+k] , mat[i*m+k].rep(),
                        val , m  );
        }
      }
      return ev;
    }

bool mmx::realroot::exclude1	(	box_rep< POL > *	b,
		Seq< POL * > &	S0
	)

Exclusion criteria (inteval arithmetic)

Definition at line 344 of file solver_mv_monomial_tests.hpp.

References mmx::assign(), DPOL, mmx::eval(), Interval< T, r >::M, and Interval< T, r >::m.

    {
      Interval<double> ev;
      Seq<Interval<double> > dom;
      dom= b->template domain<double>();
      DPOL tmp;

      for (unsigned i=0; i!=b->nbpol(); ++i)
      {
        let::assign(tmp, *S0[i] );
        tensor::eval( ev , tmp.rep(), 
                      dom , b->nbvar() );       
        if ( ev.m*ev.M > 0 )
        {   
          //   std::cout<<i<<" ev"<< 
          //    dom<<"\nf= "<< *(S0[1])<<
          //   "\n f([])= " <<ev<<std::endl;
          //  if (td!=0)
          // std::cout<<"!!!!!!----td: "<<td <<std::endl;
          return true;
        }       
        
      }
      return false;
    };

bool mmx::realroot::exclude2	(	box_rep< POL > *	b,
		polynomial< double, with< MonomialTensor > > *	J
	)

Exclusion criteria (topological degree)

Definition at line 373 of file solver_mv_monomial_tests.hpp.

References det(), eval_poly_matrix(), Interval< T, r >::M, and Interval< T, r >::m.

Referenced by solver_mv_monomial< FT, POL >::solve_system().

  {
    Interval<double> ev;
    Interval<double> * m;
    
    //Seq<POL> SS= b->system();
    unsigned dim= b->nbvar();
    
    // Check Jacobian's sign
    m= eval_poly_matrix( b->template domain<double>() ,J,dim);
    ev= det(m,dim);
    if ( ev.m*ev.M < 0 ) 
      return false;

    int td= topological_degree_2d<POL, double>(b);
    
    return ( td==0 );
  }

bool mmx::realroot::exclude3 ( box_rep< POL > * b ) [inline]

Definition at line 393 of file solver_mv_monomial_tests.hpp.

References no_variation().

    {
      Seq<POL> S= b->system();
      for (unsigned i=0; i!=b->nbpol(); ++i)
        if ( no_variation( S[i]) ) 
          return true;
      return false;
    }

void mmx::realroot::fill_data	(	T *	data,
		eenv_t *	env,
		const MPOL &	p,
		binomials< T > &	benv
	)

Definition at line 82 of file solver_mv_bernstein.hpp.

References mmx::tensor::convertm2b(), mmx::vct::fill(), and mmx::tensor::mpolfill().

Referenced by solver< C, ProjRd< MTH > >::solve_monomial().

  {
    std::fill(data,data+env->sz(),0); 
    mpolfill(data,p.rep(),*env);
    convertm2b (data,*env,benv);
  };

bool mmx::realroot::include1	(	box_rep< POL > *	b,
		polynomial< double, with< MonomialTensor > > *	J
	)

Inclusion criteria (Miranda Test)

Definition at line 203 of file solver_mv_monomial_tests.hpp.

References det(), eval_poly_matrix(), Interval< T, r >::M, Interval< T, r >::m, miranda_test(), and POL.

    {
      Interval<double> ev(0);
      Interval<double> * m;
      unsigned i,j,c; //,r ;
      POL u,l;
      Seq<POL> SS= b->system();
      unsigned dim= b->nbvar();

      bool t[dim][dim];

      m= eval_poly_matrix( b->template domain<double>() ,J,dim);

// std::cout<<"Ok " <<std::endl;     
// double dm[9];
// for(i = 0; i < 9; i++)
// {  dm[i]= m[i].M;
//   std::cout<<"dm["<<i<<"]="<<dm[i] <<std::endl;
// }
// ev= det(dm,dim);
// std::cout<<"Ok ="<<ev <<std::endl;     

// for(i = 0; i < 9; i++)
//   std::cout<<"m["<<i<<"]="<<m[i] <<std::endl;

 //   std::cout<<"SUB "<<std::endl;
 // ev= det(submatrix(m, 3, 2, 0), 2);
 //   std::cout<<"FULL "<<std::endl;

      ev=det(m,dim);

      if ( ev.m*ev.M < 0 ) 
        return false;
      
      for (i=0; i!=dim;++i)
        for (j=0; j!=dim;++j){
          t[i][j]= miranda_test(SS,i,j);
          //std::cout<<"mir("<<i<<","<<j<<")="<< t[i][j] <<std::endl;
        }

      c=0;
      for (i=0; i!=dim;++i)
        for (j=0; j!=dim;++j)
          if (t[i][j]==true) 
          {     
            c++; break;
          }
      if (c<dim) return false;
      
      c=0;
      for (i=0; i!=dim;++i)
        for (j=0; j!=dim;++j)
          if (t[j][i]==true) 
          {     
            c++; break;
          }
      if (c<dim) return false;
      
      //std::cout<<"INCLUDE. ev="<<ev<<", c="<<c <<"*dim:"<<dim <<std::endl;
      return true;
    };

bool mmx::realroot::include2	(	box_rep< POL > *	b,
		polynomial< double, with< MonomialTensor > > *	J
	)

Inclusion criteria (Jacobian+Topological Degree)

Definition at line 267 of file solver_mv_monomial_tests.hpp.

References det(), eval_poly_matrix(), Interval< T, r >::M, and Interval< T, r >::m.

  {
    Interval<double> ev;
    Interval<double> * m;
    
    //Seq<POL> SS= b->system();
    unsigned dim= b->nbvar();
    
    // Check Jacobian's sign
    m= eval_poly_matrix( b->template domain<double>() ,J,dim);
    ev= det(m,dim);
    if ( ev.m*ev.M < 0 ) 
      return false;

    int td= topological_degree_2d<POL, double>(b);

    //if ( (td==-1 || td==1) )
    //{ std::cout<<"INCLUDE. ev="<<ev<<", td="<<td <<std::endl;
    //this->print();  }
    
    return ( (td==-1 || td==1) );
  }

bool mmx::realroot::include3	(	box_rep< POL > *	b,
		polynomial< double, with< MonomialTensor > > *	J,
		Seq< POL * > &	S0
	)

Inclusion criteria based on Rump's test.

Definition at line 292 of file solver_mv_monomial_tests.hpp.

References mmx::assign(), DPOL, mmx::eval(), eval_poly_matrix(), mmx::linear::LUinverse(), Interval< T, r >::M, and Interval< T, r >::m.

        {
          double *fc, *jc, *ijc;
          Interval<double> ev, *jb, *m;
          Seq<double> c= b->point( b->middle() ) ;
          Seq<Interval<double> > bx;
          unsigned i,j,n= b->nbvar();
          DPOL tmp;

          fc= new double[n];                // system S0 evaluated on c
          for(i=0; i< n;i++ )
          {
            let::assign(tmp, *S0[i]) ;
            tensor::eval( fc[i] , tmp.rep(), c , n  );
          }
          bx= b->template domain<double>(); // box 
          jc= eval_poly_matrix( c, J, n);   // Jacobian evaluated on c
          ijc= new double[n*n];
          linear::LUinverse(ijc, jc, n);     // Inverse of above
          jb= eval_poly_matrix( bx, J, n);  // Jacobian evaluated on box
          
          m= new Interval<double>[n*n];
          for(i=0; i< n;i++ )
            for(j=0; j< n;j++ )
              m[i*n+j] += ijc[i*n+j]*jb[j*n+i];
          free(jb);

          // check if -im*fc + (I-im*jb)*bx is inside bx.
          
          for(i=0; i< n;i++ )
          {
            bx[i]= bx[i]- c[i];
            ev= bx[i];
            for(j=0; j<n; j++ )
            {
              ev+= -ijc[i*n+j]*fc[j] - m[i*n+j] *bx[i];
            }
//            std::cout<<ev<<" --- "<<bx[i]<<std::endl;
            if (  ev.m < bx[i].m ||
                  ev.M > bx[i].M )
            {
//              std::cout<<ev<<" not in "<<bx[i]<<std::endl;
                  return false;
            }
          }
//          std::cout<<"found root in "<<bx<< "(ev="<<ev<<")"<<std::endl;
            return ( true );
        }

polynomial<double,with<MonomialTensor> >* mmx::realroot::jacobian ( Seq< POL * > & S0 ) [inline]

Definition at line 42 of file solver_mv_monomial_tests.hpp.

References mmx::assign(), mmx::diff(), DPOL, POL, and Seq< C, R >::size().

Referenced by solver_mv_monomial< FT, POL >::approximate(), and solver_mv_monomial< FT, POL >::isolate().

    {
      int n=S0.size();
      DPOL * jac;//[n*n]; 
      POL df;
      jac= new DPOL [n*n];
      
      for( int i=0;i<n;i++)
      {
        for( int k=0;k<n;k++)
        {
          df= diff( *(S0[i]),k);
          let::assign( jac[i*n+k] , df);
        }
      
      }
      return jac;
    }

bool mmx::realroot::miranda_test	(	Seq< POL > &	S,
		const int	i,
		const int	j
	)		`[inline]`

Definition at line 187 of file solver_mv_monomial_tests.hpp.

References mmx::tensor::face(), no_variation(), and POL.

Referenced by include1().

  {
    POL u,l;
    
    tensor::face(l, S[i], j, 0);
    tensor::face(u, S[i], j, 1);
    //std::cout<<"l="<<l<<",u="<<u<<std::endl;

    return  ( no_variation(l)    && 
              no_variation(u)    &&
              (l[0]>0)!=(u[0]>0) );
  };

bool mmx::realroot::no_variation ( POL & p ) [inline]

True iff all coefficients of p have the same sign.

Definition at line 404 of file solver_mv_monomial_tests.hpp.

References mmx::has_sign_variation().

Referenced by exclude3(), miranda_test(), and box_rep< POL >::miranda_test().

  {
    return !has_sign_variation(p.begin(),p.end()) ;
  };

void mmx::realroot::precondition	(	Seq< T >	val,
		polynomial< double, with< MonomialTensor > > *	mat,
		int	n
	)		`[inline]`

Definition at line 63 of file solver_mv_monomial_tests.hpp.

References mmx::eval(), mmx::rep(), and Seq< C, R >::size().

    {
      int m(val.size());
      T *ev;//[m*n];
      ev= new T [m*n];

      for( int i=0;i<m;i++)
      {
        for( int k=0;k<n;k++)
        {
          tensor::eval( ev[i*m+k] , mat[i*m+k].rep(),
                        val , m  );
        }
      }
      return ev;
    }

void mmx::realroot::scan_monom_env	(	std::set< int > &	env,
		const monom &	m
	)

Definition at line 89 of file solver_mv_bernstein.hpp.

References monom< C, E >::rep(), and mmx::size().

Referenced by solver< C, ProjRd< MTH > >::solve_monomial().

  { for ( unsigned i = 0; i < (m.rep()).size(); i++ ) if ( (m.rep())[i] != 0 ) env.insert(i); };

int mmx::realroot::sgn	(	POL &	f,
		FT	a
	)		`[inline]`

Sign of f at point a>=0 (-1=inf)

Definition at line 424 of file solver_mv_monomial_tests.hpp.

References sgn().

  {
    if (a==FT(-1) ) 
      return sgn ( f[ f.rep().szs()[0]-1 ]);
    else 
      return sgn ( f(a) );
  };

int mmx::realroot::sgn ( FT a ) [inline]

Sign of a.

Definition at line 416 of file solver_mv_monomial_tests.hpp.

Referenced by parallel< system >::process(), sgn(), descartes_solver< real_t, local_method >::solve(), and topological_degree_2d().

  {
    if (a==FT(0) ) return 0;
    return ( (a>0) ? 1 : -1 );
  }

int mmx::realroot::signof	(	polynomial< double, with< MonomialTensor > > *	p,
		Seq< FT >	t
	)		`[inline]`

Definition at line 156 of file solver_mv_monomial_tests.hpp.

References mmx::eval(), and Seq< C, R >::size().

    {
      
      Interval<double> ev;
      //int dim= p->nbvar();
      
      tensor::eval( ev , p->rep() , 
                    t , t.size() );
      
      if (ev< .1e-10) return (0);
      return (ev>0?1:-1);
    }

int mmx::realroot::signof_det	(	T *	mat,
		const int &	n
	)		`[inline]`

Definition at line 102 of file solver_mv_monomial_tests.hpp.

References det().

    { // mat= double[n*n]
      double ev;
      ev= det(mat,n);

      if (ev< .1e-15) return (0);
      return (ev>0?1:-1);
    }

T* mmx::realroot::submatrix	(	T *	matrix,
		int	order,
		int	i,
		int	j
	)

Definition at line 112 of file solver_mv_monomial_tests.hpp.

Referenced by det().

    {
      int n=order-1;
      T * subm;//[n*n];
      int p, q;         // Indexes for matrix
      int a = 0, b;     // Indexes for subm
      subm = new T [n*n];
      
      for(p = 0; p < order; p++) {
        if(p==i) continue;      //Skip ith row
        
        b = 0;
        
        for(q = 0; q < order; q++) {
          if(q==j) continue;    //Skip jth column
          subm[a*order + b++] = matrix[p*order+q];
        }
        a++; //Increment row index
      }
      return subm;
    }

int topological_degree_2d ( box_rep< POL > * b )

Definition at line 434 of file solver_mv_monomial_tests.hpp.

References assert, POL, sgn(), and Seq< C, R >::size().

  {
    /*-1=inf*/
    typedef typename POL::scalar_t C;
    
    assert( b->nbvar() ==2 );
    Seq<FT> p;
    Seq<FT> tmp;
    unsigned ev1(0), ev2(0), ev3(0), ev4(0);
    FT p1[2],p2[2];
    int tdeg;
    unsigned i;
    
      POL fl0, fr0,fl1,fr1, s,
        gl0, gr0,gl1,gr1  ;

      fl0=b->lface(0,1); gl0=b->lface(1,1);
      fr0=b->rface(0,1); gr0=b->rface(1,1);
      fl1=b->lface(0,0); gl1=b->lface(1,0);
      fr1=b->rface(0,0); gr1=b->rface(1,0);
      
      s= fl0 * gl0 ;
      p<< FT(0);p << FT(0) ;
      tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
      for (i=0;i< tmp.size(); ++i )
        p<< tmp[i] << FT(0) ;
     
      s= fr1 * gr1; 
      p<< FT(-1);p << FT(0);
      tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
      for (i=0;i< tmp.size(); ++i )
        p<< FT(-1)  << tmp[i]; 

      s=  fr0 * gr0 ;
      p<<  FT(-1);p << FT(-1) ;
      tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
      for (i=0;i< tmp.size(); ++i )
        p<<  tmp[i] << FT(-1)  ;

      s= fl1 * gl1 ; 
      p<< FT(0);p << FT(-1) ;
      tmp= solver<ring<C,MonomialTensor>, CFseparate>::template solve<FT>(s);
      for (i=0;i< tmp.size(); ++i )
        p<<  FT(0) << tmp[i]  ;
      
      p<< FT(0);p << FT(0) ;

      tdeg= 0;
      for (i=0; i<p.size()-3; i+=2)
      {
        p1[0]= p[i]  ;p1[1]=p[i+1];
        p2[0]= p[i+2];p2[1]=p[i+3];
        
        if (p1[1]==FT(0) && p2[1]==FT(0)){
          ev1= sgn( fl0, p1[0] );
          ev3= sgn( fl0, p2[0] );
          ev4= sgn( gl0, p1[0] );
          ev2= sgn( gl0, p2[0] );
        }else
          if (p1[0]==FT(-1) && p2[0]==FT(-1)){
            ev1= sgn( fr1, p1[1] );
            ev3= sgn( fr1, p2[1] );
            ev4= sgn( gr1, p1[1] );
            ev2= sgn( gr1, p2[1] );
          }else
            if (p1[1]==FT(-1) && p2[1]==FT(-1)){
              ev1= sgn( fr0, p1[0] );
              ev3= sgn( fr0, p2[0] );        
              ev4= sgn( gr0, p1[0] );
              ev2= sgn( gr0, p2[0] );
            }else
              if (p1[0]==FT(0)  && p2[0]==FT(0)){
                ev1= sgn( fl1, p1[1] );
                ev3= sgn( fl1, p2[1] );
                ev4= sgn( gl1, p1[1] );
                ev2= sgn( gl1, p2[1] );
              }
              else 
                std::cout<< "TROUBLE!"<<std::endl;
        
        tdeg+=  ( (ev1)*(ev2) - (ev3)*(ev4) );
      }
      return (tdeg/8);
    }