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

Go to the documentation of this file.

/********************************************************************
 *   This file is part of the source code of the realroot library.
 *   Author(s): J.P. Pavone, GALAAD, INRIA
 *   $Id: descartes1d.hpp,v 1.1 2005/07/11 10:03:57 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_USOLVER_DESCARTES1D_HPP
#define realroot_SOLVE_SBDSLV_USOLVER_DESCARTES1D_HPP
//--------------------------------------------------------------------
#include <realroot/bernstein_bzenv.hpp>
#include <realroot/loops_vctops.hpp>
#include <realroot/loops_upoldse.hpp>
#include <realroot/loops_brnops.hpp>
//#include <arithmetix/TypeTests.h>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
template< typename Prms >
struct sign_wanted { typedef texp::false_t result_t; };

namespace solvers
{
  /* méthode locales */
  
  template<class real_t> 
  struct bsearch         /* bissection sous forme monomiale */
  {
    real_t * m_data;
    real_t * m_mons;
    unsigned m_sz;
    unsigned m_mxsz;
    
    template< class In > 
    /* initialisation: conversion sous forme monomiale */
    bsearch( const In&  bzrep,  unsigned sz ) : m_sz(sz) 
    {
      m_data = new real_t[ 2*sz ];
      std::copy( bzrep, bzrep + sz, m_data );
      m_mons = m_data + sz;
      std::copy( m_data, m_data + sz, m_mons );
      bernstein::bzenv<real_t>::_default_->toMonoms( m_mons, sz );
    };
    ~bsearch() { delete[] m_data; };
    /* bissection */
    void reach( real_t * lbzrep, real_t& a, real_t& b, const real_t& eps  ) 
    {
      real_t m;
      if ( lbzrep[m_sz-1] > lbzrep[0] ) 
        do
          if ( upoldse_::horner( m_mons, m_sz, m=(a+b)/2.0 ) < 0  ) a = m; else b = m; 
        while( b-a > eps ) ;
      else
        do
          if ( upoldse_::horner( m_mons, m_sz, m=(a+b)/2.0 ) < 0  ) b = m; else a = m;
        while( b-a > eps ); 
    };
  };

  template<class real_t> 
  struct bsearch_castel     /* bissection sous forme de Bernstein */
  {
    real_t * m_data;
    unsigned m_sz;
    template< class In >
    /* initialisation: copie de la forme de Bernstein */
    bsearch_castel( const In&  bzrep,  unsigned sz ) : m_sz(sz) 
    {
      m_data = new real_t[ sz ];
      std::copy( bzrep, bzrep + sz, m_data );
    };
    ~bsearch_castel() { delete[] m_data; };
    
    /* bissection */
    void reach( real_t * lbzrep, real_t& a, real_t& b, const real_t& eps  ) 
    {
      real_t m;
      if ( lbzrep[m_sz-1] > lbzrep[0] ) 
        do
          if ( brnops::eval( m_data, m_sz, m = (a+b)/2.0 ) < 0  ) a = m; else b = m; 
        while( b-a > eps ) ;
      else
        do
          if ( brnops::eval( m_data, m_sz, m = (a+b)/2.0 ) < 0  ) b = m; else a = m;
        while( b-a > eps ); 
    };
  };
  
  template<class real_t> 
  struct bsearch_newton       /* bissection / newton mix  */
  {
    real_t * m_data;
    real_t * m_mons;
    unsigned m_sz;
    template< class In >
    /* initialisation conversion sous forme monomiale */
    bsearch_newton( const In&  bzrep,  unsigned sz ) : m_sz(sz) 
    {
      m_data = new real_t[ 2*sz ];
      std::copy( bzrep, bzrep + sz, m_data );
      m_mons = m_data + sz;
      std::copy( m_data, m_data + sz, m_mons );
      bernstein::bzenv<real_t>::_default_->toMonoms( m_mons, sz );
    };
    ~bsearch_newton() { delete[] m_data; };
    
    /* bissection / newton mix */
    void reach( real_t * lbzrep, real_t& a, real_t& b, const real_t& eps  ) 
    {
      real_t m;
      if ( lbzrep[m_sz-1] > lbzrep[0] ) 
        do
          {
            real_t p,dp,x;
            /* évaluation du polynôme en m (p) et de sa derivée (dp) */
            upoldse_::dhorner( p, dp, m_mons, m_sz, m = (a+b)/2.0 );
            /* étape de bissection                                   */
            if ( p < 0  ) a = m; else b = m; 
            if ( dp > eps ) /* si la valeur de la derivée n'est pas trop faible */
              {
                /* on fait une itération de la méthode de newton */
                x = m - p/dp;
                /* si le résultat est dans l'intervalle de la bissection */
                if ( x > a && x < b )
                  {
                    /* évaluation du polynôme en ce point */
                    real_t xp = upoldse_::horner( m_mons, m_sz, x ); 
                    /* correction d'une des deux bornes */
                    if ( xp < 0 ) a = x; else b = x;
                    /* recherche de la seconde          */
                    if ( xp < 0 )
                      {
                        real_t step = eps;
                        while( xp < 0 && x < b )
                          {
                            x += step;
                            xp = upoldse_::horner( m_mons, m_sz, x );
                            step *= 2;
                          };
                        if ( x < b ) b = x;
                      }
                    else
                      {
                        real_t step = eps;
                        while( xp > 0 && x > a )
                          {
                            x -= step;
                            xp = upoldse_::horner( m_mons, m_sz, x );
                            step *= 2;
                          };
                        if ( x > a ) a = x;
                      };
                  };
              };
          }
        while( b-a > eps );
      else
        do
          {
            real_t p,dp,x;
            upoldse_::dhorner( p, dp, m_mons, m_sz, m = (a+b)/2.0 );
            if ( p < 0  ) b = m; else a = m; 
            if ( dp > eps )
              {
                x = m - p/dp;        
                if ( x > a && x < b )
                  {
                    real_t xp = upoldse_::horner( m_mons, m_sz, x );
                    if ( xp < 0 ) b = x; else a = x;
                    if ( xp < 0 )
                      {
                        real_t step = eps/2;
                        while( xp < 0 && x > a )
                          {
                            x -= step;
                            xp = upoldse_::horner( m_mons, m_sz, x );
                            step *= 2;
                          };
                        if ( x > a  ) a = x;
                      }
                    else
                      {
                        real_t step = eps/2;
                        while( xp > 0 && x < b )
                          {
                            x += step;
                            xp = upoldse_::horner( m_mons, m_sz, x );
                            step *= 2;
                          };
                        if ( x < b ) b = x;
                      };
                  };
              };
          }
        while( b-a > eps ); 
    };
  };

  template<class real_t> 
  struct bsearch_newton2
  {
    real_t * m_data;
    real_t * m_mons;
    unsigned m_sz;
    template< class In >
    bsearch_newton2( const In&  bzrep,  unsigned sz ) : m_sz(sz) 
    {
      m_data = new real_t[ 2*sz ];
      std::copy( bzrep, bzrep + sz, m_data );
      m_mons = m_data + sz;
      std::copy( m_data, m_data + sz, m_mons );
      bernstein::bzenv<real_t>::_default_->toMonoms( m_mons, sz );
    };
    ~bsearch_newton2() { delete[] m_data; };
    
    void reach( real_t * lbzrep, real_t& a, real_t& b, const real_t& eps  ) 
    {
      real_t m;
      if ( lbzrep[m_sz-1] > lbzrep[0] ) 
        do
          {
            //    -
            //     - 
            // -----------------
            //       -
            //        -
            //std::cout << (b-a) << std::endl;
            real_t p,dp,x;
            upoldse_::dhorner( p, dp, m_mons, m_sz, m = (a+b)/2.0 );
            if ( p < 0  ) a = m; else b = m; 
            if ( dp > eps )
              {
                real_t ex;
                int c = 0;
                x = m;
                do
                  {
                    ex = x;
                    x = x - p/dp;
                    c ++ ;
                    upoldse_::dhorner(p,dp,m_mons,m_sz,x);
                  }
                while ( c < 6 );
                if ( x > a && x < b ) 
                  {
                    if ( p < 0 ) 
                      {
                        a = x;
                        p = upoldse_::horner(m_mons,m_sz,x+eps/2);
                        if ( p > 0 ) b = x+eps;
                      }
                    else 
                      {
                        b = x;
                        p = upoldse_::horner(m_mons,m_sz,x-eps/2);
                        if ( p < 0 ) a = x-eps/2;
                      };
                  };
              };
          }
        while( b-a > eps );
      else
        do
          {
            //              std::cout << (b-a) << std::endl;
            real_t p,dp,x;
            upoldse_::dhorner( p, dp, m_mons, m_sz, m = (a+b)/2.0 );
            if ( p < 0  ) b = m; else a = m; 
            if  ( dp > eps )
              {
                real_t ex;
                int c = 0;
                x = m;
                do
                  {
                    ex = x;
                    x = x - p/dp;
                    c ++ ;
                    upoldse_::dhorner(p,dp,m_mons,m_sz,x);
                  }
                while (c < 6 );
                
                if ( x > a && x < b ) 
                  {
                    if ( p < 0 ) 
                      {
                        b = x; 
                        p = upoldse_::horner(m_mons,m_sz,x-eps/2);
                        if ( p > 0 ) a = x-eps/2;
                      }
                    else 
                      {
                        a = x;
                        p = upoldse_::horner(m_mons,m_sz,x+eps/2);
                        if ( p < 0 ) b = x-eps/2;
                      };
                  };
              }
          }
        while( b-a > eps ); 
    };
  };

  template<class real_t, class local_method = bsearch< real_t > >
  class descartes_solver
  {
  public:
    unsigned m_sz, m_ssz;
    real_t        m_prec;
    real_t *      m_stck;
    local_method * m_lmth;

      
    void reset( unsigned sz, const real_t& prec )
    {
      if ( sz <= m_sz && prec >= m_prec ) 
        {
          m_sz = sz;
          m_prec = prec;
          return;
        };
      m_sz = sz; m_prec = prec; 
      this->alloc();
    };
    
    void split( real_t * r )
    {
      real_t * l = r+m_sz+2;
      brnops::decasteljau(r,l,m_sz);
      l[m_sz]   = r[m_sz];
      l[m_sz+1] = (r[m_sz]+r[m_sz+1])/2.0;
      r[m_sz]   = l[m_sz+1];
    };
    inline real_t precision( real_t * l ) { return l[m_sz+1]-l[m_sz] ; };
    bool precision_reached( real_t * l ) { return l[m_sz+1]-l[m_sz] < m_prec; };
    
  public:  
    /* main loop */
    /* 
       input : in  is a one dimensional vector containing bernstein coefficients (implement an operator[])
       output: out is a function called on each interval containing roots        (assumed to store the results)
       option: 0 to compute all the roots, n to compute the nth root
    */
    
    

    void rockwoodcuts( real_t * nxt, real_t * prv, real_t * mid )
    {
      //      [-------------] b
      //      [---][--------] b 
      //      [---][---][---] b
      real_t i0, i1;
      brnops::rockwood_cut( i0, nxt, m_sz );
      brnops::rockwood_cut( i1, nxt+m_sz-1, m_sz, -1 );
      std::cout << i0 << ", " << (1.0-i1) << std::endl;
      i1 = (1.0-i1);
      brnops::decasteljau(nxt,prv,m_sz,i0);
      brnops::decasteljau(nxt,mid,m_sz,(i1-i0)/(1.0-i0));
      i0 = (1.0-i0)*nxt[m_sz] + i0*nxt[m_sz+1];
      i1 = (1.0-i1)*nxt[m_sz] + i1*nxt[m_sz+1];
      prv[m_sz]   = nxt[m_sz];
      prv[m_sz+1] = i0;
      mid[m_sz]   = i0;
      mid[m_sz+1] = i1;
      nxt[m_sz]   = i1;
    };
    
    inline void set_sz( unsigned sz ) { m_sz = sz;};
    
    template< class Prms > inline 
    void output( Prms& prms, real_t * stck, const texp::true_t&  )
    { prms.output( stck[m_sz], stck[m_sz+1], stck[0] > 0, stck[m_sz-1] > 0 ); };
    template< class Prms > inline 
    void output( Prms& prms, real_t * stck, const texp::false_t& )
    { prms.output( stck[m_sz], stck[m_sz+1] ); };
    
    template< class Prms > inline
    void output( Prms& prms, real_t * stck )
    { output( prms, stck, sign_wanted< Prms >::result_t() ); };
     
    template<class Prms,class In>
    unsigned solve( Prms& prms, const In& in, int option = 0  )
    {
      for ( unsigned i = 0; i < m_sz; i ++ )  m_stck[i] = in[i]; 
      m_stck[m_sz] = 0.0; m_stck[m_sz+1] = 1.0;
      unsigned    roots = 0;
      const unsigned    inc  = m_sz+2;
      real_t * stck = m_stck;
      real_t * stop = stck-inc;
      do
        {
          unsigned sgn = vctops::sgncnt(stck,m_sz);
          switch( sgn )
            {
            case 0:
              stck -= inc;
              break;
            case 1:
              roots ++;
              if ( !option || (roots == option) )
                {
                  if ( m_lmth == 0 ) m_lmth = new local_method(in,m_sz);
                  m_lmth->reach( stck, stck[m_sz], stck[m_sz+1], m_prec );
                  output( prms, stck );
                };
              stck -= inc;
              break;
            default:
              if(  precision( stck ) < m_prec  )
                { 
                  roots++; 
                  if ( !option || (roots == option) ) output(prms,stck);
                  stck -= inc;
                  break;
                };
              split(stck);
              stck += inc;
              break;
            }
        }
      while( stck != stop  );
      if ( m_lmth ) { delete m_lmth; m_lmth = 0; };
      return roots;
    };

        

    descartes_solver( unsigned sz, const real_t& eps = 1e-12 ) 
    {
      m_stck = 0;
      m_sz   = sz;
      m_prec = eps;
      m_lmth = 0;
      unsigned maxdeep = 1;
      real_t ex = 1;
      while ( ex > m_prec ) { ex/=2; maxdeep ++ ; };
      //      std::cout << maxdeep << std::endl;
      maxdeep ++;
      maxdeep *= 3;
      unsigned tsz = maxdeep*(2+m_sz);
      m_stck = new real_t[ tsz ];
    };
    
    ~descartes_solver() { 
      delete[] m_stck; 
    };
  };
};
//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //