///////////////////////////////////////////////////////////////////////////////
//  Copyright 2018 John Maddock
//  Distributed under the Boost
//  Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_
#define BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_

#ifndef BOOST_MATH_PFQ_MAX_B_TERMS
#  define BOOST_MATH_PFQ_MAX_B_TERMS 5
#endif

#include <array>
#include <cstdint>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/expm1.hpp>
#include <boost/math/special_functions/detail/hypergeometric_series.hpp>

  namespace boost { namespace math { namespace detail {

     template <class Seq, class Real>
     unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations)
     {
        BOOST_MATH_STD_USING
        unsigned N_terms = 0;

        if(aj.size() == 1 && bj.size() == 1)
        {
           //
           // For 1F1 we can work out where the peaks in the series occur,
           //  which is to say when:
           //
           // (a + k)z / (k(b + k)) == +-1
           //
           // Then we are at either a maxima or a minima in the series, and the
           // last such point must be a maxima since the series is globally convergent.
           // Potentially then we are solving 2 quadratic equations and have up to 4
           // solutions, any solutions which are complex or negative are discarded,
           // leaving us with 4 options:
           //
           // 0 solutions: The series is directly convergent.
           // 1 solution : The series diverges to a maxima before converging.
           // 2 solutions: The series is initially convergent, followed by divergence to a maxima before final convergence.
           // 3 solutions: The series diverges to a maxima, converges to a minima before diverging again to a second maxima before final convergence.
           // 4 solutions: The series converges to a minima before diverging to a maxima, converging to a minima, diverging to a second maxima and then converging.
           //
           // The first 2 situations are adequately handled by direct series evaluation, while the 2,3 and 4 solutions are not.
           //
           Real a = *aj.begin();
           Real b = *bj.begin();
           Real sq = 4 * a * z + b * b - 2 * b * z + z * z;
           if (sq >= 0)
           {
              Real t = (-sqrt(sq) - b + z) / 2;
              if (t >= 0)
              {
                 crossover_locations[N_terms] = itrunc(t);
                 ++N_terms;
              }
              t = (sqrt(sq) - b + z) / 2;
              if (t >= 0)
              {
                 crossover_locations[N_terms] = itrunc(t);
                 ++N_terms;
              }
           }
           sq = -4 * a * z + b * b + 2 * b * z + z * z;
           if (sq >= 0)
           {
              Real t = (-sqrt(sq) - b - z) / 2;
              if (t >= 0)
              {
                 crossover_locations[N_terms] = itrunc(t);
                 ++N_terms;
              }
              t = (sqrt(sq) - b - z) / 2;
              if (t >= 0)
              {
                 crossover_locations[N_terms] = itrunc(t);
                 ++N_terms;
              }
           }
           std::sort(crossover_locations, crossover_locations + N_terms, std::less<Real>());
           //
           // Now we need to discard every other terms, as these are the minima:
           //
           switch (N_terms)
           {
           case 0:
           case 1:
              break;
           case 2:
              crossover_locations[0] = crossover_locations[1];
              --N_terms;
              break;
           case 3:
              crossover_locations[1] = crossover_locations[2];
              --N_terms;
              break;
           case 4:
              crossover_locations[0] = crossover_locations[1];
              crossover_locations[1] = crossover_locations[3];
              N_terms -= 2;
              break;
           }
        }
        else
        {
           unsigned n = 0;
           for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n)
           {
              crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1;
           }
           std::sort(crossover_locations, crossover_locations + bj.size(), std::less<Real>());
           N_terms = (unsigned)bj.size();
        }
        return N_terms;
     }

     template <class Seq, class Real, class Policy, class Terminal>
     std::pair<Real, Real> hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, long long& log_scale)
     {
        BOOST_MATH_STD_USING
        Real result = 1;
        Real abs_result = 1;
        Real term = 1;
        Real term0 = 0;
        Real tol = boost::math::policies::get_epsilon<Real, Policy>();
        std::uintmax_t k = 0;
        Real upper_limit(sqrt(boost::math::tools::max_value<Real>())), diff;
        Real lower_limit(1 / upper_limit);
        long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>()) - 2;
        Real scaling_factor = exp(Real(log_scaling_factor));
        Real term_m1;
        long long local_scaling = 0;
        bool have_no_correct_bits = false;

        if ((aj.size() == 1) && (bj.size() == 0))
        {
           if (fabs(z) > 1)
           {
              if ((z > 0) && (floor(*aj.begin()) != *aj.begin()))
              {
                 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == 1 and q == 0 and |z| > 1, result is imaginary", z, pol);
                 return std::make_pair(r, r);
              }
              std::pair<Real, Real> r = hypergeometric_pFq_checked_series_impl(aj, bj, Real(1 / z), pol, termination, log_scale);
              
              #if (defined(__GNUC__) && __GNUC__ == 13)
              Real mul = pow(-z, Real(-*aj.begin()));
              #else
              Real mul = pow(-z, -*aj.begin());
              #endif
              
              r.first *= mul;
              r.second *= mul;
              return r;
           }
        }

        if (aj.size() > bj.size())
        {
           if (aj.size() == bj.size() + 1)
           {
              if (fabs(z) > 1)
              {
                 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| > 1, series does not converge", z, pol);
                 return std::make_pair(r, r);
              }
              if (fabs(z) == 1)
              {
                 Real s = 0;
                 for (auto i = bj.begin(); i != bj.end(); ++i)
                    s += *i;
                 for (auto i = aj.begin(); i != aj.end(); ++i)
                    s -= *i;
                 if ((z == 1) && (s <= 0))
                 {
                    Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| == 1, in a situation where the series does not converge", z, pol);
                    return std::make_pair(r, r);
                 }
                 if ((z == -1) && (s <= -1))
                 {
                    Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| == 1, in a situation where the series does not converge", z, pol);
                    return std::make_pair(r, r);
                 }
              }
           }
           else
           {
              Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p > q+1, series does not converge", z, pol);
              return std::make_pair(r, r);
           }
        }

        while (!termination(k))
        {
           for (auto ai = aj.begin(); ai != aj.end(); ++ai)
           {
              term *= *ai + k;
           }
           if (term == 0)
           {
              // There is a negative integer in the aj's:
              return std::make_pair(result, abs_result);
           }
           for (auto bi = bj.begin(); bi != bj.end(); ++bi)
           {
              if (*bi + k == 0)
              {
                 // The series is undefined:
                 result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
                 return std::make_pair(result, result);
              }
              term /= *bi + k;
           }
           term *= z;
           ++k;
           term /= k;
           //std::cout << k << " " << *bj.begin() + k << " " << result << " " << term << /*" " << term_at_k(*aj.begin(), *bj.begin(), z, k, pol) <<*/ std::endl;
           result += term;
           abs_result += abs(term);
           //std::cout << "k = " << k << " term = " << term * exp(log_scale) << " result = " << result * exp(log_scale) << " abs_result = " << abs_result * exp(log_scale) << std::endl;

           //
           // Rescaling:
           //
           if (fabs(abs_result) >= upper_limit)
           {
              abs_result /= scaling_factor;
              result /= scaling_factor;
              term /= scaling_factor;
              log_scale += log_scaling_factor;
              local_scaling += log_scaling_factor;
           }
           if (fabs(abs_result) < lower_limit)
           {
              abs_result *= scaling_factor;
              result *= scaling_factor;
              term *= scaling_factor;
              log_scale -= log_scaling_factor;
              local_scaling -= log_scaling_factor;
           }

           if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term)))
              break;
           if (abs_result * tol > abs(result))
           {
              // Check if result is so small compared to abs_resuslt that there are no longer any
              // correct bits... we require two consecutive passes here before aborting to
              // avoid false positives when result transiently drops to near zero then rebounds.
              if (have_no_correct_bits)
              {
                 // We have no correct bits in the result... just give up!
                 result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the result are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol);
                 return std::make_pair(result, result);
              }
              else
                 have_no_correct_bits = true;
           }
           else
              have_no_correct_bits = false;
           term0 = term;
        }
        //std::cout << "result = " << result << std::endl;
        //std::cout << "local_scaling = " << local_scaling << std::endl;
        //std::cout << "Norm result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(result) * exp(boost::multiprecision::mpfr_float_50(local_scaling)) << std::endl;
        //
        // We have to be careful when one of the b's crosses the origin:
        //
        if(bj.size() > BOOST_MATH_PFQ_MAX_B_TERMS)
           policies::raise_domain_error<Real>("boost::math::hypergeometric_pFq<%1%>(Seq, Seq, %1%)",
              "The number of b terms must be less than the value of BOOST_MATH_PFQ_MAX_B_TERMS (" BOOST_MATH_STRINGIZE(BOOST_MATH_PFQ_MAX_B_TERMS)  "), but got %1%.",
              Real(bj.size()), pol);

        unsigned crossover_locations[BOOST_MATH_PFQ_MAX_B_TERMS];

        unsigned N_crossovers = set_crossover_locations(aj, bj, z, crossover_locations);

        bool terminate = false;   // Set to true if one of the a's passes through the origin and terminates the series.

        for (unsigned n = 0; n < N_crossovers; ++n)
        {
           if (k < crossover_locations[n])
           {
              for (auto ai = aj.begin(); ai != aj.end(); ++ai)
              {
                 if ((*ai < 0) && (floor(*ai) == *ai) && (*ai > static_cast<decltype(*ai)>(crossover_locations[n])))
                    return std::make_pair(result, abs_result);  // b's will never cross the origin!
              }
              //
              // local results:
              //
              Real loop_result = 0;
              Real loop_abs_result = 0;
              long long loop_scale = 0;
              //
              // loop_error_scale will be used to increase the size of the error
              // estimate (absolute sum), based on the errors inherent in calculating
              // the pochhammer symbols.
              //
              Real loop_error_scale = 0;
              //boost::multiprecision::mpfi_float err_est = 0;
              //
              // b hasn't crossed the origin yet and the series may spring back into life at that point
              // so we need to jump forward to that term and then evaluate forwards and backwards from there:
              //
              unsigned s = crossover_locations[n];
              std::uintmax_t backstop = k;
              long long s1(1), s2(1);
              term = 0;
              for (auto ai = aj.begin(); ai != aj.end(); ++ai)
              {
                 if ((floor(*ai) == *ai) && (*ai < 0) && (-*ai <= static_cast<decltype(*ai)>(s)))
                 {
                    // One of the a terms has passed through zero and terminated the series:
                    terminate = true;
                    break;
                 }
                 else
                 {
                    int ls = 1;
                    Real p = log_pochhammer(*ai, s, pol, &ls);
                    s1 *= ls;
                    term += p;
                    loop_error_scale = (std::max)(p, loop_error_scale);
                    //err_est += boost::multiprecision::mpfi_float(p);
                 }
              }
              //std::cout << "term = " << term << std::endl;
              if (terminate)
                 break;
              for (auto bi = bj.begin(); bi != bj.end(); ++bi)
              {
                 int ls = 1;
                 Real p = log_pochhammer(*bi, s, pol, &ls);
                 s2 *= ls;
                 term -= p;
                 loop_error_scale = (std::max)(p, loop_error_scale);
                 //err_est -= boost::multiprecision::mpfi_float(p);
              }
              //std::cout << "term = " << term << std::endl;
              Real p = lgamma(Real(s + 1), pol);
              term -= p;
              loop_error_scale = (std::max)(p, loop_error_scale);
              //err_est -= boost::multiprecision::mpfi_float(p);
              p = s * log(fabs(z));
              term += p;
              loop_error_scale = (std::max)(p, loop_error_scale);
              //err_est += boost::multiprecision::mpfi_float(p);
              //err_est = exp(err_est);
              //std::cout << err_est << std::endl;
              //
              // Convert loop_error scale to the absolute error
              // in term after exp is applied:
              //
              loop_error_scale *= tools::epsilon<Real>();
              //
              // Convert to relative error after exp:
              //
              loop_error_scale = fabs(expm1(loop_error_scale, pol));
              //
              // Convert to multiplier for the error term:
              //
              loop_error_scale /= tools::epsilon<Real>();

              if (z < 0)
                 s1 *= (s & 1 ? -1 : 1);

              if (term <= tools::log_min_value<Real>())
              {
                 // rescale if we can:
                 long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2);
                 term -= scale;
                 loop_scale += scale;
              }
               if (term > 10)
               {
                  int scale = itrunc(floor(term));
                  term -= scale;
                  loop_scale += scale;
               }
               //std::cout << "term = " << term << std::endl;
               term = s1 * s2 * exp(term);
               //std::cout << "term = " << term << std::endl;
               //std::cout << "loop_scale = " << loop_scale << std::endl;
               k = s;
               term0 = term;
               long long saved_loop_scale = loop_scale;
               bool terms_are_growing = true;
               bool trivial_small_series_check = false;
               do
               {
                  loop_result += term;
                  loop_abs_result += fabs(term);
                  //std::cout << "k = " << k << " term = " << term * exp(loop_scale) << " result = " << loop_result * exp(loop_scale) << " abs_result = " << loop_abs_result * exp(loop_scale) << std::endl;
                  if (fabs(loop_result) >= upper_limit)
                  {
                     loop_result /= scaling_factor;
                     loop_abs_result /= scaling_factor;
                     term /= scaling_factor;
                     loop_scale += log_scaling_factor;
                  }
                  if (fabs(loop_result) < lower_limit)
                  {
                     loop_result *= scaling_factor;
                     loop_abs_result *= scaling_factor;
                     term *= scaling_factor;
                     loop_scale -= log_scaling_factor;
                  }
                  term_m1 = term;
                  for (auto ai = aj.begin(); ai != aj.end(); ++ai)
                  {
                     term *= *ai + k;
                  }
                  if (term == 0)
                  {
                     // There is a negative integer in the aj's:
                     return std::make_pair(result, abs_result);
                  }
                  for (auto bi = bj.begin(); bi != bj.end(); ++bi)
                  {
                     if (*bi + k == 0)
                     {
                        // The series is undefined:
                        result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
                        return std::make_pair(result, result);
                     }
                     term /= *bi + k;
                  }
                  term *= z / (k + 1);

                  ++k;
                  diff = fabs(term / loop_result);
                  terms_are_growing = fabs(term) > fabs(term_m1);
                  if (!trivial_small_series_check && !terms_are_growing)
                  {
                     //
                     // Now that we have started to converge, check to see if the value of
                     // this local sum is trivially small compared to the result.  If so
                     // abort this part of the series.
                     //
                     trivial_small_series_check = true;
                     Real d;
                     if (loop_scale > local_scaling)
                     {
                        long long rescale = local_scaling - loop_scale;
                        if (rescale < tools::log_min_value<Real>())
                           d = 1;  // arbitrary value, we want to keep going
                        else
                           d = fabs(term / (result * exp(Real(rescale))));
                     }
                     else
                     {
                        long long rescale = loop_scale - local_scaling;
                        if (rescale < tools::log_min_value<Real>())
                           d = 0;  // terminate this loop
                        else
                           d = fabs(term * exp(Real(rescale)) / result);
                     }
                     if (d < boost::math::policies::get_epsilon<Real, Policy>())
                        break;
                  }
               } while (!termination(k - s) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || terms_are_growing));

               //std::cout << "Norm loop result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(loop_result)* exp(boost::multiprecision::mpfr_float_50(loop_scale)) << std::endl;
               //
               // We now need to combine the results of the first series summation with whatever
               // local results we have now.  First though, rescale abs_result by loop_error_scale
               // to factor in the error in the pochhammer terms at the start of this block:
               //
               std::uintmax_t next_backstop = k;
               loop_abs_result += loop_error_scale * fabs(loop_result);
               if (loop_scale > local_scaling)
               {
                  //
                  // Need to shrink previous result:
                  //
                  long long rescale = local_scaling - loop_scale;
                  local_scaling = loop_scale;
                  log_scale -= rescale;
                  Real ex = exp(Real(rescale));
                  result *= ex;
                  abs_result *= ex;
                  result += loop_result;
                  abs_result += loop_abs_result;
               }
               else if (local_scaling > loop_scale)
               {
                  //
                  // Need to shrink local result:
                  //
                  long long rescale = loop_scale - local_scaling;
                  Real ex = exp(Real(rescale));
                  loop_result *= ex;
                  loop_abs_result *= ex;
                  result += loop_result;
                  abs_result += loop_abs_result;
               }
               else
               {
                  result += loop_result;
                  abs_result += loop_abs_result;
               }
               //
               // Now go backwards as well:
               //
               k = s;
               term = term0;
               loop_result = 0;
               loop_abs_result = 0;
               loop_scale = saved_loop_scale;
               trivial_small_series_check = false;
               do
               {
                  --k;
                  if (k == backstop)
                     break;
                  term_m1 = term;
                  for (auto ai = aj.begin(); ai != aj.end(); ++ai)
                  {
                     term /= *ai + k;
                  }
                  for (auto bi = bj.begin(); bi != bj.end(); ++bi)
                  {
                     if (*bi + k == 0)
                     {
                        // The series is undefined:
                        result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol);
                        return std::make_pair(result, result);
                     }
                     term *= *bi + k;
                  }
                  term *= (k + 1) / z;
                  loop_result += term;
                  loop_abs_result += fabs(term);

                  if (!trivial_small_series_check && (fabs(term) < fabs(term_m1)))
                  {
                     //
                     // Now that we have started to converge, check to see if the value of
                     // this local sum is trivially small compared to the result.  If so
                     // abort this part of the series.
                     //
                     trivial_small_series_check = true;
                     Real d;
                     if (loop_scale > local_scaling)
                     {
                        long long rescale = local_scaling - loop_scale;
                        if (rescale < tools::log_min_value<Real>())
                           d = 1;  // keep going
                        else
                           d = fabs(term / (result * exp(Real(rescale))));
                     }
                     else
                     {
                        long long rescale = loop_scale - local_scaling;
                        if (rescale < tools::log_min_value<Real>())
                           d = 0;  // stop, underflow
                        else
                           d = fabs(term * exp(Real(rescale)) / result);
                     }
                     if (d < boost::math::policies::get_epsilon<Real, Policy>())
                        break;
                  }

                  //std::cout << "k = " << k << " result = " << result << " abs_result = " << abs_result << std::endl;
                  if (fabs(loop_result) >= upper_limit)
                  {
                     loop_result /= scaling_factor;
                     loop_abs_result /= scaling_factor;
                     term /= scaling_factor;
                     loop_scale += log_scaling_factor;
                  }
                  if (fabs(loop_result) < lower_limit)
                  {
                     loop_result *= scaling_factor;
                     loop_abs_result *= scaling_factor;
                     term *= scaling_factor;
                     loop_scale -= log_scaling_factor;
                  }
                  diff = fabs(term / loop_result);
               } while (!termination(s - k) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || (fabs(term) > fabs(term_m1))));

               //std::cout << "Norm loop result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(loop_result)* exp(boost::multiprecision::mpfr_float_50(loop_scale)) << std::endl;
               //
               // We now need to combine the results of the first series summation with whatever
               // local results we have now.  First though, rescale abs_result by loop_error_scale
               // to factor in the error in the pochhammer terms at the start of this block:
               //
               loop_abs_result += loop_error_scale * fabs(loop_result);
               //
               if (loop_scale > local_scaling)
               {
                  //
                  // Need to shrink previous result:
                  //
                  long long rescale = local_scaling - loop_scale;
                  local_scaling = loop_scale;
                  log_scale -= rescale;
                  Real ex = exp(Real(rescale));
                  result *= ex;
                  abs_result *= ex;
                  result += loop_result;
                  abs_result += loop_abs_result;
               }
               else if (local_scaling > loop_scale)
               {
                  //
                  // Need to shrink local result:
                  //
                  long long rescale = loop_scale - local_scaling;
                  Real ex = exp(Real(rescale));
                  loop_result *= ex;
                  loop_abs_result *= ex;
                  result += loop_result;
                  abs_result += loop_abs_result;
               }
               else
               {
                  result += loop_result;
                  abs_result += loop_abs_result;
               }
               //
               // Reset k to the largest k we reached
               //
               k = next_backstop;
           }
        }

        return std::make_pair(result, abs_result);
     }

     struct iteration_terminator
     {
        iteration_terminator(std::uintmax_t i) : m(i) {}

        bool operator()(std::uintmax_t v) const { return v >= m; }

        std::uintmax_t m;
     };

     template <class Seq, class Real, class Policy>
     Real hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, long long& log_scale)
     {
        BOOST_MATH_STD_USING
        iteration_terminator term(boost::math::policies::get_max_series_iterations<Policy>());
        std::pair<Real, Real> result = hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, term, log_scale);
        //
        // Check to see how many digits we've lost, if it's more than half, raise an evaluation error -
        // this is an entirely arbitrary cut off, but not unreasonable.
        //
        if (result.second * sqrt(boost::math::policies::get_epsilon<Real, Policy>()) > abs(result.first))
        {
           return boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that fewer than half the bits in the result are correct, last result was %1%", Real(result.first * exp(Real(log_scale))), pol);
        }
        return result.first;
     }

     template <class Real, class Policy>
     inline Real hypergeometric_1F1_checked_series_impl(const Real& a, const Real& b, const Real& z, const Policy& pol, long long& log_scale)
     {
        std::array<Real, 1> aj = { a };
        std::array<Real, 1> bj = { b };
        return hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, log_scale);
     }

  } } } // namespaces

#endif // BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_