//-----------------------------------------------------------------------------
// ackermann.h
//-----------------------------------------------------------------------------

#include <mutex>
#include <stack>
#include <stdexcept>
#include <unordered_map>
#include <vector>

#ifdef __has_include
#if __has_include(<boost/multiprecision/cpp_int.hpp>)
#include <boost/multiprecision/cpp_int.hpp>
#endif // __has_include
#endif // __has_include(<boost/multiprecision/cpp_int.hpp>)

//-----------------------------------------------------------------------------

namespace ackermann
{
    static constexpr auto M = 4, N = 12; // m in [0, 4), n in [0, 12)
    static constexpr auto f_3_11 = 16381;
} // namespace ackermann

//-----------------------------------------------------------------------------

template<typename T> constexpr T ackermann1(T m, T n)
{
    if (m == 0)
        return n + 1;
    if (n == 0)
        return ackermann1(m - 1, T(1));
    return ackermann1(m - 1, ackermann1(m, n - 1));
}

//-----------------------------------------------------------------------------

template<typename T> constexpr T ackermann2(T m, T n)
{
    static std::unordered_map<T, std::unordered_map<T, T>> solution;

    try
    {
        return solution.at(m).at(n);
    }
    catch (std::out_of_range &)
    {
        if (m == 0)
            solution[0][n] = n + 1;
        else if (n == 0)
            solution[m][0] = ackermann2(m - 1, T(1));
        else
            solution[m][n] = ackermann2(m - 1, ackermann2(m, n - 1));
    }

    return solution[m][n];
}

//-----------------------------------------------------------------------------

template<typename T> constexpr T ackermann3(T m, T n)
{
    static std::unordered_map<T, std::unordered_map<T, T>> solution;

    auto m_found = solution.find(m);
    if (m_found != solution.end())
    {
        auto n_found = m_found->second.find(n);
        if (n_found != m_found->second.end())
            return n_found->second;
    }

    if (m == 0)
        solution[m][n] = n + 1;
    else if (n == 0)
        solution[m][n] = ackermann3(m - 1, T(1));
    else
        solution[m][n] = ackermann3(m - 1, ackermann3(m, n - 1));

    return solution[m][n];
}

//-----------------------------------------------------------------------------

template<typename T> constexpr T ackermann4(T m, T n)
{
    static T solution[ackermann::M][1 << (ackermann::M * ackermann::M - 2)];

    if (!solution[m][n])
    {
        if (m == 0)
            solution[0][n] = n + 1;
        else if (n == 0)
            solution[m][0] = ackermann4(m - 1, T(1));
        else
            solution[m][n] = ackermann4(m - 1, ackermann4(m, n - 1));
    }

    return solution[m][n];
}

//-----------------------------------------------------------------------------

template<typename T> constexpr T ackermann5(T m, T n)
{
    static std::once_flag flag;
    std::call_once(
        flag,
        [&] { ackermann4(T(ackermann::M - 1), T(ackermann::N - 1)); });

    return ackermann4(m, n);
}

//-----------------------------------------------------------------------------

template<typename T> constexpr T ackermann6(T m, T n)
{
    std::stack<T, std::vector<T>> stack({m});
    do
    {
        m = stack.top();
        stack.pop();
        if (m == 0)
        {
            ++n;
        }
        else
        {
            stack.push(m - 1);
            if (n == 0)
            {
                ++n;
            }
            else
            {
                stack.push(m);
                --n;
            }
        }
    } while (!stack.empty());
    return n;
}

//-----------------------------------------------------------------------------

#ifdef __has_include
#if __has_include(<boost/multiprecision/cpp_int.hpp>)

using bigint = boost::multiprecision::cpp_int;

bigint ipow(bigint base, bigint exp)
{
    bigint result = 1;
    while (exp)
    {
        if (exp & 1)
            result *= base;
        exp >>= 1;
        base *= base;
    }
    return result;
}

template<typename T> bigint ackermann7(T m, T n)
{
    static bigint (*ack)(unsigned, bigint) = [](unsigned m,
                                                bigint n) -> bigint
    {
        switch (m)
        {
            case 0: return n + 1;
            case 1: return n + 2;
            case 2: return 3 + 2 * n;
            case 3: return 5 + 8 * (ipow(2, n) - 1);
            default:
                return n == 0 ? ack(m - 1, 1) : ack(m - 1, ack(m, n - 1));
        }
    };
    return ack(m, n);
}

//-----------------------------------------------------------------------------

template<typename T> bigint ackermann8(T m, bigint n) // stack overflow
{
    static std::unordered_map<T, std::unordered_map<T, bigint>> solution{
        {0, {{0, 1}, {1, 2}, {2, 3}, {3, 4}, {4, 5}}},
        {1, {{0, 2}, {1, 3}, {2, 4}, {3, 5}, {4, 6}}},
        {2, {{0, 3}, {1, 5}, {2, 7}, {3, 9}, {4, 11}}},
        {3, {{0, 5}, {1, 13}, {2, 29}, {3, 61}, {4, 125}}}};

    T t = static_cast<T>(n);
    auto m_found = solution.find(m);
    if (m_found != solution.end())
    {
        auto n_found = m_found->second.find(t);
        if (n_found != m_found->second.end())
            return n_found->second;
    }

    if (m == 0)
        solution[m][t] = n + 1;
    else if (n == 0)
        solution[m][t] = ackermann8(m - 1, 1);
    else
        solution[m][t] = ackermann8(m - 1, ackermann8(m, n - 1));

    return solution[m][t];
}

#endif // __has_include
#endif // __has_include(<boost/multiprecision/cpp_int.hpp>)

//-----------------------------------------------------------------------------