polynomial_algo.hpp

Go to the documentation of this file.

/*
 * Copyright(c):
 * Signal Image and Communications (SIC) Department 
 * http://www.sic.sp2mi.univ-poitiers.fr/ 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - XLIM  Institute UMR CNRS 7252  http://www.xlim.fr/ 
 * 
 * and 
 * 
 * D2 Fluid, Thermic and Combustion 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - PPRIME Institute -  UPR CNRS 3346  http://www.pprime.fr
 *  - ISAE-ENSMA http://www.ensma.fr
 *
 * Contributor(s):
 * The SLIP team,
 * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
 * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>, 
 * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>, 
 * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>, 
 * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>, 
 * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
 *
 * Description:
 * The Simple Library of Image Processing (SLIP) is a new image processing 
 * library. It is written in the C++ language following as much as possible 
 * the ISO/ANSI C++ standard. It is consequently compatible with any system  
 * satisfying the ANSI C++ complience. It works on different Unix , Linux ,  
 * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that  
 * was created by the Signal, Image and Communications (SIC) departement of  
 * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic  
 * and Combustion departement of the P', UPR 3346 CNRS Institute of the  
 * University of Poitiers.
 * 
 * The SLIP Library source code has been registered to the APP (French Agency 
 * for the Protection of Programs) by the University of Poitiers and CNRS, 
 * under  registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.

 * http://www.sic.sp2mi.univ-poitiers.fr/slip/
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use, 
 * modify and/ or redistribute the software under the terms of the CeCILL-C 
 * license as circulated by CEA, CNRS and INRIA at the following URL 
 * http://www.cecill.info.
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 * 
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the 
 * same conditions as regards security. 
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 */



#ifndef SLIP_POLYNOMIAL_ALGO_HPP
#define SLIP_POLYNOMIAL_ALGO_HPP

#include "statistics.hpp"
#include "Array2d.hpp"
#include "arithmetic_op.hpp"
//#include "MultivariatePolynomial.hpp"
#include <iterator>
#include <iostream>
#include <algorithm>
#include <numeric>
#include <cmath>


namespace slip
{

 
  /* @{ */
  template <typename Poly, typename Poly1, typename Poly2, typename Poly3>
  inline Poly
  legendre_nd_next(std::size_t k, Poly1 x, Poly2 Pk, Poly3 Pkm1)
  {
    typedef typename Poly::mapped_type T;
   
    return ((((Poly(x) * Poly(Pk))* T(2 * k + 1)) - (Poly(Pkm1)* T(k))) / T(k + 1));
    
  }

  template <typename Poly, typename Poly1, typename Poly2, typename Poly3>
  inline Poly
  legendre_next(std::size_t k, Poly1 x, Poly2 Pk, Poly3 Pkm1)
  {
    typedef typename Poly::value_type T;
   
    return ((((Poly(x) * Poly(Pk))* T(2 * k + 1)) - (Poly(Pkm1)* T(k))) / T(k + 1));
    
  }

  /* @} */


  /* @{ */

 template <typename Poly, 
                   typename Poly1, 
                   typename Poly2, 
                   typename Poly3>
 inline Poly
 chebyshev_next(std::size_t k, Poly1 x, Poly2 Pk, Poly3 Pkm1)
  {
    typedef typename Poly::value_type T;
   
    return ((Poly(x) * Poly(Pk))* T(2)) - Poly(Pkm1);
    
  }

  template <typename Poly, typename Poly1, typename Poly2, typename Poly3>
  inline Poly
  chebyshev_nd_next(std::size_t k, Poly1 x, Poly2 Pk, Poly3 Pkm1)
  {
    typedef typename Poly::mapped_type T;
   
    return ((Poly(x) * Poly(Pk))* T(2)) - Poly(Pkm1);
    
  }

