Как ускорить алгоритм поиска числа каталана?

Вы написали нормальный код, который у меня работает за 0.5с.

$ g++ -std=c++11 -pedantic -Wall -Wextra -Werror -O2 original.cpp 
$ time echo 10000 1000000000 | ./a.out 
500322272

real  0m0.477s
user  0m0.476s
sys   0m0.000s

Сложность вашего варианта O(n^2). Если считать через биномиальные коэффициенты рекуррентно, то будет такая же сложность.

Сперва мне показалось что биномиальные коэффициенты можно вычислить через факториалы за линейное время. Но там деление, которое не уложить в модулярную арифметику, а без модулей всё сильно замедляется и толку от линейного числа умножений нет, если каждое умножение требует более чем линейного времени.

В конце концов я вернулся к вашему варианту, убрал переполнения, изменил порядок вычисления чисел, сократил внутренний цикл в два раза. Время снизилось в полтора раза. Этого хватит?

// g++ -std=c++11 -pedantic -Wall -Wextra -Werror -O2 catalan.cpp

#include <iostream>
#include <vector>

int main()
{
    int n;
    uint64_t m;
    if (!(std::cin >> n >> m)) {
        return 1;
    }
    std::vector<uint64_t> catalan(n + 1);
    catalan[0] = 1;
    for (int k = 1; k <= n; ++k) {
        uint64_t s = 0;
        if ((k - 1) % 2 == 0) {
            const int i = (k - 1) / 2;
            s = catalan[i] * catalan[i] % m;
        }
        for (int i = 0, j = k - 1; i < j; ++i, --j)
        {
            s = (s + 2 * catalan[i] * catalan[j]) % m;
        }
        catalan[k] = s;
    }
    std::cout << catalan[n] << "\n";
}

$ g++ -std=c++11 -pedantic -Wall -Wextra -Werror -O2 catalan.cpp 
$ time echo 10000 1000000000 | ./a.out 
287516640

real  0m0.289s
user  0m0.288s
sys   0m0.000s

Stanislav Volodarskiy немного слукавил, когда сказал, что делить по модулю нельзя. На самом деле деление по модулю вполне определено для делителей взаимно простых с самим модулем как умножение на обратный к делителю элемент в кольце.

Для того чтобы разделить N на 7 по модулю 20, если N % 20 = 5 (например, N = 105), нужно найти 7^(-1) = 3 (mod 20) и просто вычислить ((N % mod) * D^(-1)) % mod, что при подстановке тождественно (5 * 3) % 20 = 15. Проверим результат реальным делением 105 на 7: 105 / 7 = 15, 15 % 20 = 15. Таким образом, мы действительно "разделили" остаток R на некоторое число D и получили новый остаток R', причем это произошло корректно для каждого числа в кольце, соответствующего данному остатку R, которое делится на D.

При этом, можно утверждать, что делить на не взаимно простые делители нельзя, так как результат либо не будет существовать, либо будет не единственным. Доказывать это, как и свой прошлый абзац я не буду, послушайте какой-нибудь нормальный курс модульной арифметики, если хотите подробностей.

Таким образом, мы можем делить заданное остатком число на любое взаимно простое с модулем. Чтобы делить заданное остатком число вообще на любое, просто будем хранить степени простых, входящих в разложение модуля, вдобавок к его остатку (точнее к особому остатку).

Если конкретно, то давайте представим наше число N как (A^a * B^b * ... * D^d) * R, где A, B, ..., D - простые, входящие в разложение модуля, а R - произведение всех остальных степеней простых. Тогда остаток от N мы можем вычислить просто перемножив остатки степеней A^a, B^b, ... D^d между собой и с остатком от R. Получается, что поддерживая степени a, b, ..., d и остаток от произведения R, мы можем в любой момент получить остаток от текущего значения N. Очевидно, что можно поддерживать степени при умножении и делении и R при умножении. Но так как R мы будем делить только на некоторое D, в разложение которого не входят простые из разложения M, то D - взаимно просто с модулем, а значит, для деления мы можем просто умножить R на обратное к D, что не составляет никаких проблем.

То есть число 105 по модулю 20 мы представили бы как [{2: 0, 5: 1}, (105 / (2^0 * 5^1)) % 20]. Таким образом, при делении на 15, мы можем сначала разделить наше число N = 105 разделить на 5, которое входит в разложение модуля, а потом разделить остаток на 3, которое уже взаимно просто с модулем. Тогда мы получаем число [{2: 0, 5: 1 - 1}, (1 * 3^(-1)) % 20] = 7 (mod 20).

Если использовать данный способ, то в общем случае асимптотики операций при их оптимальной реализации такие:

