#include <system_descartes1d.hpp>

Public Member Functions

template<class In > bsearch_newton (const In &bzrep, unsigned sz)
~bsearch_newton ()
void reach (real_t *lbzrep, real_t &a, real_t &b, const real_t &eps)

Public Attributes

real_t * m_data
real_t * m_mons
unsigned m_sz

Detailed Description

template<class real_t>
struct mmx::solvers::bsearch_newton< real_t >

Definition at line 88 of file system_descartes1d.hpp.

Constructor & Destructor Documentation

bsearch_newton	(	const In &	bzrep,
		unsigned	sz
	)		`[inline]`

Definition at line 95 of file system_descartes1d.hpp.

References mmx::sparse::copy(), bsearch_newton< real_t >::m_data, bsearch_newton< real_t >::m_mons, and bzenv< X >::toMonoms().

                                                     : 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 ( ) [inline]

Definition at line 103 of file system_descartes1d.hpp.

References bsearch_newton< real_t >::m_data.

{ delete[] m_data; };

Member Function Documentation

void reach	(	real_t *	lbzrep,
		real_t &	a,
		real_t &	b,
		const real_t &	eps
	)		`[inline]`

Definition at line 106 of file system_descartes1d.hpp.

References mmx::upoldse_::dhorner(), mmx::upoldse_::horner(), bsearch_newton< real_t >::m_mons, and bsearch_newton< real_t >::m_sz.

    {
      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 ); 
    };