Pade.h

Go to the documentation of this file.

#ifndef AMPGEN_PADE_H
#define AMPGEN_PADE_H
 
#include <array>
#include <functional>
#include <vector>
#include <iostream>
 
namespace AmpGen {
  enum Strategy { linear, quadratic, cubic, quartic };
  namespace detail {
    std::vector<double> solve_pade(const std::function<double(const double &)> &fcn, const double &min, const double &max, const unsigned &N,
                                   const Strategy &strat = Strategy::linear);
  }
 
  template <unsigned N, class T = double> class Pade {
  public:
    Pade(const std::vector<double> &r, const double &min, const double &max) : min(min), max(max) {
      for(unsigned i = 0; i <= N; ++i) co_f[i] = r[i];
      for(unsigned i = 0; i < N; ++i) co_g[i] = r[i + (N + 1)];
      range = 1. / (max - min);
    }
 
    Pade(const std::function<double(const double &)> &fcn, const double &min, const double &max, const Strategy &strat = Strategy::linear)
        : m_function(fcn), min(min), max(max) {
      auto r = detail::solve_pade(fcn, min, max, N, strat);
      for(unsigned i = 0; i <= N; ++i) co_f[i] = r[i];
      for(unsigned i = 0; i < N; ++i) co_g[i] = r[i + (N + 1)];
      range = 1. / (max - min);
    }
    template <typename T2> T2 operator()(const T2 &s) const {
      T2 x = (s - min) * range;
      T2 f = 0.;
      T2 g = 1.;
      T2 acc = 1.;
      for(unsigned i = 0; i < N; ++i) {
        f += co_f[i] * acc;
        acc *= x;
        g += co_g[i] * acc;
      }
      return (f + co_f[N] * acc) / g;
    }
    void print() const {
      for(int i = 0; i != N + 1; ++i) std::cout << co_f[i] << std::endl;
      for(int i = 0; i != N; ++i) std::cout << co_g[i] << std::endl;
    }
 
  private:
    std::function<double(const double &)> m_function;
    std::array<T, N + 1> co_f;
    std::array<T, N> co_g;
    T min;
    T max;
    T range;
  };
}
 
#endif