  template<typename T, typename RandomAccessIterator>
  inline
  void eval_chebyshevII_basis(const T& x, 
                                                   const std::size_t n,
                                                   RandomAccessIterator first, RandomAccessIterator last)
  {
    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(n + 1));
    if(n == 0)
      {
                first[0] = T(1);
      }
    else if (n == 1)
      {
                 first[0] = T(1);
                 first[1] = 2*x;
      }
    else
      {

                T twox = T(2) * x;
                first[0] = T(1);
                first[1] = twox;

                
                for (std::size_t j = 2; j <= n; ++j)
                  {
                    first[j] = (twox * first[j-1]) - first[j-2];
                  }
                 
      }
   
  }

  
  template<typename T, typename RandomAccessIterator>
  inline
  void eval_discrete_chebyshev_basis(const T& x, 
                                                                     const std::size_t n,
                                                                     RandomAccessIterator first, RandomAccessIterator last)
  {
    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(n + 1));
    T N =  static_cast<T>(last-first);
    if(n == 0)
      {
                first[0] = T(1);
      }
    else if (n == 1)
      {
                 first[0] = T(1);
                 first[1] = T(2)*x - N + T(1);
      }
    else
      {

                T twox = T(2)*x - N + T(1);
                first[0] = T(1);
                first[1] = twox;
                T N2 = N * N;
                
                for (std::size_t j = 2; j <= n; ++j)
                  {
                    T jf = static_cast<T>(j);
                    T jf_m_1 = jf - static_cast<T>(1);
                    T alpha = (jf_m_1 + jf_m_1 + static_cast<T>(1)) / jf;
                    T jf_m_12 = static_cast<T>(jf_m_1 * jf_m_1);
                    T beta = (static_cast<T>(jf_m_1)*(N2-jf_m_12))/jf;
                    
                    first[j] = alpha*(twox * first[j-1]) - beta*first[j-2];
                  }
                 
      }
   
  }

 /* @} */  

  /* @{ */
 template <typename Poly, 
                   typename Poly1, 
                   typename Poly2, 
                   typename Poly3>
 inline Poly
 hermite_next(std::size_t k, Poly1 x, Poly2 Pk, Poly3 Pkm1)
  {
    typedef typename Poly::value_type T;
   
    return ((Poly(x) * Poly(Pk))* T(2)) - Poly(Pkm1)* T(2*k);
    
  }

  template <typename Poly, typename Poly1, typename Poly2, typename Poly3>
  inline Poly
  hermite_nd_next(std::size_t k, Poly1 x, Poly2 Pk, Poly3 Pkm1)
  {
    typedef typename Poly::mapped_type T;
   
    return ((Poly(x) * Poly(Pk))* T(2)) - Poly(Pkm1)* T(2*k);
    
  }

  template<typename T, typename RandomAccessIterator>
  inline
  void eval_hermite_basis(const T& x, 
                                                  const std::size_t n,
                                                  RandomAccessIterator first, RandomAccessIterator last)
  {
    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(n + 1));
    if(n == 0)
      {
                first[0] = T(1);
      }
    else if (n == 1)
      {
                 first[0] = T(1);
                 first[1] = 2*x;
      }
    else
      {

                T twox = T(2) * x;
                first[0] = T(1);
                first[1] = twox;

                
                for (std::size_t j = 2; j <= n; ++j)
                  {
                    first[j] = (twox * first[j-1]) - T(2*j-2)*first[j-2];
                  }
                 
      }
   
  }
/* @} */  

  /* @{ */
 template <typename Poly, 
                   typename Poly1, 
                   typename Poly2, 
                   typename Poly3>
 inline Poly
 laguerre_next(std::size_t k, Poly1 x, Poly2 Pk, Poly3 Pkm1)
  {
    typedef typename Poly::value_type T;
   
    return (((T(2*k+1) - Poly(x)) * Poly(Pk)) - Poly(Pkm1)* T(k))/T(k+1);
    
  }



