Недавно я столкнулся с задачей, в которой необходимо было найти $$$N$$$-ное ($$$1 \leq N \leq 10^4$$$) число Каталана по заданному натуральному модулю $$$M$$$ ($$$2 \leq M \leq 10^9$$$).
Учитывая, что вычисление модуля весьма тяжелая операция, несложно понять, что квадратный алгоритм, основанный на главной формуле $$$C_{0}=1$$$ и $$$C_n=\sum_{i=0}^{n-1}C_{i}C_{n-1-i}$$$ для $$$n \geq 1$$$ для решения этой задачи подходит не очень. А с остальными формулами проблема очень проста: в них необходимо вычислять биномиальные коэффициенты по непростому модулю, для чего нужно по этому самому модулю научиться делить.
(На самом деле, как оказалось немного позднее, делить по модулю тут вовсе не обязательно — есть и более хитрые способы, которые позволяют добиться асимптотики $$$O(Nlog(logN))$$$ без какого-либо деления, но не решение задачи является главной темой поста. Если оно вам так же интересно, то держите ссылку: Как ускорить алгоритм поиска числа каталана?).
После недолгого гугления я нашел упоминание о таком делении на e-maxx'e в статье "Длинная арифметика" в блоке "Длинная арифметика в факторизованном виде", однако почти ничего о реализации там сказано не было. Тогда я начал искать информацию о делении по непростому модулю тут, на CF, но снова потерпел неудачу — либо я не умею искать, либо тут действительно нет ни одной статьи как минимум на русском языке, посвященной этому вопросу, из-за чего я и пишу эту заметку.
После некоторого насилия над своим разленившимся мозгом я понял, что для вычисления значения нашего числа по некоторому модулю необходимо и достаточно хранить разложение этого числа на простые, входящие в разложение модуля, а также остаток произведения всех остальных простых.
Другими словами: пусть мы желаем найти остаток числа $$$N$$$ при делении на $$$M = {p_1}^{a_1} * {p_2}^{a_2} * \ldots * {p_k}^{a_k}$$$. Тогда мы представляем наше число в виде $$$N = ({p_1}^{b_1} * {p_2}^{b_2} * \ldots * {p_k}^{b_k}) * ({q_1}^{c_1} * {q_2}^{c_2} * \ldots * {q_m}^{c_m}) = {p_1}^{b_1} * {p_2}^{b_2} * \ldots * {p_k}^{b_k} * R$$$. Тогда понятно, что $$$N \mod M = ({p_1}^{b_1} \mod {M}) * ({p_2}^{b_2} \mod {M}) * \ldots * ({p_k}^{b_k} \mod {M}) * ({R} \mod {M})$$$.
Теперь заметим, что этого достаточно для возможности умножения и деления нашего числа $$$N$$$ по любому натуральному модулю $$$M$$$. Для умножения все очевидно: нужные степени увеличатся, остаток нормально вычислится. А про деление скажу чуть подробнее.
Пусть мы делим наше число $$$N$$$ на некоторое $$$D = G * D'$$$, где $$$G = gcd(D, M)$$$. Мы знаем, что $$$G$$$ также равен $$${p_1}^{d_1} * \ldots * {p_k}^{d_k}$$$, где $$$p_1, \ldots, p_k$$$ — простые в разложении $$$M$$$, и что $$$D'$$$ взаимно просто с модулем. Тогда для деления на $$$D$$$ мы сначала отнимем от степеней простых в $$$N$$$ соответственные значения $$$d_1, \ldots, d_k$$$, а потом просто "разделим" $$$R$$$ в составе $$$N$$$ на $$$D'$$$ (умножим на $$$D'^{-1} \mod M$$$), так как обратный элемент существует для любого числа, взаимно простого с модулем и деление для этого взаимно простого числа работает по непростому модулю так же, как и по простому.
Теперь разберем это на маленьком примере: Число $$$105$$$ по модулю $$$20$$$ мы представили бы как [{2: 0, 5: 1}, (105 / (2^0 * 5^1)) % 20]. Таким образом, при делении на $$$15$$$, мы можем сначала разделить наше число $$$N = 105$$$ разделить на $$$15$$$, которое входит в разложение модуля, а потом разделить остаток на $$$3$$$, которое уже взаимно просто с модулем. Тогда мы получаем число [{2: 0, 5: 1 - 1}, (1 * 3^(-1)) % 20] = [{2: 0, 5: 0}, 7], которое соответствует $$$7 \mod 20$$$. И, действительно, $$$105 / 15 = 7, 7 \mod 20 = 7$$$.
Если использовать данный способ, то в общем случае асимптотики операций при их оптимальной реализации такие:
$$$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()$$$: Получение остатка рассматриваемого числа $$$N$$$ по рассматриваемому модулю $$$M$$$: $$$O(log(M) * log(log(N)))$$$. То есть, если принять сумму степеней в разложении $$$N$$$ за $$$S$$$, то это будет $$$O(log(M) * log(S))$$$. Если мы умножаем только на 32-битные натуральные $$$K$$$ раз, то $$$O(log(M) * log(K))$$$.
Теперь собственно реализация, которую я так жаждал найти. Сильно я ее не оптимизировал, но надеюсь, что сам принцип понятен:
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)
{
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;
}
// В моей реализации его асимптотика составляет O(K + log(M)), где K - количество выполненных умножений
// Однако используя бинарное возведение в степень можно улучшить до O(log(M) * log(K))
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;
}
};
P.S. Сама задача о числах Каталана легко решается и без этой громоздкой структуры с помощью хранения разложения вообще всех множителей/делителей на все простые — даже при плохой реализации получится $$$O(N \sqrt{N})$$$
