\Section{Длинная арифметика}{11 декабря 2017}{Сергей Копелиович}

Мы займёмся целыми беззнаковыми числами.
Целые со знаком -- ещё отдельно хранить знак.
Вещественные -- то же, но ещё хранить экспоненту: $12.345 = 12345e{-3}$, мы храним $12345$ и $-3$.

\down
Удобно хранить число в ``массиве цифр'', младшие цифры в младших ячейках. \\
Во примерах ниже мы выбираем систему счисления $\t{BASE} = 10^k$, $k \to \max \colon$
нет переполнений.\\
Пусть есть длинное число $a$. При оценки времени работы будем использовать обозначения: \\
$|a| = n$ -- битовая длина числа и $\mfrac{n}{k}$ -- длина числа, записанного в системе $10^k$.
Помните, $\max k \approx 9$.

\down
Если мы ленивы и уверены, что в процессе вычислений не появятся числа длиннее $N$, \\
наш выбор -- \t{int[N];}, иначе обычно используют \t{vector<int>} и следят за длиной числа.

\down
Примеры простейших операций:
\begin{code}
const int N = 100, BASE = 1e9, BASE_LEN = 9;
void add( int *a, int *b ) { // `\red{сложение за $\O(n/k)$}`
	for (int i = 0; i + 1 < N; i++) // `+1, чтобы точно не было range check error`
		if ((a[i] += b[i]) >= BASE)
			a[i] -= BASE, a[i + 1]++;
}
int divide( int *a, int k ) { // `\red{деление на короткое за $\O(n/k)$}`, `делим со старших разрядов`
	long long carry = 0; // `перенос с более старшего разряда, он же остаток`
	for (int i = N - 1; i >= 0; i--) {
		carry = carry * BASE + a[i]; // `максимальное значение \t{carry} < \t{BASE}$^2$`
		a[i] = carry / k, carry %= k;
	}
	return carry; // `как раз остаток`
}
int mul_slow( int *a, int *b, int *c ) { // `\red{умножение за $(n/k)^2$}`
	fill(c, c + N, 0);
	for (int i = 0; i < N; i++)
		for (int j = 0; i + j < N; j++)
			`\color{dkred}{c[i + j] += a[i] * b[j]}`; // `здесь почти наверняка произойдёт переполнение`
	for (int i = 0; i + 1 < N; i++) // `\color{dkred}{сначала умножаем, затем делаем переносы}`
		c[i + 1] += c[i] / BASE, c[i] %= BASE;
}
void out( int *a ) { // `\red{вывод числа за $\O(n/k)$}`
	int i = 0;
	while (i && !a[i]) i--;
	printf("%d", a[i--]);
	while (i >= 0) printf("%0*d", BASE_LEN, a[i--]); // `воспользовались таки BASE\_LEN!`
}
\end{code}
Чтобы в строке \t{19} не было переполнения, нужно выбирать \t{BASE} так, что $\t{BASE}^2\,\t{N}$ помещалось
в тип данных. Например, хорошо сочетаются $\t{BASE} = 10^8, \t{N} = 10^3$, тип -- \t{unsigned long long}.

\pagebreak
\vspace*{-1.5em}
\Subsection{Бинарная арифметика}

Пусть у нас реализованы простейшие процедуры: ``\t{+}, \t{-}, \t{*2}, \t{/2}, \t{\scalebox{.9}{\%{}}2}, \t{>}, 
$\begingroup{\displaystyle{\mathrm{\ge}}}\endgroup$, \t{isZero}''.

Давайте выразим через них ``\t{*}, \t{\textbackslash}, \t{gcd}''. Обозначим $|a| = n, |b| = m$.

