6.6. Быстрое возведение в степень и умножение

Пусть нам надо число \(a\) возвести в степень \(p\) или умножить на число \(p\) (примеры кода для умножения и возведения буду приводить параллельно; число \(p\) должно быть натуральным, а \(a\) может быть произвольным). Для задачи возведения в степень будем считать, что мы умеем только умножать, для задачи умножения будем считать, что мы умеет только складывать (на самом деле задачу про умножение я привожу тут только для понятности — мне кажется, её понять может быть проще. В реальности я не представляю случая, когда вам пришлось бы так умножать). Основной текст будет про возведение в степень, только примеры буду приводить параллельно для обеих задач.

Конечно, можно сделать это тупо, написав что-нибудь вроде:

// возведение в степень
ans:=1; {в этой переменной будем хранить ответ}
for i:=1 to p do ans:=ans*a;
// умножение
ans:=0; {в этой переменной будем хранить ответ}
for i:=1 to p do ans:=ans+a;

К сожалению, этот цикл выполняется за \(O(p)\). Если \(p\) маленькое, то можно (и нужно!) не мучиться и написать этот простой цикл (ведь любое усложнение программы повышает вероятность ошибки). Если же \(p\) большое, или же \(p\) «среднее», но возводить в степень придётся много раз, то лучше делать все за \(O(\log p)\) (ясно, что, если \(p\) совсем маленькое, то в любом случае надо писать по-тупому). Этот метод на самом деле интуитивно понятен. Если вам нужно посчитать \(6^8\), то можно просто \(6\) возвести в квадрат, потом ещё раз в квадрат, потом ещё раз в квадрат. Эту идею мы будем использовать в нашем алгоритме.

Пусть нам надо посчитать \(a^{43}\). Разложим \(43\) по степеням двойки: \(43 = 32 + 8 + 2 + 1\). Тогда \(a^{43} = a^1 \cdot a^2 \cdot a^8 \cdot a^{32}\). Будем в переменной \(t\) считать \(a^1\), потом \(a^2\), \(a^4\), \(a^8\), \(a^{16}\) и т.д. (просто последовательно возводя \(t\) в квадрат) и постепенно умножать на некоторые из них (на те, которые нам нужны) текущий результат. Как проверить, какие степени нам нужны? Очень просто, побитовыми операциями. Степень \(a^{2^i}\) нам нужна тогда и только тогда, когда \(p \text{ and } 2^i \neq 0\).

Аналогично для умножения: например, \(43a=32a+8a+2a+1a\), в переменной \(t\) будем последовательно считать \(1a\), \(2a\), \(4a\), \(8a\) и т.д. путём последовательного сложения \(t\) с собой; необходимые числа будем добавлять к ответу.

// возведение в степень
ans:=1;
t:=a;
cp:=1; {текущая степень, т.е. t=a^(2^cp)}
while 1 shl cp<p do begin
  if p and (1 shl cp)<>0 then ans:=ans*t;
  t:=t*t;
  inc(cp);
end;
// умножение
ans:=0;
t:=a;
cp:=1; {t=a*(2^cp)}
while 1 shl cp<p do begin
  if p and (1 shl cp)<>0 then ans:=ans+t;
  t:=t+t;
  inc(cp);
end;

(Немного про условие \(1 \text{ shl } cp<p\): ясно, что нам надо действовать до максимальной степени двойки, входящей в \(p\), а это как раз и даёт такое условие)

Можно фактически эту же идею реализовать и по-другому. Посмотрим на последнюю цифру числа \(p\) в двоичной записи (т.е. на остаток от деления \(p\) на 2), и, если это 1, то умножим текущий результат на \(a\). Отбросим последнюю цифру \(p\) (просто поделим на 2, сохранив только целую часть). Посмотрим опять на последнюю цифру (она раньше была второй). Если это 1, то нужно умножить на \(a^2\). Опять отбросим последнюю цифру и посмотрим на ту, что теперь стала последней. Если она 1, то ответ надо умножить на \(a^4\) и т.д., пока число \(p\) не кончится. В отличие от предыдущего способа, мы не будем считать текущую степень отдельно, а будем менять \(p\), постепенно деля его пополам, так, чтобы нужная нам цифра всякий раз оказывалась последней. Вот код:

