/Users/mourrain/Devel/mmx/realroot/include/realroot/solver_mv_monomial_tests.hpp

Go to the documentation of this file.

/********************************************************************
 *   Author(s): A. Mantzaflaris, GALAAD, INRIA
 ********************************************************************/

# pragma once
/********************************************************************/

# include <realroot/ring.hpp>
# include <realroot/polynomial_tensor.hpp>

# include <realroot/polynomial.hpp>
# include <realroot/homography.hpp>
# include <realroot/sign_variation.hpp>
# include <realroot/Interval.hpp>
//# include <linalg/matrix.hpp>
# include <stack>

# include <realroot/solver_ucf.hpp>
# include <realroot/solver_uv_bspline.hpp>


# define TMPL  template<class POL>
# define DPOL polynomial<double,with<MonomialTensor> >
# define BOX box_rep<POL>

# define INFO 1
# define TODOUBLE as<double>

//====================================================================
namespace mmx {
//====================================================================
namespace realroot {

  template<class POL, class FT>  int topological_degree_2d(BOX* b);


//Jacobian 

  TMPL
  inline DPOL * jacobian(Seq<POL*>& S0)
    {
      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;
    }


  template<class T> inline 
  void precondition (Seq<T> val, DPOL * mat, int n)
    {
      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;
    }




  template<class T> inline 
  T * eval_poly_matrix(Seq<T> val, DPOL * mat, int n)
    {
      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;
    }

  template<class T> inline 
  int signof_det(T *mat,const int& n)
    { // mat= double[n*n]
      double ev;
      ev= det(mat,n);

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

  template<class T>
  T* submatrix(T *matrix, int order, int i, int j)
    {
      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;
    }

  // Compute matrix determinant, recursively
  template<class T>  
  T det(T *mat, int order)
    {
      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;
    }


//
  template<class FT> 
  inline int signof(DPOL * p, Seq<FT> t)
    {
      
      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);
    }
  
    
  /*
    template<class C> 
    inline bool is_root(Seq<C> t)
        {
            assert( t.size()==dim );
            C ev;
            for ( unsigned i=0; i!=S.size() ; ++i )
            {
                tensor::eval( ev , S[i].rep() , t , dim );
                if ( ev != 0 ) 
                    return false;
            }
            return true;
        }
  */    

  TMPL
  bool inline miranda_test(Seq<POL>& S,const int i,const int j)
  {
    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) );
  };
  
  // assumes a square system for now
  template<class POL>
    bool include1(BOX* b, DPOL * J) 
    {
      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;
    };
  
  TMPL
  bool include2(BOX* b,DPOL * J) 
  {
    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) );
  }

  template<class POL>
  bool include3(BOX* b,DPOL * J, Seq<POL*>& S0) 
        {
          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 );
        }


  template<class POL>
  bool exclude1(BOX*b, Seq<POL *>& S0) 
    {
      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;
    };


  TMPL
  bool exclude2(BOX* b,DPOL * J) 
  {
    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 );
  }

    TMPL
    inline bool exclude3(BOX* b) 
    {
      Seq<POL> S= b->system();
      for (unsigned i=0; i!=b->nbpol(); ++i)
        if ( no_variation( S[i]) ) 
          return true;
      return false;
    }
  
  TMPL
  inline bool no_variation(POL & p)
  {
    return !has_sign_variation(p.begin(),p.end()) ;
  };
  

//====================================================================
// ---------------------- Topological degree in 2d -------------------
//====================================================================

  template<class FT>
  inline int sgn(FT a )
  {
    if (a==FT(0) ) return 0;
    return ( (a>0) ? 1 : -1 );
  }
  
  template<class POL, class FT>   
  int inline sgn(POL & f, FT  a )
  {
    if (a==FT(-1) ) 
      return sgn ( f[ f.rep().szs()[0]-1 ]);
    else 
      return sgn ( f(a) );
  };
  
  //template<class POL, class FT>
  template<class POL, class FT>
  int topological_degree_2d(BOX* b)
  {
    /*-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);
    }
  

} //namespace realroot
} //namespace mmx 

# undef DPOL