#ifndef PI_SPIGOT_2019_05_12_H
  #define PI_SPIGOT_2019_05_12_H

  #include <algorithm>
  #include <cstdint>
  #include <iterator>
  #include <type_traits>
  #include <vector>

  #include <ALCryptoLib/cfg/al_crypto_cfg.h>

  #if(__cplusplus >= 201703L)
  namespace al::crypto::math::constants::pi {
  #else
  namespace al { namespace crypto { namespace math { namespace constants { namespace pi {
  #endif

  namespace detail {

  template<const std::uint32_t n>
  constexpr auto pow10() -> std::uint32_t
  {
    static_assert(n <= 9U, "Error: n exceeds its limit of 9 or less.");

    return pow10<n - 1U>() * UINT32_C(10);
  }

  template<>
  constexpr auto pow10<UINT32_C(0)>() -> std::uint32_t
  {
    return UINT32_C(1);
  }

  } // namespace detail

  template<const std::uint32_t result_digit,
           const std::uint32_t loop_digit = UINT32_C(9),
           typename allocator_type = std::allocator<void*>>
  struct pi_spigot
  {
  public:
    using output_value_type = std::uint8_t;
    using input_value_type  = std::uint32_t;

  private:
    // Support up to one million one thousand and one decimal digits.
    static_assert(result_digit <= 1001001UL,
                  "Error: result_digit exceeds its limit of 1,001,001");

    static_assert((loop_digit >= 4U) && (loop_digit <= 9U),
                  "Error: loop_digit is not within its limits of 4...9");

    using output_allocator_type = typename std::allocator_traits<allocator_type>::template rebind_alloc<output_value_type>;
    using input_allocator_type  = typename std::allocator_traits<allocator_type>::template rebind_alloc<input_value_type>;

  public:
    using output_container_type = std::vector<output_value_type, output_allocator_type>;
    using input_container_type  = std::vector<input_value_type, input_allocator_type>;

    using input_iterator_type  = typename input_container_type::iterator;
    using output_iterator_type = typename output_container_type::iterator;

    static constexpr auto input_terms() -> std::uint_fast16_t
    {
      return
        static_cast<std::uint_fast16_t>
        (
          (static_cast<std::uint_fast16_t>(10U * loop_digit) / 3U) + 1U
        );
    }

    static constexpr auto input_scale(std::uint32_t x) -> std::uint32_t
    {
      return static_cast<std::uint32_t>(x * input_terms()) / loop_digit;
    }

    static constexpr auto input_size() -> std::uint32_t
    {
      return input_scale(result_digit);
    }

    static constexpr auto output_size() -> std::uint32_t
    {
      return result_digit;
    }

    static AL_CRYPTO_CFG_INLINE void calculate(input_iterator_type  input_first,
                                               input_iterator_type  input_last,
                                               output_iterator_type output_first,
                                               std::uintmax_t*      operation_count_ptr = nullptr);
  };


  template<const std::uint32_t result_digit,
           const std::uint32_t loop_digit,
           typename allocator_type>
  AL_CRYPTO_CFG_INLINE void pi_spigot<result_digit,
                                      loop_digit,
                                      allocator_type>::calculate
    (input_iterator_type  input_first,
     input_iterator_type  input_last,
     output_iterator_type output_first,
     std::uintmax_t* operation_count_ptr)
  {
    // Use pi_spigot::calculate() to calculate
    // result_digit decimal digits of pi.

    // The caller is responsible for providing both
    // input memory for the internal calculation details
    // as well as output memory for the result of pi.

    // Table of calculation characteristics
    //   digits(*)   operation_count   time[s]
    //  ------------------------------------
    // 1,000,001(9)  191,360,030,868  2,100
    //   500,001(9)   47,840,521,610    520
    //   250,001(9)   11,960,387,349    130
    //   100,001(9)    1,913,780,868     21
    //    50,001(9)      478,496,610      5.3
    //    25,001(9)      119,649,849      1.4
    //    10,001(9)       19,155,868      0.23
    //     5,001(9)        4,794,110      0.055
    //     1,001(9)          193,368      0.002

    //    50,001(8)      527,446,878
    //    25,001(8)      131,887,503
    //    10,001(8)       21,114,378

    //    54,932(4)    1,320,263,154
    //    50,001(4)    1,093,875,003
    //    25,001(4)      273,500,003
    //    10,001(4)       43,775,003
    //     5,001(4)       10,950,003
    //     2,001(4)        1,755,003
    //     1,001(4)          440,003
    //       501(4)          110,628
    //       201(4)           18,003
    //       101(4)            4,628
    //        51(4)            1,222
    //        21(4)              228
    //        11(4)               72

    // (*) Here, the number in parentheses such as
    // (9), (8) or (4) means calculating groups
    // of 9, 8 or 4 digits per loop, corresponding
    // to the template parameter loop_digit.

    // Initialize the input.
    std::fill(input_first,
              input_last,
              static_cast<input_value_type>(detail::pow10<loop_digit>() / 5U));

    std::uint32_t c = UINT32_C(0);

    std::uintmax_t operation_count = UINTMAX_C(0);

    for(auto j  = static_cast<std::uint32_t>(0U);
             j  < result_digit;
             j += loop_digit)
    {
      std::uint64_t d = UINT64_C(0);

      auto i =
        static_cast<std::int32_t>
        (
          input_scale(result_digit - j) - 1U
        );

      for( ; i >= 0; --i)
      {
        d +=
          std::uint64_t(input_first[static_cast<std::uint32_t>(i)])
        * detail::pow10<loop_digit>();

        const auto b =
          static_cast<std::uint32_t>(2U * static_cast<std::uint32_t>(i)) + 1U;

        input_first[static_cast<std::uint32_t>(i)] =
          static_cast<std::uint32_t>(d % b);

        d /= b;

        if(i > 1)
        {
          d *= static_cast<std::uint32_t>(i);
        }

        ++operation_count;
      }

      // Parse groups of pi, where each group has
      // loop_digit digits. If loop_digit is 4,
      // for instance, then successive groups
      // of digits have a form like
      // 3141, 5926, ..., etc.

      const auto digits =
        static_cast<std::uint32_t>
        (
            c
          + static_cast<std::uint32_t>(d / detail::pow10<loop_digit>())
        );

      c = static_cast<std::uint32_t>(d % detail::pow10<loop_digit>());

      const auto n =
        (std::min)
        (
          loop_digit,
          static_cast<std::uint32_t>(result_digit - j)
        );

      auto scale10 = detail::pow10<loop_digit>();

      for(auto   index_scale = static_cast<std::size_t>(0U);
                 index_scale < static_cast<std::size_t>(n);
               ++index_scale)
      {
        scale10 /= UINT8_C(10);

        output_first[static_cast<std::uint32_t>(j + index_scale)] =
          static_cast<output_value_type>
          (
            static_cast<std::uint32_t>(digits / scale10) % UINT8_C(10)
          );
      }
    }

    // Retrieve the operation count if requested.
    if(operation_count_ptr != nullptr)
    {
      *operation_count_ptr = operation_count;
    }
  }

  #if(__cplusplus >= 201703L)
  } // namespace al::crypto::math::constants::pi
  #else
  } // namespace pi
  } // namespace constants
  } // namespace math
  } // namespace crypto
  } // namespace al
  #endif

#endif // PI_SPIGOT_2019_05_12_H