// возведение в степень
ans:=1;
t:=a;
while p>0 do begin
  if p mod 2 = 1 then ans:=ans*t;
  p:=p div 2;
  t:=t*t;
end;
// умножение
ans:=0;
t:=a;
while p>0 do begin
  if p mod 2 = 1 then ans:=ans+t;
  p:=p div 2;
  t:=t+t;
end;

Код совсем несложный, главное — чётко понимать, что вы делаете и что тут происходит :) Понять корректность этого кода можно и следующим образом. Заметим, что в каждый момент у нас будет верно, что \((ans\cdot t^p)\) равно требуемому ответу (так называемый «инвариант цикла»). Действительно, изначально \(ans=1\), \(t=a\), а \(p\) ещё не менялось, поэтому \(ans\cdot t^p=t^p\) как раз и есть то, что мы и хотим получить. Далее, на очередном шаге есть два варианта: если \(p\) чётное, то \(ans\cdot t^p=ans\cdot (t^2)^{(p/2)}\), т.е. можно просто поделить \(p\) на 2 и возвести \(t\) в квадрат, при этом инвариант не нарушится. Если же \(p\) нечётное, то «отщепим» от \(p\) единицу: \(ans\cdot t^p=(ans\cdot t)\cdot (t^2)^{(p/2)}\) (деление пополам тут имеется ввиду в смысле div 2, т.е. с сохранением только целой части). В конце концов, когда \(p\) станет равно 0, получится, что искомый ответ равен \(ans\cdot t^p=ans\cdot t^0=ans\), т.е. он уже лежит в переменной \(ans\). Фактически, мы как бы «перелили» множители из \(t^p\) в \(ans\).

Нетрудно видеть, что цикл работает за \(O(\log p)\), т.к. каждый раз число \(p\) уменьшается хотя бы в два раза (ровно в два раза, когда на конце нолик, и чуть больше, чем в два раза, когда единичка).

Дальнейшие рассуждения имеют смысл в первую очередь именно для возведения в степень, а не умножения. Как уже было сказано выше, данный код имеет смысл писать, когда число \(p\) большое. Но тогда \(a^p\), скорее всего, не влезет ни в какую память, даже в extended, int64 и т.д. Конечно, может стоять задача вычисления степени длинной арифметикой, но тут — внимание! — приведённый алгоритм не столь эффективен, т.к. требует умножения длинного на длинное, в то время как тупой алгоритм требует лишь умножения длинного на короткое (можете попробовать оценить сложности обоих алгоритмов, учитывая, что в ответе будет \(O(p)\) цифр; у меня получилась сложность \(O(p^2)\) для тупого алгоритма и \(O(p^2\log p)\) для «продвинутого», т.е. «продвинутый» в этом случае — при использовании длинной арифметики — даже хуже простого).

Но нередко бывает нужно вычислить ответ только по некоторому (большому, но лезущему в longint и т.п.) модулю \(inf\), т.е. вычислить \(a^p \bmod inf\) (от слова infinity — бесконечность:) ). Это, пожалуй, и есть как раз самый основной случай применения этого алгоритма.

ans:=1;
t:=a;
while p>0 do begin
  if p mod 2 = 1 then ans:=(ans*t) mod inf;
  p:=p div 2;
  t:=(t*t) mod inf;
end;

Тут, конечно, надо быть осторожным: если хранить все переменные в лонгинте, а \(inf\sim 10^9\) (т.е. порядка \(10^9\)), то во время умножения \(ans \cdot t\) или \(t \cdot t\) может произойти переполнение; надо подумать, что с ним делать.