\down
Умножение будет полностью изоморфно бинарному возведению в степень.
\begin{code}
num mul(num a, num b) {
	if (isZero(b)) return 1; // `если храним число, как \t{vector}, то \t{isZero} за $\O(1)$`
	num res = mul(mul(mul2(a), div2(b));
	if (mod2(b) == 1) add(res, a); // `функция \t{mod2} всегда за $\O(1)$`
	return res;
}
\end{code}

% \pagebreak
% \vspace*{-1.5em}

Глубина рекурсии равна $m$. В процессе появляются числа не более $(n{+}m)$ бит длины $\SO$\\
каждая операция выполняется за $\O(\mfrac{n+m}{k}) \SO$ суммарное время работы $\O((n+m)\mfrac{m}{k})$.

Если большее умножать на меньшее, то $\O(\max(n,m)\min(n,m)/k)$.

%Время работы $\O((n+m)m)$, где $n$ -- битовая длина $a$, $m$ -- битовая длина $b$.

\down
Деление в чём-то похоже... деля $a$ на $b$, мы будем пытаться вычесть из $a$ числа $b, 2b, 4b, \dots$.  
\begin{code}
pair<num, num> div(num a, num b) { // `найдём для удобства и частное, и остаток`
	num c = 1, res = 0;
	while (b < a) // `$(n{-}m)$ раз`
		mul2(b), mul2(c);
	while (!isZero(c)) { // `Этот цикл сделает $\approx n{-}m$ итераций`
		if (a >= b) // `$\O(n)$, так как длины $a$ и $b$ убывают от $n$ до $1$`
			sub(a, b), add(res, c); `$\O(n)$`
		div2(b), div2(c); `$\O(n)$`
	}
	return {res, a};
}
\end{code}
Глубина рекурсии равна $n{-}m$. Все операции за $\O(\mfrac{n}{k}) \SO$ суммарное время $\O((n - m)\mfrac{n}{k})$.
%Время работы $\O((n{-}m{+}1)n)$, где $n$ -- битовая длина $a, \ m$ -- битовая длина $b$.

\down
Наибольший общий делитель сделаем самым простым Евклидом ``с вычитанием''.

Добавим только одну оптимизацию: если числа чётные, надо сразу их делить на два...

\begin{code}
num gcd(num a, num b) {
	int pow2 = 0;
	while (mod2(a) == 0 && mod2(b) == 0)
		div2(a), div2(b), pow2++;
	while (!isZero(b)) { 
		while (mod2(a) == 0) div2(a);
		while (mod2(b) == 0) div2(b);
		if (a < b) swap(a, b);
		a = sub(a, b); // `одно из чисел станет чётным`
	}
	while (pow2--) mul2(a);
	return a;
}
\end{code}
Шагов главного цикла не больше $n{+}m$. Все операции выполняются за $\max(n,m)/k$.\\
Отсюда суммарное время работы: $\O(\max(n,m)^2/k)$.

\Subsection{Деление многочленов за $\O(n \log^2 n)$}

Коэффициенты многочлена $A(x)\colon A[0]$ -- младший, $A[\deg A]$ -- старший. $\gamma(A) = \deg A - 1$.

\down
{\bf Задача:} даны $A(x), B(x) \in \R[x]$, найти $Q(x), R(x) \colon \deg R < \deg B \wedge A(x) = B(x)Q(x) + R(x)$.

\down
Сперва простейшее решение за $\O(\deg A \cdot \deg B)$, призванное побороть страх перед делением:
\begin{code}
pair<F*,F*> divide( int n, F *a, int m, F *b ) { // `$\deg A = n, \ \deg B = m$, $\mathbb{F}$ -- поле`
	F q[n-m+1];
	for (int i = n - m; i >= 0; i--) { // `выводим коэффициенты частного в порядке убывания`
		q[i] = a[i + m] / b[m]; // `\t{m} -- степень $\SO \t{b[m]} \not= 0$`.
		for (int j = 0; j <= m; j--) // `конечно, вычитать имеет смысл, только если $\t{q[i]} \not=0$`
			a[i + j] -= b[j] * q[i]; // `можно соптимизить, перебирать только ненулевые \t{b[j]}`
	}
	return {q, a}; // `в \t{a} как раз остался остаток`
}
\end{code}

Теперь перейдём к решению за $\O(n \log^2 n)$.

Зная $Q$, мы легко найдём $R$, как $A(x) - B(x)Q(x)$ за $\O(n \log n)$. Сосредоточимся на поиске $Q$.

Пусть $\deg A = \deg B = n$, тогда $Q(x) = \mfrac{a_n}{b_n}$. То есть, $Q(x)$ можно найти за $\O(1)$. 

Из этого мы делаем вывод, что $Q$ зависит не обязательно от всех коэффициентов $A$ и $B$.

\begin{Lm}
$\deg A = m, \deg B = n \SO \deg Q = m - n$, и $Q$ зависит только \\
от $m{-}n{+}1$ старших коэффициентов $A$ и $m{-}n{+}1$ коэффициентов $B$.
\end{Lm}

%\vspace*{-1.5em}
\begin{proof}
%\hspace*{0.4em}\\
Рассмотрим деление в столбик, шаг которого: $A\ \t{-=}\ \alpha x^i B$. \ $\alpha = \mfrac{A[i + \deg B]}{B[\deg B]}$.

\vspace*{0.1em}
Поскольку $i + \deg B \ge \deg B = n$, младшие $n$ коэффициентов $A$ не будут использованы.
% $\deg R < n \SO$ у $A$ и $BQ$ должны совпадать $m{-}n{+}1$ старший коэффициент.
% В этом сравнении участвуют только $m{-}n{+}1$  старших коэффициентов $A$.
% При домножении $B$ на $x^{\deg Q}$, сравнятся как раз $m{-}n{+}1$ старших коэффициентов $A$ и $B$.
% При домножении $B$ на меньшие степени $x$, в сравнении будут участвовать лишь какие-то первые из этих $m{-}n{+}1$ коэффициентов.
\end{proof}

{\bf Теперь будем решать задачу:} \\
Даны $A,B \in \R[x] \colon \gamma(A) = \gamma(B) = n$, найти $C \in \R[x] \colon \gamma(C) = n$, \\
что у $A$ и $BC$ совпадает $n$ старших коэффициентов.

% \pagebreak
% \vspace*{-1.5em}
\begin{code}
int* Div( int n, int *A, int *B ) // `n -- степень двойки (для удобства)`
  C = Div(n/2, A + n/2, B + n/2) // `нашли старших n/2 коэффициентов ответа`
  A`\color{black}{'}` = Subtract(n, A, n + n/2 - 1, Multiply(C, B))
  D = Div(n/2, A`\color{black}{'}`, B + n/2) // `сейчас A' состоит из n/2 не нулей и n/2 нулей`
  return concatenate(D, C) // `склеили массивы коэффициентов; вернули массив длины ровно n`
\end{code}

Здесь \t{Subtract} -- хитрая функция. Она знает длины многочленов, которые ей передали, и сдвигает вычитаемый многочлен так, чтобы старшие коэффициенты совместились.

\down
{\bf Время работы:} $T(n) = 2T(n/2) + \O(n \log n) = \O(n \log^2 n)$. Здесь $\O(n \log n)$ -- время умножения.

\Subsection{Деление чисел}

Оптимально использовать метод Ньютона, внутри которого все умножения -- FFT. \\
Тогда мы получим асимптотику $\O(n \log n)$. Об этом можно будет узнать на третьем курсе.

Сегодня лучшими результатами будут $\O((n/k)^2)$ и $\O(n \log^2 n)$.

\down
{\bf Простейшие методы} (оценка времени деление числа битовой длины $2n$ на число длины $n$).

\down
1. Бинпоиск по ответу: $n^3/k^2$ при простейшем умножении, $n^2 \log n$ при Фурье внутри.\\
2. Двоичное деление: $n^2/k$ времени.\\
3. Деление в столбик: $n^2/k^2$ времени. На нём остановимся подробнее.

\pagebreak
\vspace*{-1.5em}
\Subsection{Деление чисел за $\O((n/k)^2)$}

Делить будем в столбик. У нас уже было деление многочленов за квадрат.
Если мы научимся вычислять за $\O(n/k)$ старшую цифру частного, мы сможем воспользоваться им без изменений.
Пусть даны числа $a, b, |a| = n, |b| = m$.

\begin{Lm}
Старшая цифра $\mfrac{a}{b}$ отличается от $x = \mfrac{a_na_{n-1}}{b_mb_{m-1}}$ не более чем на $1$.
\end{Lm}

\vspace*{-1.5em}
\begin{proof}
$
\mfrac{a_na_{n-1}}{b_mb_{m-1}+\sfrac{1}{base}} \le \mfrac{a}{b} \le \mfrac{a_na_{n-1}+\sfrac{1}{base}}{b_mb_{m-1}}
\SO
|\mfrac{a}{b} - x| \le (\mfrac{a_na_{n-1}+\sfrac{1}{base}}{b_mb_{m-1}} {-} x) + (x {-} \mfrac{a_na_{n-1}}{b_mb_{m-1}+\sfrac{1}{base}})
\t{:=}\ y
$

\down
Заметим, $b_m \neq 0 \SO b_mb_{m-1} \ge base$. Продолжаем преобразования:\\
$y \le \mfrac{1}{base}{\cdot}\mfrac{1}{b_mb_{m-1}} + \mfrac{a_na_{n-1}}{base(b_mb_{m-1})^2} = 
\mfrac{1}{base{\cdot}b_mb_{m-1}}(1 + \mfrac{a_na_{n-1}}{b_mb_{m-1}}) \le \mfrac{1}{base^2}(1 + \mfrac{base^2}{base}) \le 1$.
\end{proof}

\up
\THE{Алгоритм деления:}

Длина частного, т.е. $\frac{n{-}m}{k}$, раз вычисляем $\alpha$ -- приближение старшей цифры частного за $\O(1)$,
затем умножением за $\O(\frac{n}{k})$ вычитаем $(\alpha{-}1) b\,10^{ki}$ из $a$ и не более чем двумя вычитаниями
$b\,10^{ki}$ доводим дело до конца. Важно было начать с $\alpha{-}1$, чтобы не уйти в минус при вычитании.