EhrhartPolynomial.h

 
#pragma once
 
#if defined(EhrhartPolynomial_RECURSES)
#error Recursive header files inclusion detected in EhrhartPolynomial.h
#else // defined(EhrhartPolynomial_RECURSES)
 
#define EhrhartPolynomial_RECURSES
 
#if !defined EhrhartPolynomial_h
 
#define EhrhartPolynomial_h
 
// Inclusions
#include <iostream>
#include <vector>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/NumberTraits.h"
#include "DGtal/kernel/CInteger.h"
#include "DGtal/kernel/CSpace.h"
#include "DGtal/math/LagrangeInterpolation.h"
#include "DGtal/geometry/volumes/BoundedLatticePolytope.h"
 
namespace DGtal
{
  template <typename TSpace, typename TInteger>
  class EhrhartPolynomial
  {
    BOOST_CONCEPT_ASSERT(( concepts::CSpace< TSpace > ));
    BOOST_CONCEPT_ASSERT(( concepts::CInteger< TInteger > ));
 
    // ----------------------- public types -----------------------------------
  public:
    typedef EhrhartPolynomial< TSpace, TInteger > Self;
    typedef TSpace                               Space;
    typedef TInteger                             Integer;
    typedef BoundedLatticePolytope<Space>        LatticePolytope;
    typedef LagrangeInterpolation< Integer >     Lagrange;
    typedef typename Lagrange::Polynomial        Polynomial;
    typedef std::size_t                          Size;
    
    // ----------------------- Standard services ------------------------------
  public:
 
    ~EhrhartPolynomial() = default;
 
    EhrhartPolynomial()
      : myE(), myD( NumberTraits< Integer >::ZERO )
    {}
 
    EhrhartPolynomial( const LatticePolytope& polytope )
    {
      init( polytope );
    }
                          
    EhrhartPolynomial( const EhrhartPolynomial & other ) = default;
 
    EhrhartPolynomial( EhrhartPolynomial&& other ) = default;
 
    EhrhartPolynomial & operator=( const EhrhartPolynomial & other ) = default;
 
    EhrhartPolynomial & operator=( EhrhartPolynomial&& other ) = default;
 
    inline Polynomial numerator() const
    {
      return myE;
    }
 
    inline Integer denominator() const
    {
      return myD;
    }
 
    void init( const LatticePolytope& polytope )
    {
      std::vector< Integer > X( Space::dimension + 1 );
      std::vector< Integer > Y( Space::dimension + 1 );
      X[ 0 ] = NumberTraits< Integer >::ZERO;
      Y[ 0 ] = NumberTraits< Integer >::ONE;
      for ( Integer i = 1; i <= Space::dimension; ++i )
        {
          LatticePolytope Q = i * polytope;
          Integer nb = Q.count();
          X[ i ] = i;
          Y[ i ] = nb;
        }
      // Compute polynomial (degree d)
      Lagrange L( X );
      myE = L.polynomial( Y );
      myD = L.denominator();
    }
 
    // @param t the dilation factor for this polytope P
    // @return the number of lattice points of t * P
    Integer count( Integer t ) const
    {
      return myE( t ) / myD;
    }
 
    // @param t the dilation factor for this polytope P
    // @return 0 if everything is correct.
    Integer remainder( Integer t ) const
    {
      return myE( t ) % myD;
    }
    
    // @param t the dilation factor for this polytope P
    // @return the number of lattice points interior to t * P
    // @note Use reciprocity formula (-1)^d E(-t)
    Integer countInterior( Integer t ) const { 
      return myE.degree() % 2 == 0 
        ?  ( myE( -t ) / myD ) 
        :  - ( myE( -t ) / myD ); 
    }
 
    // @param t the dilation factor for this polytope P
    // @return 0 if everything is correct.
    Integer remainderInterior( Integer t ) const { 
      return myE.degree() % 2 == 0 
        ?  ( myE( -t ) % myD ) 
        :  - ( myE( -t ) % myD ); 
    }
      
    // ----------------------- Interface --------------------------------------
  public:
 
    void selfDisplay( std::ostream & that_stream ) const
    {
      that_stream << "[EhrhartPolynomial num=" << numerator()
                  << " den=" << denominator() << " ]" << std::endl;
    }
 
    bool OK() const
    {
      return myD != NumberTraits< Integer >::ZERO;
    }
  
    // ------------------------- Datas ----------------------------------------
  protected:
 
    Polynomial myE;
    Integer    myD;
  
  };
 
  template <typename TSpace, typename TInteger>
  std::ostream&
  operator<<( std::ostream & that_stream, 
              const EhrhartPolynomial< TSpace, TInteger > & that_object_to_display )
  {
    that_object_to_display.selfDisplay( that_stream );
    return that_stream;
  }
  
  
} // namespace DGtal
 
 
// Includes inline functions/methods if necessary.
 
//                                                                           //
 
#endif // !defined EhrhartPolynomial_h
 
#undef EhrhartPolynomial_RECURSES
#endif // else defined(EhrhartPolynomial_RECURSES)