LagrangeInterpolation.h

 
#pragma once
 
#if defined(LagrangeInterpolation_RECURSES)
#error Recursive header files inclusion detected in LagrangeInterpolation.h
#else // defined(LagrangeInterpolation_RECURSES)
#define LagrangeInterpolation_RECURSES
 
#if !defined LagrangeInterpolation_h
#define LagrangeInterpolation_h
 
// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"
#include "DGtal/math/MPolynomial.h"
#include "DGtal/kernel/NumberTraits.h"
#include "DGtal/kernel/CEuclideanRing.h"
#include <vector>
 
namespace DGtal
{
  template <typename TEuclideanRing>
  class LagrangeInterpolation
  {
 
    // ----------------------- public types -----------------------------------
  public:
    typedef LagrangeInterpolation< TEuclideanRing > Self;
    typedef TEuclideanRing Ring;
    BOOST_CONCEPT_ASSERT(( concepts::CEuclideanRing< Ring > ) );
 
    typedef DGtal::MPolynomial< 1, Ring > Polynomial;
    typedef std::size_t                   Size;
    
    // ----------------------- Standard services ------------------------------
  public:
 
    ~LagrangeInterpolation() = default;
 
    LagrangeInterpolation()
      : myDeterminant( NumberTraits< Ring >::ZERO )
    {}
 
    LagrangeInterpolation( const std::vector< Ring >& xvalues )
    {
      init( xvalues );
    }
                          
    LagrangeInterpolation( const LagrangeInterpolation & other ) = default;
 
    LagrangeInterpolation( LagrangeInterpolation&& other ) = default;
 
    LagrangeInterpolation & operator=( const LagrangeInterpolation & other ) = default;
 
    LagrangeInterpolation & operator=( LagrangeInterpolation&& other ) = default;
 
    inline Size size() const
    {
      return myX.size();
    }
 
    inline Size degree() const
    {
      return size() - 1;
    }
 
    inline Ring denominator() const
    {
      return myDeterminant;
    }
    
    void init( const std::vector< Ring >& xvalues )
    {
      myX           = xvalues;
      // Compute Vandermonde determinant
      myDeterminant = NumberTraits< Ring >::ONE;
      for ( Dimension i = 0; (i+1) < size(); ++i )
        for ( Dimension j = i + 1; j < size(); ++j )
          myDeterminant *= myX[ j ] - myX[ i ];
      // compute Lagrange polynomial basis.
      myLagrangeBasis.clear();
      myLagrangeBasis.reserve( size() );
      for ( Dimension j = 0; j < size(); ++j )
        {
          Ring c = myDeterminant;
          for ( Dimension m = 0; m < size(); ++m )
            if ( m != j ) 
              c /= myX[ j ] - myX[ m ];
          Polynomial P = c; // constant
          for ( Dimension m = 0; m < size(); ++m )
            if ( m != j ) 
              P *= mmonomial<Ring>( 1 ) - myX[ m ] * mmonomial<Ring>( 0 );
          myLagrangeBasis.push_back( P );
        }
    }
 
    Polynomial polynomial( const std::vector< Ring >& yvalues )
    {
      Polynomial P;
      if ( yvalues.size() != size() ) return P;
      for ( Dimension j = 0; j < size(); ++j )
        P += yvalues[ j ] * myLagrangeBasis[ j ];
      return P;
    }
    
    // ----------------------- Accessors ------------------------------
  public:
 
    Polynomial basis( Size i ) const
    {
      ASSERT( 0 <= i && i < size() );
      return myLagrangeBasis[ i ];
    }
      
    // ----------------------- Interface --------------------------------------
  public:
 
    void selfDisplay( std::ostream & that_stream ) const
    {
      that_stream << "[LagrangeInterpolation det=" << myDeterminant << std::endl;
      for ( Size i = 0; i < size(); i++ )
        that_stream << "l_" << i << "=" << basis( i ) << std::endl;
      that_stream << "]";
    }
 
    bool OK() const
    {
      return myDeterminant != NumberTraits< Ring >::ZERO;
    }
  
    // ------------------------- Datas ----------------------------------------
  protected:
 
    std::vector< Ring > myX;
    Ring                myDeterminant;
    std::vector<Polynomial> myLagrangeBasis;
  
  };
 
  template <typename TEuclideanRing>
  std::ostream&
  operator<<( std::ostream & that_stream, 
        const LagrangeInterpolation< TEuclideanRing > & that_object_to_display )
  {
    that_object_to_display.selfDisplay( that_stream );
    return that_stream;
  }
  
  
} // namespace DGtal
 
 
// Includes inline functions/methods if necessary.
 
//                                                                           //
 
#endif // !defined LagrangeInterpolation_h
 
#undef LagrangeInterpolation_RECURSES
#endif // else defined(LagrangeInterpolation_RECURSES)