template <typename Poly, typename Poly1, typename Poly2, typename Poly3>
inline Poly
laguerre_nd_next(std::size_t k, Poly1 x, Poly2 Pk, Poly3 Pkm1)
  {
    typedef typename Poly::mapped_type T;
   
    return (((T(2*k+1) - Poly(x)) * Poly(Pk)) - Poly(Pkm1)* T(k))/T(k+1);
    
  }


  template<typename T, typename RandomAccessIterator>
  inline
  void eval_laguerre_basis(const T& x, 
                                                   const std::size_t n,
                                                   RandomAccessIterator first, RandomAccessIterator last)
  {
    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(n + 1));
    if(n == 0)
      {
                first[0] = T(1);
      }
    else if (n == 1)
      {
                 first[0] = T(1);
                 first[1] = -x + T(1);
      }
    else
      {

                first[0] = T(1);
                first[1] = -x + T(1);

                
                for (std::size_t j = 2; j <= n; ++j)
                  {
                    first[j] = -(((x - T(2*j-1))* first[j-1]) + T(j-1)*first[j-2])/T(j);
                  }
                 
      }
   
  }
  /* @} */ 

  /* @{ */
  template<typename T, typename RandomAccessIterator>
  inline
  T eval_horner(RandomAccessIterator first, RandomAccessIterator last, 
                                const T& x)
  {
    assert(last != first);
    T s = T(0);
    if(x == T(0))
      {
                s = *first;
      }
    else if(x == T(1))
      {
                s = std::accumulate(first,last,T(0));
      }
    else if(x == T(-1))
      {
                typedef typename std::iterator_traits<RandomAccessIterator>::difference_type _Distance;
                _Distance size = (last - first);
                s = *first;
                for(_Distance i = 1; i < size; ++i)
                  {
                    if(i%2 == 0)
                      {
                                s = s + first[i];
                      }
                    else
                      {
                                s = s - first[i];
                      }
                  }
      }
    else
      {
                typedef typename std::iterator_traits<RandomAccessIterator>::difference_type _Distance;
                
                _Distance size = last - first;
                s = first[size - 1];
                for(_Distance i = (size - 2); i >= 0; --i)
                  {
                    s = first[i] + x * s;
                  }
      }
    return s;
  }  

  template<typename T, typename RandomAccessIterator>
  inline
  void eval_power_basis(const T& x, 
                                                const std::size_t n,
                                                RandomAccessIterator first, RandomAccessIterator last)
  {
    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(n + 1));
    if(n == 0)
      {
                *first = T(1);
      }
    else if (n == 1)
      {
                 *first++ = T(1);
                 *first++ = x;
      }
    else
      {
                *first++ = T(1);

                if(x == T(0))
                  {
                    std::fill(first,last,T(0));
                  }
                else if (x == T(1))
                  {
                    std::fill(first,last,T(1));
                  }
                else if (x == T(-1))
                  {
                    
                    *first = T(-1);
                    std::size_t i = 0;
                    while(++first != last)
                      {
                                if(i%2 == 0)
                                  {
                                    *first = T(1);
                                  }
                                else
                                  {
                                    *first = T(-1);
                                  }
                                i++;
                      }
                  }
                else
                  {
                    T pow = x;
                    *first = pow;
                    
                    while(++first != last)
                      {
                                pow = pow * x;
                                *first = pow;
                      }
                  }
      }
   
  }
  
  template<typename T, typename RandomAccessIterator>
  inline
  void eval_legendre_basis(const T& x, 
                                                   const std::size_t n,
                                                   RandomAccessIterator first, RandomAccessIterator last)
  {
    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(n + 1));
    if(n == 0)
      {
                first[0] = T(1);
      }
    else if (n == 1)
      {
                 first[0] = T(1);
                 first[1] = x;
      }
    else
      {
                first[0] = T(1);
                first[1] = x;

                T twox = T(2) * x;
                T f2 = x;
                T d = T(1);
                T f1 = T(0);

                for (std::size_t j = 2; j <= n; ++j)
                  {
                    f1 = d++;
                    f2 += twox;
                    first[j] = (f2 * first[j-1] - f1 * first[j-2]) / d;
                  }
                
      }
   
  }

  template<typename T, typename RandomAccessIterator>
  inline
  void eval_chebyshev_basis(const T& x, 
                                                   const std::size_t n,
                                                   RandomAccessIterator first, RandomAccessIterator last)
  {
    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(n + 1));
    if(n == 0)
      {
                first[0] = T(1);
      }
    else if (n == 1)
      {
                 first[0] = T(1);
                 first[1] = x;
      }
    else
      {
                first[0] = T(1);
                first[1] = x;

                T twox = T(2) * x;
                for (std::size_t j = 2; j <= n; ++j)
                  {
                    first[j] = (twox * first[j-1]) - first[j-2];
                  }
                
      }
   
  }


  template<class Vector, typename RandomAccessIterator>
  inline
  void eval_power_2d_basis(const Vector& x, 
                                                   const std::size_t n,
                                                   RandomAccessIterator first, RandomAccessIterator last)
  {
    assert(x.size() == 2);
    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(((n + 1)*(n+2))/2)); 
    slip::Array2d<typename std::iterator_traits<RandomAccessIterator>::value_type> Eval(x.size(),n+1);
   
    for(std::size_t i = 0; i < Eval.rows(); ++i)
      {
                slip::eval_power_basis(x[i],n,Eval.row_begin(i),Eval.row_end(i));
      }

    for(std::size_t i = 0; i < (n+1); ++i)
      {
                for(std::size_t j = 0; j < ((n+1) - i); ++j)
                  {
                    *first++ = Eval[0][j] * Eval[1][i];
                  }
      }
  }


 
  template<class Vector, typename RandomAccessIterator>
  inline
  void eval_power_nd_basis(const Vector& x, 
                                                   const std::size_t order,
                                                   RandomAccessIterator first, RandomAccessIterator last)
  {
    std::size_t dim = x.size(); 
    slip::Array<std::size_t> V(dim+order);
    slip::iota(V.begin(),V.end(),1,1);
    std::size_t result_size = (std::accumulate(V.begin()+order,V.end(),1,std::multiplies<int>()))/(std::accumulate(V.begin(),V.begin()+dim,1,std::multiplies<int>()));

    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(result_size));
    
    
    //array of all the powers combinations sort by lexicographical order
     slip::Array2d<std::size_t> powers(slip::power(order+1,dim),dim);
    for(std::size_t d = 0; d < powers.dim2(); ++d)
      {
                std::size_t step = slip::power((order+1),d);
                for(std::size_t ord = 0, ord2 = 0; ord <= (powers.dim1() - step); ord+=step,++ord2)
                  {
                    std::size_t tmp = ord2%(order+1);
                    for(std::size_t k = 0; k < step; ++k)
                      {
                                powers[ord+k][d] = tmp;
                      }
                    
                  }
      }
   
    typedef typename std::iterator_traits<RandomAccessIterator>::value_type value_type;
    //Pre-computation of the powers of x
    slip::Array2d<value_type> Eval(x.size(),order+1);
    for(std::size_t i = 0; i < Eval.rows(); ++i)
      {
                slip::eval_power_basis(x[i],order,Eval.row_begin(i),Eval.row_end(i));
      }
   
    //Evaluation of the monomials
    for(std::size_t ord = 0; ord < powers.dim1(); ++ord)
      {
                //              if(std::accumulate(powers.row_begin(ord),powers.row_end(ord),value_type()) <= order)
                if(std::accumulate(powers.row_begin(ord),powers.row_end(ord),static_cast<std::size_t>(0)) <= order)
                  {
                   
                    value_type res = value_type(1);
                    for(std::size_t d = 0; d < dim; ++d)
                      {
                                res = res * Eval[d][powers[ord][d]];
                      }
                    *first++ = res;
                  }
      }


  }

  template<typename MultivariatePolynomial, 
                   typename RandomAccessIterator>
  inline