• Init(N, M): Инициализация объекта числом N по модулю M: O(sqrt(M) + Mult(N)).
• Mult(N): Умножение рассматриваемого числа на N по рассматриваемому модулю M: O(log(N) + log(M))
• Div(N): Деление рассматриваемого числа на N по рассматриваемому модулю M: O(log(N) + log(M))
• Get(): Получение остатка рассматриваемого числа по рассматриваемому модулю M после K операций умножения: O(log(M) * log(K)) или около того

Так как существует формула C(n) = C(n - 1) * (4n - 2) / (n + 1), то общая сложность алгоритма составит O(N * (log(N) + log(M)) + sqrt(M)), если ничего не путаю

Код:

#include "iostream"
#include "vector"
#include "algorithm"

#define ll long long

using namespace std;

struct Number
{
    int M;
    vector<int> M_factors;
    vector<int> M_degrees;
    vector<int> X_degrees;
    ll remainder;
    int phi;

    void factor(int n, vector<int>& factors, vector<int>& degrees)
    {
        int x = n;
        for (int i = 2; i * i <= n; i++)
        {
            if (x % i == 0)
            {
                factors.push_back(i);
                degrees.push_back(0);

                while (x % i == 0)
                {
                    degrees[degrees.size() - 1]++;
                    x /= i;
                }
            }
        }

        if (x > 1)
        {
            factors.push_back(x);
            degrees.push_back(1);
        }
    }

    int phi_function(vector<int>& factors, vector<int>& degrees)
    {
        int res = 1;
        for (int i = 0; i < factors.size(); i++)
        {
            int z = 1;
            for (int j = 0; j < degrees[i]; j++)
                z *= factors[i];
            res *= z - z / factors[i];
        }
        return res;
    }

    int reverse_element(ll x)
    {
        ll res = 1, pow = phi - 1;

        while (pow) //101
        {
            if (pow & 1)
                res = (res * x) % M;
            x = (x * x) % M;
            pow >>= 1;
        }

        return res;
    }

    Number(int mod)
    {
        M = mod;
        factor(mod, M_factors, M_degrees);
        X_degrees = vector<int>(M_factors.size(), 0);
        remainder = 1;
        phi = phi_function(M_factors, M_degrees);
    }

    Number operator *=(int x)
    {
        for (int i = 0; i < M_factors.size(); i++)
        {
            while (x % M_factors[i] == 0)
            {
                X_degrees[i]++; 
                x /= M_factors[i];
            }
        }
            
        remainder = (remainder * x) % M;

        return *this;
    }

    Number operator /=(int x)
    {
        for (int i = 0; i < M_factors.size(); i++)
        {
            while (x % M_factors[i] == 0)
            {
                X_degrees[i]--;
                x /= M_factors[i];
            }
        }

        remainder = (remainder * reverse_element(x)) % M;

        return *this;
    }

    // Как оказалось, в итоге я не могу оптимизировать get() так, как изначально предполагал. 
    // В моей реализации его асимптотика составляет O(N*log(N)), где N - количество выполненных умножений
    // Однако используя бинарное возведение в степень можно улучшить до O(log(M) * log(N)), где N - количество выполненных умножений 
    int get() 
    {
        ll res = remainder;
        for (int i = 0; i < M_factors.size(); i++)
            for (int j = 0; j < X_degrees[i]; j++)
                res = (res * M_factors[i]) % M;
        return res;
    }
};

int main() 
{
    int n, m;
    cin >> n >> m;

    Number x(m);

    for (int i = 1; i <= n; i++)
    {
        x *= 4 * i - 2;
        x /= i + 1;
    }

    cout << x.get() << '\n';

    return 0;
}

Спасибо Станиславу за проверку и измерение скорости работы программы: n=10000, m=1000000000 -> 0.010s

Считать будем по формуле C(0) = 1, C(n) = (2(2n-1) / (n+1)) * C(n-1). Все множители из формулы разложим на простые (решето Эратосфена) и запомним их степени. При умножении чисел степени их простых складываются, при делении вычитаются.

После всех умножений и делений должны остаться только неотрицательные степени. Простые возводим в эти степени и перемножаем по модулю. Это первый раз где модуль используется.

Сложность оценить точно не берусь. Надеюсь что где-то в районе O(nlogn):

// g++ -std=c++11 -pedantic -Wall -Wextra -Werror -O2 catalan-primes.cpp

#include <algorithm>
#include <cassert>
#include <iostream>
#include <vector>

std::vector<int> make_sieve(int n) {
    std::vector<int> sieve(n);
    for (int i = 0; i < n; ++i) {
        sieve[i] = i;
    }
    for (int i = 2; i * i < n; ++i) {
        if (sieve[i] == i) { // is prime
            for (int j = i * i; j < n; j += i) {
                sieve[j] = i;
            }
        }
    }
    return sieve;
}

void add(const std::vector<int> &sieve, int n, std::vector<int> &powers) {
    while (n > 1) {
        const int p = sieve[n];
        ++powers[p];
        n /= p;
    }
}

