apfel/ome/laurent_polynomial.h

0001 /******************************************************************************
0002  * This file is part of libome                                                *
0003  * Copyright (C) 2025 Arnd Behring, Kay Schoenwald                            *
0004  * SPDX-License-Identifier: GPL-3.0-or-later                                  *
0005  ******************************************************************************/
0006
0007 /**
0008  * \file
0009  * \brief Evaluatable (and nestable) Laurent polynomials
0010  */
0011
0012 #ifndef LIBOME_LAURENT_POLYNOMIAL_H
0013 #define LIBOME_LAURENT_POLYNOMIAL_H
0014
0015 #include <utility>
0016 #include <type_traits>
0017 #include <initializer_list>
0018 #include <iterator>
0019 #include <vector>
0020 #include <cmath>
0021 #include <numeric>
0022 #include <apfel/ome/traits.h>
0023
0024 namespace apfel
0025 {
0026   namespace ome
0027   {
0028     /**
0029      * \brief View on an evaluatable Laurent polynomial
0030      *
0031      * \details
0032      * This class provides a view on a Laurent polyonmial (i.e., it does not store
0033      * the Laurent polynomial itself, but it only refers to one). It can be nested
0034      * (see documentation of \ref laurent_polynomial) and it can be evaluated.
0035      *
0036      * \tparam Tnum Numerical data type. This is the type of the polynomial
0037      *         variable \f$x\f$, the type of the coefficients or the return type
0038      *         of evaluating the coefficients \f$c_i(\dots)\f$ and the type that
0039      *         the polynomial evaluates to.
0040      * \tparam Tcoeff Type of the coefficients \f$c_i\f$ of the polynomial. This
0041      *         can either be a numerical type that can be implicitly converted to
0042      *         Tnum or a callable that evaluates to Tnum when invoked.
0043      * \tparam Tcoefffargs Pack of types that specifies the types of the arguments
0044      *         when invoking Tcoeff. If Tcoeff is a numerical type (not a
0045      *         callable) this pack should be empty.
0046      */
0047     template<typename Tnum, typename Tcoeff, typename... Tcoeffargs>
0048     class laurent_polynomial_view
0049     {
0050     public:
0051       /// Type alias for the numerical type template parameter
0052       using numeric_type = Tnum;
0053       /// Type alias for the coefficient type template parameter
0054       using coefficient_type = Tcoeff;
0055       /// Boolean type alias for whether the coefficient_type has a view
0056       using coefficient_has_view = has_view<coefficient_type>;
0057       /// Boolean type alias indicating that this class has an eval_plus_int method
0058       using has_eval_plus_int = std::true_type;
0059       /// Type alias for the view of the coefficient type (if it exists, otherwise it is void)
0060       using coefficient_view_type = typename get_view_type<coefficient_type>::type;
0061       /// Type alias for the coefficient iterators
0062       using const_reverse_iterator = typename std::vector<coefficient_type>::const_reverse_iterator;
0063
0064       /**
0065        * \brief Construct from const iterators
0066        *
0067        * \param crbegin Reverse iterator to last coefficient
0068        * \param crend Reverse iterator to behind first coefficient
0069        * \param min_power Exponent of the lowest-order monomial (\f$n_0\f$).
0070        */
0071       laurent_polynomial_view(
0072         const_reverse_iterator crbegin,
0073         const_reverse_iterator crend,
0074         int min_power = 0)
0075         : coeffs_crbegin_(crbegin),
0076           coeffs_crend_(crend),
0077           min_power_(min_power) {};
0078
0079       /**
0080        * \brief Swap function
0081        */
0082       friend
0083       void swap(laurent_polynomial_view& first, laurent_polynomial_view& second)
0084       {
0085         using std::swap;
0086         swap(first.coeffs_crbegin_,second.coeffs_crbegin_);
0087         swap(first.coeffs_crend_,second.coeffs_crend_);
0088         swap(first.min_power_,second.min_power_);
0089       };
0090
0091       /**
0092        * \brief Evaluate Laurent polynomial
0093        *
0094        * \details
0095        * Evaluate the polynomial for a given set of values of the variables. The
0096        * first argument is the value for the variable \f$x\f$. If there are more
0097        * arguments, they are passed on to the evaluation of each coefficient.
0098        *
0099        * \param x Value for the polynomial variable \f$x\f$.
0100        * \param rest Pack of function arguments that are passed to the
0101        *        evaluation of the coefficients \f$c_i(\dots)\f$.
0102        *
0103        * \return Numerical value for the Laurent polynomial evaluated at the
0104        *         given point.
0105        */
0106       numeric_type operator()(numeric_type x, Tcoeffargs... rest) const
0107       {
0108         using std::pow;
0109         const numeric_type prefactor = (min_power_ != 0) ? pow(x, min_power_) : static_cast<numeric_type>(1);
0110
0111         if constexpr(sizeof...(rest) > 0)
0112           {
0113             if(std::next(coeffs_crbegin_) == coeffs_crend_)
0114               {
0115                 return(prefactor*(*coeffs_crbegin_)(rest...));
0116               }
0117             else
0118               {
0119                 return(prefactor*std::accumulate(
0120                          coeffs_crbegin_,
0121                          coeffs_crend_,
0122                          static_cast<numeric_type>(0),
0123                          [&x, &rest...] (const numeric_type& acc, const coefficient_type& next)
0124                 {
0125                   return(x*acc + next(rest...));
0126                 }
0127                        ));
0128               }
0129           }
0130         else
0131           {
0132             if(std::next(coeffs_crbegin_) == coeffs_crend_)
0133               {
0134                 return(prefactor*(*coeffs_crbegin_));
0135               }
0136             else
0137               {
0138                 return(prefactor*std::accumulate(
0139                          coeffs_crbegin_,
0140                          coeffs_crend_,
0141                          static_cast<numeric_type>(0),
0142                          [&x] (const numeric_type& acc, const coefficient_type& next)
0143                 {
0144                   return(x*acc + next);
0145                 }
0146                        ));
0147               }
0148           }
0149       };
0150
0151       /**
0152        * \brief Evaluate polynomial with precomputed monomials
0153        *
0154        * \details
0155        * Evaluate the Laurent polynomial, but instead of the \f$x^i\f$, the
0156        * monomials are taken from the input vector subst. Effectively this method computes the inner product between subst and the vector of polynomial coefficients.
0157        *
0158        * \param subst Vector of precomputed monomials. The caller has to
0159        *        ensure that it has the same length as the vector of
0160        *        polynomial coefficients
0161        * \param rest Pack of function arguments that are passed to the
0162        *        evaluation of the coefficients \f$c_i(\dots)\f$.
0163        *
0164        * \return Numerical value for the Laurent polynomial evaluated with the
0165        *         given monomials.
0166        */
0167       numeric_type eval_subst(const std::vector<numeric_type>& subst, Tcoeffargs... rest) const
0168       {
0169         if constexpr(sizeof...(rest) > 0)
0170           {
0171             return(std::transform_reduce(
0172                      coeffs_crbegin_,
0173                      coeffs_crend_,
0174                      subst.crbegin(),
0175                      static_cast<numeric_type>(0),
0176                      std::plus<>(),
0177                      [&rest...] (const coefficient_type& c, const numeric_type& s)
0178             {
0179               return(c(rest...) * s);
0180             }
0181                    ));
0182           }
0183         else
0184           {
0185             return(std::transform_reduce(
0186                      coeffs_crbegin_,
0187                      coeffs_crend_,
0188                      subst.crbegin(),
0189                      static_cast<numeric_type>(0)
0190                    ));
0191           }
0192       };
0193
0194       /**
0195        * \brief Evaluate polynomial, but call eval_plus_int on wrapped type
0196        *
0197        * \details
0198        * This is meant to be used in conjunction with func_plusfunc_omx, which
0199        * allows to calculate the "extra" integral of plus functions over
0200        * \f$[0,x]\f$ analytically. This method exposes that functionality if
0201        * func_plusfunc_omx is wrapped into further laurent_polynomials.
0202        *
0203        * \param x Value for the polynomial variable \f$x\f$.
0204        * \param rest Pack of function arguments that are passed to the
0205        *        evaluation of the coefficients \f$c_i(\dots)\f$.
0206        *
0207        * \return Numerical value for the Laurent polynomial evaluated for the
0208        *         given values of the variables, but with eval_plus_int called
0209        *         on the wrapped callable.
0210        */
0211       numeric_type eval_plus_int(numeric_type x, Tcoeffargs... rest) const
0212       {
0213         using std::pow;
0214         const numeric_type prefactor = (min_power_ != 0) ? pow(x, min_power_) : static_cast<numeric_type>(1);
0215
0216         if(std::next(coeffs_crbegin_) == coeffs_crend_)
0217           {
0218             return(prefactor * coeffs_crbegin_->eval_plus_int(rest...));
0219           }
0220         else
0221           {
0222             return(prefactor*std::accumulate(
0223                      coeffs_crbegin_,
0224                      coeffs_crend_,
0225                      static_cast<numeric_type>(0),
0226                      [&x, &rest...] (const numeric_type& acc, const coefficient_type& next)
0227             {
0228               return(x*acc + next.eval_plus_int(rest...));
0229             }
0230                    ));
0231           }
0232       };
0233
0234       /**
0235        * \brief Truncate the Laurent polynomial to a given order
0236        *
0237        * \details
0238        * Truncation means that the polynomial is only evaluated up to a certain
0239        * order, which can by lower than the highest available term. If the
0240        * requested truncation order is lower than the minimal exponent
0241        * (\f$n_0\f$), the resulting truncated polynomial will always evaluate to
0242        * zero. If the requested truncation order is higher than the maximal
0243        * exponent, no truncation will be performed.
0244        *
0245        * \param truncation_order Exponent of the last order to include
0246        *
0247        * \return A laurent_polynomial_view of the truncated polynomial
0248        */
0249       laurent_polynomial_view truncate(int truncation_order) const
0250       {
0251         const int requested_orders = truncation_order - min_power_ + 1;
0252         const int available_orders = coeffs_crend_ - coeffs_crbegin_;
0253
0254         if(requested_orders <= 0)
0255           {
0256             return(laurent_polynomial_view(coeffs_crbegin_,
0257                                            coeffs_crbegin_, min_power_));
0258           }
0259         else if(requested_orders > available_orders)
0260           {
0261             return(laurent_polynomial_view(coeffs_crbegin_,
0262                                            coeffs_crend_, min_power_));
0263           }
0264         else
0265           {
0266             return(laurent_polynomial_view(
0267                      coeffs_crbegin_ + (available_orders - requested_orders),
0268                      coeffs_crend_,
0269                      min_power_
0270                    ));
0271           }
0272       };
0273
0274       /**
0275        * \brief Clone the view on this Laurent polynomial
0276        *
0277        * \return A laurent_polynomial_view on this polynomial
0278        */
0279       laurent_polynomial_view get_view() const
0280       {
0281         return(laurent_polynomial_view(*this));
0282       };
0283
0284       /**
0285        * \brief Get the minimum exponent \f$n_0\f$ of the polynomial
0286        *
0287        * \return Minimum exponent
0288        */
0289       int min_power() const
0290       {
0291         return(min_power_);
0292       };
0293
0294       /**
0295        * \brief Get the maximum exponent \f$N\f$ of the polynomial
0296        *
0297        * \details
0298        * If the polynomial is empty (i.e. it doesn't have any coefficients)
0299        * it returns min_power - 1
0300        *
0301        * \return Maximum exponent
0302        */
0303       int max_power() const
0304       {
0305         return(min_power_ + (coeffs_crend_ - coeffs_crbegin_ - 1));
0306       };
0307
0308       /**
0309        * \brief Access to the coefficients of the Laurent polynomial
0310        *
0311        * \param pow Index of the coefficient to access
0312        *
0313        * \return Constant reference to the coefficient of \f$x^\mathtt{pow}\f$.
0314        */
0315       const coefficient_type& operator[](int pow) const
0316       {
0317         const int requested_position = pow - min_power_;
0318         const int size = coeffs_crend_ - coeffs_crbegin_;
0319         if(requested_position < 0 || requested_position >= size)
0320           return zero_value_;
0321         else
0322           return *(coeffs_crbegin_ + ((size-1) - requested_position));
0323       };
0324
0325     protected:
0326       /// Reverse iterator pointing to the beginning of the coefficient vector
0327       const_reverse_iterator coeffs_crbegin_;
0328       /// Reverse iterator pointing to the end of the coefficient vector
0329       const_reverse_iterator coeffs_crend_;
0330
0331       /**
0332        * \brief Default constructor
0333        */
0334       laurent_polynomial_view(int min_power = 0)
0335         : min_power_(min_power) {};
0336
0337     private:
0338       int min_power_;
0339
0340       inline const static coefficient_type zero_value_{};
0341     };
0342
0343     /**
0344      * \brief Potentially nested, evaluatable Laurent polynomial
0345      *
0346      * \details
0347      * This class implements a Laurent polynomial \f$p(x, \dots)\f$ (i.e. it can
0348      * have monomials with negative degree).
0349      * \f[
0350      *    p(x, \dots) = \sum_{n=n_0}^N x^i c_i(\dots)
0351      * \f]
0352      * The polynomial can be nested, which means that its coefficients \f$c_i\f$
0353      * can be callables. In particular, they can be laurent_polynomials. This
0354      * allows to implement multivariate polynomials. Coefficients are specified
0355      * at the time of construction and cannot be modified.
0356      *
0357      * \tparam Tnum Numerical data type. This is the type of the polynomial
0358      *         variable \f$x\f$, the type of the coefficients or the return type
0359      *         of evaluating the coefficients \f$c_i(\dots)\f$ and the type that
0360      *         the polynomial evaluates to.
0361      * \tparam Tcoeff Type of the coefficients \f$c_i\f$ of the polynomial. This
0362      *         can either be a numerical type that can be implicitly converted to
0363      *         Tnum or a callable that evaluates to Tnum when invoked.
0364      * \tparam Tcoefffargs Pack of types that specifies the types of the arguments
0365      *         when invoking Tcoeff. If Tcoeff is a numerical type (not a
0366      *         callable) this pack should be empty.
0367      */
0368     template<typename Tnum, typename Tcoeff, typename... Tcoeffargs>
0369     class laurent_polynomial
0370       : public laurent_polynomial_view<Tnum, Tcoeff, Tcoeffargs...>
0371     {
0372     public:
0373       /// Type alias for the underlying view from which we inherit
0374       using view_type = laurent_polynomial_view<Tnum, Tcoeff, Tcoeffargs...>;
0375       // Make type aliases from base class usable here
0376       using typename view_type::numeric_type;
0377       using typename view_type::coefficient_type;
0378
0379
0380       /**
0381        * \brief Construct a null polynomial (always evaluates to zero)
0382        */
0383       laurent_polynomial()
0384         : view_type(0),
0385           coefficients_()
0386       {
0387         this->coeffs_crbegin_ = coefficients_.crbegin();
0388         this->coeffs_crend_ = coefficients_.crend();
0389       };
0390
0391       /**
0392        * \brief Construct Laurent polynomial from a vector of coefficients
0393        *
0394        * \param coefficients Vector of coefficients \f$c_i(\dots)\f$ of the
0395        *        polynomial. The first element is the coefficient of the
0396        *        lowest-order monomial \f$x^{n_0}\f$ and subsequent elements are
0397        *        the coefficients of the higher powers.
0398        * \param min_power Exponent of the lowest-order monomial (\f$n_0\f$).
0399        */
0400       laurent_polynomial(const std::vector<coefficient_type>& coefficients,
0401                          int min_power = 0)
0402         : view_type(min_power),
0403           coefficients_(coefficients)
0404       {
0405         this->coeffs_crbegin_ = coefficients_.crbegin();
0406         this->coeffs_crend_ = coefficients_.crend();
0407       };
0408
0409       /**
0410        * \brief Helper to construct a Laurent polynomial directly from an
0411        *        initalizer list.
0412        *
0413        * \details
0414        * See \ref laurent_polynomial(const std::vector<coefficient_type>&, int) for details.
0415        */
0416       laurent_polynomial(std::initializer_list<coefficient_type> coefficients,
0417                          int min_power = 0)
0418         : view_type(min_power),
0419           coefficients_(coefficients)
0420       {
0421         this->coeffs_crbegin_ = coefficients_.crbegin();
0422         this->coeffs_crend_ = coefficients_.crend();
0423       };
0424
0425       /**
0426        * \brief Copy constructor
0427        */
0428       laurent_polynomial(const laurent_polynomial& other)
0429         : view_type(other),
0430           coefficients_(other.coefficients_)
0431       {
0432         this->coeffs_crbegin_ = coefficients_.crbegin();
0433         this->coeffs_crend_ = coefficients_.crend();
0434       };
0435
0436       /**
0437        * \brief Swap function
0438        */
0439       friend
0440       void swap(laurent_polynomial& first, laurent_polynomial& second)
0441       {
0442         // enable proper ADL
0443         using std::swap;
0444
0445         // swap base class members
0446         swap(static_cast<view_type&>(first),static_cast<view_type&>(second));
0447
0448         // swap local data members
0449         swap(first.coefficients_,second.coefficients_);
0450
0451         // update iterators
0452         first.coeffs_crbegin_ = first.coefficients_.crbegin();
0453         first.coeffs_crend_   = first.coefficients_.crend();
0454         second.coeffs_crbegin_ = second.coefficients_.crbegin();
0455         second.coeffs_crend_   = second.coefficients_.crend();
0456       };
0457
0458       /**
0459        * \brief Assignment operator
0460        */
0461       laurent_polynomial& operator=(laurent_polynomial other)
0462       {
0463         swap(*this, other);
0464         return(*this);
0465       };
0466
0467       /**
0468        * \brief Move operator
0469        */
0470       laurent_polynomial(laurent_polynomial&& other)
0471         : laurent_polynomial()
0472       {
0473         swap(*this, other);
0474       };
0475
0476     private:
0477       std::vector<coefficient_type> coefficients_;
0478     };
0479
0480   }
0481 }
0482
0483 #endif