void eval_multi_poly_power_basis(const MultivariatePolynomial& x, 
                                                                 const std::size_t n,
                                                                 RandomAccessIterator first, 
                                                                 RandomAccessIterator last)
  {
    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(n + 1));
    if(n == 0)
      {
                *first = MultivariatePolynomial(0);
      }
    else if (n == 1)
      {
                 *first++ = MultivariatePolynomial(0);
                 *first++ = x;
      }
    else
      {
                *first++ = MultivariatePolynomial(0);

                MultivariatePolynomial pow = x;
                *first = pow;
                
                while(++first != last)
                  {
                    pow = pow * x;
                    *first = pow;
                  }
                
      }
   
  }

template<class Vector, typename RandomAccessIterator>
inline
void eval_multi_poly_power_nd_basis(const Vector& x, 
                                                                    const std::size_t order,
                                                                    RandomAccessIterator first, RandomAccessIterator last)
  {
    std::size_t dim = x.size(); 
    slip::Array<std::size_t> V(dim+order);
    slip::iota(V.begin(),V.end(),1,1);
    std::size_t result_size = (std::accumulate(V.begin()+order,V.end(),1,std::multiplies<int>()))/(std::accumulate(V.begin(),V.begin()+dim,1,std::multiplies<int>()));

    assert((last - first) == typename std::iterator_traits<RandomAccessIterator>::difference_type(result_size));
    
    
    //array of all the powers combinations sort by lexicographical order
     slip::Array2d<std::size_t> powers(slip::power(order+1,dim),dim);
    for(std::size_t d = 0; d < powers.dim2(); ++d)
      {
                std::size_t step = slip::power((order+1),d);
                for(std::size_t ord = 0, ord2 = 0; ord <= (powers.dim1() - step); ord+=step,++ord2)
                  {
                    std::size_t tmp = ord2%(order+1);
                    for(std::size_t k = 0; k < step; ++k)
                      {
                                powers[ord+k][d] = tmp;
                      }
                    
                  }
      }
   
    typedef typename std::iterator_traits<RandomAccessIterator>::value_type value_type;
    //Pre-computation of the powers of x
    slip::Array2d<value_type> Eval(x.size(),order+1);
    for(std::size_t i = 0; i < Eval.rows(); ++i)
      {
                slip::eval_multi_poly_power_basis(x[i],order,Eval.row_begin(i),Eval.row_end(i));
      }
   
    //Evaluation of the monomials
    for(std::size_t ord = 0; ord < powers.dim1(); ++ord)
      {
                if(std::accumulate(powers.row_begin(ord),powers.row_end(ord),static_cast<std::size_t>(0)) <= order)
                  {
                   
                    value_type res = value_type(0);
                    for(std::size_t d = 0; d < dim; ++d)
                      {
                                res = res * Eval[d][powers[ord][d]];
                      }
                    *first++ = res;
                  }
      }


  }