void negate(std::vector<int> &powers) {
    for (int &v : powers) {
        v = -v;
    }
}

int main() {
    unsigned n;
    uint64_t m;
    while (std::cin >> n >> m) {
        std::vector<int> sieve = make_sieve(std::max(3u, 2 * n));

        std::vector<int> powers(2 * n, 0);

        // denominators
        for (unsigned i = 1; i <= n; ++i) {
            add(sieve, i + 1, powers);
        }
        negate(powers);

        // numerators
        for (unsigned i = 1; i <= n; ++i) {
            add(sieve, 2, powers);
            add(sieve, 2 * i - 1, powers);
        }

        uint64_t product = 1;

        for (unsigned i = 0; i < 2 * n; ++i) {
            assert(powers[i] >= 0);
            if (powers[i] > 0) {
                for (int j = 0; j < powers[i]; ++j) {
                    product = (product * i) % m;
                }
            }
        }
        std::cout << product << '\n';
    }
}

$ g++ -std=c++11 -pedantic -Wall -Wextra -Werror -O2 catalan-primes.cpp
$ time echo 10000 1000000000 | ./a.out 
287516640

real  0m0.003s
user  0m0.004s
sys   0m0.000s

Среди всей этой математики в ответах не хватает двух алгоритмов, один более известный https://ru.wikipedia.org/wiki/Факториал#Разложение_на_простые_множители, второй менее https://ru.wikipedia.org/wiki/Качающийся_факториал, которые позволяют эффективно следовать совету "Считайте числа Каталана через биномиальные коэффициенты" от @Stanislav Volodarskiy

#include <cassert>
#include <cstdint>
#include <vector>

#if !defined(prime_t)
typedef int32_t prime_t;
#endif
typedef int32_t module_t;

template <typename T, class Pcont, class Fcont>
void factorize(const Pcont &primes, T n,  Fcont &factors)
{
    assert(primes.back() >= n - n/3); // Проверка на дурацкие ошибки

    auto pf = factors.begin();
    for(const auto p : primes) {
        if(n <= 1)
            break;

        T f = 0;

        while(n%p == 0) {
            n /= p;
            f++;
        }

        *pf++ += f;
    }
    assert(n == 1);
}

template <typename T, class Pcont, class Fcont>
void factorial_factorize(const Pcont &primes, T n,  Fcont &factors)
{
    assert(primes.back() >= n - n/3); // Проверка на дурацкие ошибки

    auto pf = factors.begin();
    for(const auto p : primes) {
        if(n < p)
            break;

        T x = n;
        T f = 0;

        do {
            x = x/p;
            f += x;
        } while(x >= p);

        *pf++ += f;
    }
}

// https://ru.wikipedia.org/wiki/Качающийся_факториал
template <typename T, class Pcont, class Fcont>
void swinging_factorial_factorize(const Pcont &primes, T n,  Fcont &factors)
{
    assert(primes.back() >= n - n/3); // Проверка на дурацкие ошибки

    auto pf = factors.begin();
    for(const auto p : primes) {
        if(n < p)
            break;

        T x = n;
        T f = 0;

        do {
            x = x/p;
            f += (x&1);
        } while(x >= p);

        *pf++ += f;
    }
}

template <class Pcont, typename TP, typename TM, typename DT = int64_t>
TM katlan_num(const Pcont &ps, TP n, TM m)
{
        // C(n) = binomial(2n, n) - binomial(2n, n - 1) =
        //      = (2*n)!/( (n!)^2 * (n + 1) )
        //      = swinging_factorial(2*n)/(n + 1)
    std::vector<TP> f(ps.size(), 0);

    #if 0
            // (n!)^2 * (n + 1)
        factorial_factorize(ps, n, f);
        for(auto &&x : f) {
            x = 2*x;
        }
        factorize(ps, n + 1, f);
            // 1/(n!)^2 * (n + 1)
        for(auto &&x : f) {
            x = -x;
        }
            // (2*n)!/( (n!)^2 * (n + 1) )
        factorial_factorize(ps, 2*n, f);
    #else
        factorize(ps, n + 1, f);
        for(auto &&x : f) {
            x = -x;
        }
        swinging_factorial_factorize(ps, 2*n, f);
    #endif

    DT product = 1;
    for(size_t i = 0; i < ps.size(); i++) {
        auto p = ps[i];

        assert(f[i] >= 0);
        for(auto j = f[i]; j > 0; --j) {
            product = (product*p)%m;
        }
    }

    return product;
}

#include <algorithm>
#include <array>
#include <cassert>
#include <initializer_list>
#include <vector>

#if __cplusplus+0 <= 201703 || defined(DONT_CE_PRIMES)
    #define PRIMES_CONSTEXP_
