Рассмотрим такую задачу: даны три целых числа $$$x$$$, $$$y$$$ и $$$m$$$, $$$0 \leqslant x,y < m < 2^{32}$$$, найти $$$xy\bmod m$$$. По-хорошему хотелось бы просто перемножить эти два числа, а потом применить операцию остатка: ~~~~~ uint32_t prod(const uint32_t x, const uint32_t y, const uint32_t m) { return x * y % m; } ~~~~~ Как вы, возможно, догадываетесь, это решение неверное. Всё дело в том, что в такой процедуре возможно переполнение: операция x * y выполняется в типе uint32_t, и на самом деле промежуточным результатом выполнения этой операции будет не $$$xy$$$, а $$$xy\bmod2^{32}$$$. Если после этого взять результат по модулю $$$m$$$, он может отличаться от правильного:
Выход прост — перемножать необходимо в большем типе: ~~~~~ uint64_t prod_uint64(const uint64_t x, const uint64_t y, const uint64_t m) { return x * y % m; } ~~~~~ Если так делать, то, поскольку $$$xy<2^{64}$$$, это произведение точно не переполнится, и после взятия результата по модулю получится правильный ответ.
Вопрос: а что делать, если $$$x$$$, $$$y$$$ и $$$m$$$ могут быть больше? Предлагаю следующее.
- Бинарное умножение. Так же, как и бинарное возведение в степень, существует бинарное умножение: чтобы посчитать $$$xy$$$, посчитаем $$$x\left\lfloor\frac y2\right\rfloor$$$, сложим это число само с собой и, возможно, добавим ещё $$$x$$$. Это потратит $$$\mathcal O(\log y)$$$ действий, но среди них не будет ничего, кроме сложения и вычитания! ~~~~~ uint64_t sum(const uint64_t x, const uint64_t y, const uint64_t m) { uint64_t ans = x + y; if (ans < x || ans >= m) ans -= m; return ans; } uint64_t prod_binary(const uint64_t x, const uint64_t y, const uint64_t m) { if (y <= 1) return y ? x : 0; uint64_t ans = prod_binary(x, y >> 1, m); ans = sum(ans, ans, m); if (y & 1) ans = sum(ans, x, m); return ans; } ~~~~~
- Умножение с помощью
int128. Чтобы перемножить два 32-битных числа, нужна 64-битная промежуточная переменная. А чтобы перемножить два 64-битных числа, нужна 128-битная переменная! В современных 64-битных компиляторах C++ (кроме разве что Microsoft® Visual C++®) есть специальный тип__int128, позволяющий осуществлять операции над 128-битными числами. ~~~~~ int64_t prod_uint128(const uint64_t x, const uint64_t y, const uint64_t m) { return (unsigned __int128)x * y % m; } ~~~~~ - Умножение с помощью вещественного типа. Что такое $$$xy\bmod m$$$? Это на самом деле $$$xy-cm$$$, где $$$c=\left\lfloor\frac{xy}m\right\rfloor$$$. Давайте тогда попробуем посчитать $$$c$$$, а отсюда найдём $$$xy\bmod m$$$. При этом заметим, что нам не требуется находить $$$c$$$ прямо точно. Что будет, если мы случайно посчитаем, скажем, $$$c-4$$$? Тогда при подсчёте остатка мы вычислим $$$xy-(c-4)m=xy-cm+4m=xy\bmod m+4m$$$. На первый взгляд, это не то, что нам надо. Но если $$$m$$$ не слишком велико и $$$xy\bmod m+4m$$$ не вылезло из 64-битного типа, то после этого можно по-честному взять остаток и получить ответ.
Получается такой код: ~~~~~ uint64_t prod_double(const uint64_t x, const uint64_t y, const uint64_t m) { uint64_t c = (double)x * y / m; int64_t ans = int64_t(x * y — c * m) % int64_t(m); if (ans < 0) ans += m; return ans; } ~~~~~ ~~~~~ uint64_t prod_long_double(const uint64_t x, const uint64_t y, const uint64_t m) { uint64_t c = (long double)x * y / m; int64_t ans = int64_t(x * y — c * m) % int64_t(m); if (ans < 0) ans += m; return ans; } ~~~~~doubleдостаточно точный для этой задачи, если $$$x$$$, $$$y$$$ и $$$m$$$ меньше $$$2^{57}$$$.long doubleже хватает на числа, меньшие $$$2^{63}$$$, однако стоит помнить, чтоlong doubleдля этого должен быть 80-битным, а это верно не на всех компиляторах: например, в Microsoft® Visual C++®long double— то же самое, чтоdouble.
Обратите внимание, что этот способ неприменим, если $$$m>2^{63}$$$: в этом случаеansнельзя хранить вint64_t, так как, возможно, $$$\mathtt{ans}\geqslant 2^{63}$$$ и произойдёт переполнение, из-за которого выполнится ветка(ans < 0)и мы получим неверный ответ.
Видно, что Microsoft® Visual C++® отстаёт от остальных компиляторов в развитии в наличии технических средств для умножения больших чисел по модулю, поэтому, если мы хотим, чтобы функция работала на всех компиляторах быстро, ей необходима какая-то свежая идея. К счастью, именно такую идею в 1960-м году изобрёл Анатолий Карацуба. - Умножение Карацубы. Идея изначально использовалась для быстрого умножения длинных чисел. А именно, пусть $$$x$$$ и $$$y$$$ — два неотрицательных целых числа, меньших $$$N^2$$$. Поделим их с остатком на $$$N$$$: $$$x=Nx_1+x_0$$$, $$$y=Ny_1+y_0$$$. Тогда искомое $$$xy$$$ можно найти как $$$N^2x_1y_1+Nx_1y_0+Nx_0y_1+x_0y_0=N\cdot\bigl(N\cdot x_1y_1+\left(x_0+x_1\right)\left(y_0+y_1\right)-x_1y_1-x_0y_0\bigr)+x_0y_0$$$. Как видно, это преобразование свело умножение $$$x$$$ и $$$y$$$ 1) к $$$\mathcal O(1)$$$ сложениям и вычитаниям чисел, не превосходящих $$$N^4$$$; 2) к трём умножениям чисел, не превосходящих $$$2N$$$ (а именно $$$x_0y_0$$$, $$$x_1y_1$$$ и $$$\left(x_0+x_1\right)\left(y_0+y_1\right)$$$); 3) к двум умножениям на $$$N$$$ чисел, не превосходящих $$$2N^2$$$.
Пункт 1) практически всегда очень прост. В случае с длинными числами пункт 3) также прост: можно взять $$$N$$$, равное степени двойки, и тогда он осуществляется как обычный двоичный сдвиг (оператор<<в C++). Поэтому по существу Карацуба свёл одно умножение чисел, меньших $$$N^2$$$, к трём умножениям чисел, меньших $$$2N$$$. Если эти умножения также свести методом Карацубы, по мастер-теореме асимптотика этого метода составит $$$\Theta\left(\log^{\log_23}N\right)$$$ вместо наивного $$$\Theta\left(\log^2N\right)$$$.
Но нам не потребуется использовать рекурсию, ведь, если длину $$$x$$$ и $$$y$$$ уменьшить вдвое, мы уже сможем воспользоватьсяprod_uint64илиprod_double. Сложность в нашем случае составляет пункт 3): подобрать такое $$$N$$$, чтобы, во-первых, оно было меньше $$$2^{32}$$$ либо в крайнем случае чуть-чуть больше, во-вторых, чтобы на него можно было быстро умножать числа порядка $$$N^2$$$. Оба требования выполнятся, если взять $$$N=\mathrm{round}\left(\sqrt m\right)$$$: действительно, тогда при $$$m<2^{64}$$$ верно $$$N<2^{32}$$$, а $$$\left|m_0\right|=\left|m-N^2\right|\leqslant N<2^{32}$$$; тогда $$$xN=\left(x_1N+x_0\right)N=x_1N^2+x_0N\equiv x_0N-x_1m_0\pmod m$$$, и оба умножения здесь осуществляются над числами порядка $$$N$$$.
Внимательный читатель заметит, что здесь мы обзавелись серьёзной проблемой: извлечение квадратного корня из целого числа. Если вы умеете применять в этой задаче умножение Карацубы, обходя данную проблему (в том числе достаточно быстро находить квадратный корень, написать более быстрый или более короткий код, чем у меня), напишите, пожалуйста, в комментариях!
Поскольку находить произведение $$$\left(x_0+x_1\right)\left(y_0+y_1\right)$$$ оказалось крайне неприятно (напомню,prod_doubleне работает при $$$m>2^{63}$$$), я всё же решил просто вычислить $$$x_0y_1$$$ и $$$x_1y_0$$$ по-честному — таким образом, это не метод Карацубы в истинном смысле, так как я трачу четыре умножения чисел порядка $$$N$$$. ~~~~~ uint64_t dif(const uint64_t x, const uint64_t y, const uint64_t m) { uint64_t ans = x — y; if (ans > x) ans += m; return ans; } bool check_ge_rounded_sqrt(const uint64_t m, const uint64_t r) { return ((r >= 1ull << 32) || r * (r + 1) >= m); } bool check_le_rounded_sqrt(const uint64_t m, const uint64_t r) { return (r == 0 || ((r <= 1ull << 32) && r * (r — 1) < m)); } bool check_rounded_sqrt(const uint64_t m, const uint64_t r) { return check_ge_rounded_sqrt(m, r) && check_le_rounded_sqrt(m, r); } uint64_t rounded_sqrt(const uint64_t m) { uint64_t r = floorl(.5 + sqrtl(m)); if (!check_ge_rounded_sqrt(m, r)) while (!check_ge_rounded_sqrt(m, ++r)); else if (!check_le_rounded_sqrt(m, r)) while (!check_le_rounded_sqrt(m, --r)); return r; } uint64_t prod_karatsuba_aux(const uint64_t x, const uint64_t N, const int64_t m0, const uint64_t m) { uint64_t x1 = x / N; uint64_t x0N = (x — N * x1) * N; if (m0 >= 0) return dif(x0N, x1 * (uint64_t)m0, m); else return sum(x0N, x1 * (uint64_t)-m0, m); } uint64_t prod_karatsuba(const test& t) { uint64_t x = t.x, y = t.y, m = t.modulo; uint64_t N = rounded_sqrt(t.modulo); int64_t m0 = m — N * N; uint64_t x1 = t.x / N; uint64_t x0 = t.x — N * x1; uint64_t y1 = t.y / N; uint64_t y0 = t.y — N * y1; uint64_t x0y0 = sum(x0 * y0, 0, m); uint64_t x0y1 = sum(x0 * y1, 0, m); uint64_t x1y0 = sum(x1 * y0, 0, m); uint64_t x1y1 = sum(x1 * y1, 0, m); return sum(prod_karatsuba_aux(sum(prod_karatsuba_aux(x1y1, N, m0, m), sum(x0y1, x1y0, m), m), N, m0, m), x0y0, m); } ~~~~~
Видно, что на самом деле единственное, что даёт нам тут метод Карацубы — что если вы найдёте большое число $$$N$$$, на которое вы умеете быстро умножать, то вы сможете умножать любые два числа по модулю. Фактически, если бы модуль $$$m$$$ был фиксированным, и было бы много запросов умножения по этому фиксированному модулю, то метод Карацубы был бы молниеносным, поскольку самая затратная операция в нём — это квадратный корень. Хотелось бы, таким образом, взять, например, $$$N=2^{32}$$$ и сделать всё то же самое, что в прошлом пункте, но без квадратного корня. Увы, но я не придумал, как умножать на $$$2^{32}$$$. Можно было бы написать примерно такую функцию: ~~~~~ uint64_t prod_double_small(const uint64_t x, const uint64_t y, const uint64_t m) { uint64_t c = (double)x * y / m; uint64_t ans = (x * y — c * m) % m; return ans; } ~~~~~ Она вычисляет произведение по модулю при условии, что uint64_t c = (double)x * y / m посчиталось абсолютно точно. Но гарантировать, что оно точно посчитается, не представляется возможным, поскольку $$$\frac{xy}m$$$ вполне может оказаться на $$$10^{-18}$$$ меньше, чем какое-то целое число, и типа double не хватит, чтобы это заметить. Именно эту проблему обходит функция prod_karatsuba_aux. Если вы смогли её обойти хитрее, добро пожаловать в комментарии.
Ниже приведены три таблицы по разным компиляторам (да славится имя MikeMirzayanov, ведь именно благодаря нецелевому использованию его Polygon'а я это сделал), в каждой таблице строки соответствуют различным функциям, столбцы соответствуют максимальной разрешённой битности $$$x$$$, $$$y$$$ и $$$m$$$. Если указано CE, значит, с данным компилятором программа не скомпилируется, а если WA — может выдать неверный ответ. В противном случае указано время работы функции на Intel® Core™ i3-8100 CPU @ 3.60GHz. Погрешность приблизительно равна одной-двум наносекундам, на самых медленных функциях может доходить до десяти наносекунд.
Microsoft® Visual C++® 2010
Метод 32 битов 57 битов 63 битов 64 битов prod_uint64 7 ns WA WA WA prod_binary 477 ns 847 ns 889 ns 870 ns prod_uint128 CE CE CE CE prod_double 66 ns 95 ns WA WA prod_long_double 66 ns 98 ns WA WA prod_karatsuba 128 ns 125 ns 138 ns 139 ns GNU G++17
Метод 32 битов 57 битов 63 битов 64 битов prod_uint64 4 ns WA WA WA prod_binary 455 ns 774 ns 841 ns 845 ns prod_uint128 CE CE CE CE prod_double 26 ns 36 ns WA WA prod_long_double 29 ns 20 ns 19 ns WA prod_karatsuba 82 ns 81 ns 91 ns 88 ns GNU G++17 (64 bit)
Метод 32 битов 57 битов 63 битов 64 битов prod_uint64 8 ns WA WA WA prod_binary 313 ns 550 ns 604 ns 630 ns prod_uint128 17 ns 34 ns 30 ns 30 ns prod_double 23 ns 22 ns WA WA prod_long_double 23 ns 24 ns 23 ns WA prod_karatsuba 65 ns 65 ns 69 ns 66 ns
Поэтому базовый рецепт такой: если доступен unsigned __int128, то использовать его, если доступен 80-битный long double, то его тоже должно всегда хватать, а в противном случае надо, если хватает double, использовать double, иначе применить метод Карацубы.
При желании можете попробовать применить эти идеи на задачах специального контеста.