/* @} */
  
  template<typename T,typename RandomAccessIterator>
  struct EvalBasis 
  {
   virtual void operator()(const T& x, 
                    const std::size_t n,
                    RandomAccessIterator first, 
            RandomAccessIterator last ) const = 0;
    virtual ~EvalBasis(){} 
    
   };

 template<typename T,typename RandomAccessIterator>
 struct EvalPowerBasis: public EvalBasis<T,RandomAccessIterator>
 {
   void operator()(const T& x, 
                     const std::size_t n,
                     RandomAccessIterator first, 
                     RandomAccessIterator last ) const
    {
      slip::eval_power_basis(x,n,first,last);
    }
  
 };
   
  template<typename T,typename RandomAccessIterator>
  struct EvalLegendreBasis: public EvalBasis<T,RandomAccessIterator>
  {
    void operator()(const T& x, 
                     const std::size_t n,
                     RandomAccessIterator first, 
                     RandomAccessIterator last ) const
    {
      slip::eval_legendre_basis(x,n,first,last);
    }
  };

 template<typename T,typename RandomAccessIterator>
  struct EvalChebyshevBasis: public EvalBasis<T,RandomAccessIterator>
  {
    void operator()(const T& x, 
                     const std::size_t n,
                     RandomAccessIterator first, 
                     RandomAccessIterator last ) const
    {
      slip::eval_chebyshev_basis(x,n,first,last);
    }
  };

  template<typename T,typename RandomAccessIterator>
  struct EvalPower2dBasis: public EvalBasis<T,RandomAccessIterator>
 {
   void operator()(const T& x, 
                     const std::size_t n,
                     RandomAccessIterator first, 
                     RandomAccessIterator last ) const
    {
      slip::eval_power_2d_basis(x,n,first,last);
    }
  
 };

 template<typename T,typename RandomAccessIterator>
  struct EvalPowerNdBasis: public EvalBasis<T,RandomAccessIterator>
 {

   void operator()(const T& x, 
                                   const std::size_t n,              
                                   RandomAccessIterator first, 
                                   RandomAccessIterator last ) const
    {
      slip::eval_power_nd_basis(x,n,first,last);
    }

   std::size_t dim_;
  
 };


}//::slip

#endif //SLIP_POLYNOMIAL_ALGO_HPP