#else
    #define PRIMES_CONSTEXP_ constexpr
    #define CE_PRIMES 1
#endif

template<class Pcont, class Scont>
PRIMES_CONSTEXP_ Pcont
primes_impl_(size_t n, size_t reserve)
{
    Pcont primes;
    primes.reserve(reserve);
    primes.push_back(2);
    primes.push_back(3);

    Scont sieve((n + 1)/2);
    auto pp = sieve.begin() + 1;
    size_t p = 3;                       // 2*(pp - sieve.begin()) + 1;
    size_t jp2 = 4;                     // p*p/2;

    while(jp2 < sieve.size()) {
        do {
            sieve[jp2] = true;
            jp2 += p;
        } while(jp2 < sieve.size());

        pp = std::find(pp + 1, sieve.end(), false);

        p = 2*(pp - sieve.begin()) + 1;
        primes.push_back(p);
        jp2 = p*(p/2);
    }

    while((pp = std::find(pp + 1, sieve.end(), false)) != sieve.end()) {
        primes.push_back(2*(pp - sieve.begin()) + 1);
    }

    return primes;
}

#ifdef CE_PRIMES
    template<typename T, size_t N>
    class ce_vector_array_ : public std::array<T, N>
    {
        size_t m_size;

    public:
        constexpr ce_vector_array_(size_t n = 0) :
            std::array<T, N>{},
            m_size{n}
        { ; }

        constexpr ce_vector_array_(std::initializer_list<T> l) :
            std::array<T, N>{},
            m_size{0}
        { 
            for(auto i : l) {
                (*this)[m_size++] = i;
            }
        }

        constexpr void push_back(const T value) {
            (*this)[m_size++] = value;
        }

        constexpr void reserve(size_t n) {
            assert(n <= (*this).size());
        }

        constexpr size_t size() const {
            return m_size;
        }
    };


    template<int N, typename T=unsigned>
    constexpr auto
    ce_primes(void)
    {
        static_assert(N >= 3, "N >= 3");
        constexpr const auto collected{primes_impl_<
                                        ce_vector_array_<T, N>,
                                        ce_vector_array_<bool, (N + 1)/2>
                                    >(N, 0) };
        std::array<T, collected.size()> primes{};

        for(size_t i = 0; i < collected.size(); i++) {
            primes[i] = collected[i];
        }
        return primes;
    }
#endif

//#include <bit>
#include <cmath>

template<typename T=unsigned>
std::vector<T>
primes(T n)
{
    // Чебышев: 0.92129 <= m <= pi(x)/(x/ln(x)) <= M <= 1.10555
    // Конечно, на бесконечности.
    size_t reserve = 2 + 1.2*n/std::log(n);
                                        // По первым 10⁹, с "запасом" 6..14%
    //size_t reserve = 6 + 1.8*n/std::bit_width((size_t)n); // constexpr, C++20
    //size_t reserve = 6 + 1.8*n/std::log2p1((size_t)n); // n4849, clang
                                        // По первым 10⁹, с "запасом" 3..19%
    auto ps = primes_impl_<std::vector<T>, std::vector<bool> >(n, reserve);

    assert(ps.size() <= reserve);
    assert(n < 200 || reserve <= 1.2*ps.size());

    ps.shrink_to_fit();
    
    return ps;
}

#include <iostream>

int main()
{
    prime_t n;
    module_t m;
    #ifdef USE_CE_PRIMES
        auto ps = ce_primes<32767, prime_t>();
    #endif

    while(std::cin >> n >> m) {
        #ifndef USE_CE_PRIMES
            auto ps = primes(2*n);
        #endif
        std::cout << katlan_num(ps, n, m) << '\n';
    }

    return 0;
}

Асимптотическая сложность выходит n*ln(ln(n)), а практическая вообще неприличная, типа, 51 мкс, если таблицу простых вычислять динамически, каждый раз:

$ clang++-mp-11 -std=c++20 -Wall -g -DNDEBUG -Ofast -march=native t_katlan_Serge3leo.cpp
$ yes '10000 1000000' | head -10000 | (time ./a.out) | sort -u
real    0m0.713s
user    0m0.508s
sys 0m0.011s
516640

И 27 мкс, если простые числа получать во время компиляции (C++20 constexpr):

$ clang++-mp-11 -std=c++20 -Wall -g -DUSE_CE_PRIMES -Dprime_t=int16_t -DNDEBUG -Ofast -march=native t_katlan_Serge3leo.cpp
$ yes '10000 1000000' | head -10000 | (time ./a.out) | sort -u
real    0m0.437s
user    0m0.265s
sys 0m0.010s
516640

P.S. Однако, хотя он и работает гораздо шустрее, но учитывая весьма существенную вероятность того, что katlan_num(ps, n, m) кратен m, возможно это и не самый лучший способ.