/Users/mourrain/Devel/mmx/realroot/include/realroot/loops_brnops.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: brnops.hpp,v 1.1 2005/07/11 10:03:56 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_LOOPS_BRNOPS_H
#define realroot_SOLVE_SBDSLV_LOOPS_BRNOPS_H
//--------------------------------------------------------------------
#include<algorithm>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------

/* boucle unidimensionelle de base pour la representation de bernstein */
namespace brnops
{
  /* algorithme de de Casteljau sur place
   * ------------------------------------
   *     r[i] = (1-t)*r[i] + t*r[i+1]
   */
  
  template<typename real_t>
  void decasteljau( real_t * r, unsigned sz, const real_t& t, int str = 1 )
  {
    real_t *er, *p;
      for ( er = r + (sz-1)*str; er != r; er -= str )
        for ( p = r; p != er; p += str )
          *p = (1.0-t)**p+t**(p+str);
  };
  
  /* algorithme de de Casteljau pour les bissections, t = 1/2 
   * --------------------------------------------------------
   *       r[i] = (r[i]+r[i+1])/2
   */
  
  template<typename real_t>
  void decasteljau( real_t * r, unsigned sz, int str = 1 )
  {
    real_t *er, *p;
    for ( er = r + (sz-1)*str; er != r; er -= str )
      for ( p = r; p != er; p += str )
        *p = (*p+*(p+str))/2.0;
  };
  
  /* évaluation en t    */
  template<typename real_t>
  real_t eval( real_t const * const src, unsigned sz, const real_t& t, int st = 1 )
  {
    real_t tmp[sz];
    std::copy( src, src + sz, tmp );
    decasteljau( tmp, sz, t );
    return tmp[0];
  };
    
  /* évaluation en (1/2) */
  template<typename real_t>
  real_t eval( real_t const * const src, unsigned sz, int st = 1 )
  {
    real_t tmp[sz];
    std::copy( src, src + sz, tmp );
    decasteljau( tmp, sz  );
    return tmp[0];
  };
  
  
  /* algorithme de de Casteljau, version subdivision
   * -----------------------------------------------
   * entrée:
   * ------
   * r = coefficients du polynome
   *
   * sortie:
   * ------
   * r = partie droite [t, 1]
   * l = partie gauche [0, t]
   */
  
  template<typename real_t>
  void decasteljau( real_t * r, real_t * l, unsigned sz, const real_t& t, int str = 1, int stl = 1 )
  {
    real_t * er, * p;
    for ( er = r + (sz-1)*str; er != r; er -= str, l += stl )
      {
        *l = *r;
        for ( p = r; p != er; p += str )
          *p = (1.0-t)**p+t**(p+str);
      };
    *l = *r;
  };
  
  /* algorithme de de Casteljau, version subdivision (1/2)
   * -----------------------------------------------
   * entrée:
   * ------
   * r = coefficients du polynome
   *
   * sortie:
   * ------
   * r = partie droite [t, 1]
   * l = partie gauche [0, t]
   */
  
  template<typename real_t>
  void decasteljau( real_t * r, real_t * l, unsigned sz, int str = 1, int stl = 1 )
  {
    real_t * er, * p;
    er = r + (sz-1)*str;
    for ( er = r + (sz-1)*str; er != r; er -= str, l += stl )
      {
        *l = *r;
        for ( p = r; p != er; p += str )
          *p = (*p+*(p+str))/2.0;
        };
    *l = *r;
    };
  
  /* restriction à gauche
   * --------------------
   * calcul de la representation pour [t,1]
   */
  
  template<typename real_t>
  void lrestrict( real_t * data, int sz, const real_t& t, int st ) 
  { decasteljau(data,sz,t,st); };
  
  /* restriction à droite
   * --------------------
   * calcul de la representation pour [0,t]
   */
  
  template<typename real_t>
  void rrestrict( real_t * data, int sz,  const real_t& t, int st ) 
  { decasteljau( data + (sz-1)*st, sz, (real_t)(real_t(1)-t), -st ); };
  
    
    
  template<typename real_t>
  void hodograph( real_t * dst, real_t const * const src, unsigned sz, int st)
  {
      int p = 0;
      for ( unsigned i = 0; i < sz-1; i ++, p += st ) dst[p] = src[p+st]-src[p];
    };
  
  template<typename real_t>
    void diff( real_t * dst, real_t const * const src, unsigned sz, int sta = 1, int stout = 1)
    {
      int pa = 0;
      int po = 0;
      for ( unsigned i = 0; i < sz-1; i ++, pa += sta, po += stout )
        dst[po] = (sz-1)*(src[pa+sta]-src[pa]);
    }; 
    
    
    template<typename real_t>
    bool rockwood_cut( real_t& t, real_t const * b, unsigned sz, int st = 1 )
    {
      for ( unsigned i = 0; i < sz-1; i ++ )
        {
          if ( b[i*st]*b[(i+1)*st] <= 0 )
            {
              real_t d = b[i*st]-b[(i+1)*st];
              t = (i*d+b[i])/(sz*d);
            return true;
            };
        };
      return false;
    };
  };
//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //