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 //