Бързо степенуване

Fast Exponentiation

В тази тема ще разгледаме няколко задачи, в които се налага да вдигаме число или матрица на голяма степен. Ще разгледаме един алгоритъм (степенуване чрез вдигане на квадрат), който позволява да степенуваме с логаритмична сложност.

Автор: Александър Георгиев

Защо ни е бързо степенуване?

Предполагам всеки стигнал до тук би трябвало да знае най-простия начин за вдигане на число на степен.

int
 pow(int
 num,
 int power) {
    int
 ret =
 1
;
    for
 (int
 i =
 0
;
 i <
 power;
 i+
+)
        ret =
 ret *
 num;
    return
 ret;
}

В някои случаи, обаче, този алгоритъм би могъл да бъде бавен - например когато степента е много голяма. Голямата степен е възможна, тъй като често състезателните задачи изискват отговора по модул. Кога обаче би ни се наложило да вдигаме число на толкова голяма степен, че горният алгоритъм да е бавен?

Деление по модул

В темата за модулна аритметика споменахме, че деление по прост модул изисква вдигане на степен равна на модула минус две. Когато модулът е голям - примерно 1,000,000,007 или 1,000,000,009 - за да разделите на дадено число се налага да го вдигнете на степен, съответно, 1,000,000,005 или 1,000,000,007. Със стандартния алгоритъм това би отнело няколко секунди само за едно единствено деление!

Вдигане на матрица на степен

В темата за намиране на брой пътища с дадена дължина в граф показахме, че това може да стане чрез вдигане на матрицата на инцидентност на графа на степен дължината, която искаме. Нищо не пречи (даже напротив - често се случва) от нас да искат да намерим броя пътища с дължина, да кажем, 1,000,000,000. Това означава, че трябва да вдигнем матрицата на инцидентност на графа на тази степен - отново голямо число. Имайки предвид, че алгоритъмът за умножение на матрици е O(N³), това прави общата сложност за вдигане на матрица на степен O(log(P)∙N³). Тъй като това не е малко, бързото вдигане на степен е жизнено важно при тези задачи.

Резултат по модул

Дори малки числа, вдигнати на относително ниска степен, лесно надхвърлят лимитите на вградените типове. Например, 42¹³ = 1,265,437,718,438,866,624,512, което не се побира дори в 64-битов тип. Затова много често (почти винаги), отговорът се иска по модул някакво число. В горния пример, да кажем, ако изберем модул 1,000,000,007, бихме получили:

42¹³ = 1265437718438866624512 % 1000000007 = 802657452.

В тази тема примерите ще са точно такива - навсякъде ще извършваме операциите по някакъв предварително дефиниран модул MOD. Разбира се,тривиално можете да промените кода да не извършва операциите по модул, като просто го разкарате.

Степенуване чрез вдигане на квадрат

!	Забележете, че когато изчислявате сложността на бързо вдигане на степен, трябва да вземете в предвид и цената за умножение. Например вдигането на `P`-та степен на матрица с размер `M` на `M` би било със сложност `O(log(P)∙M³)`.

Най-популярният (поради простотата си) алгоритъм за бързо степенуване е този чрез вдигане на квадрат. Той се базира на техниката разделяй и владей и както почти всички такива алгоритми има сложност O(log(P)), където P е степента, на която вдигаме числото.

Идеята е следната: нека имаме някакво число N, което искаме да вдигнем на степен P. Ще се възползваме от това, че:

Ако P е четно, то N^P = N^P/2 * N^P/2;
Ако P е нечетно, то N^P = N^P/2 * N^P/2 * N;

Забележете, че N^P/2 се ползва два пъти, но можем да изчислим само веднъж, да запазим в променлива, и да я вдигнем на квадрат (от където идва и името на алгоритъма). Самото N^P/2 можем да изчислим рекурсивно ползвайки същия алгоритъм. Тъй като степента пада на половина на всяка стъпка е що-годе очевидно как стигаме до сложността O(log(P)).

Рекурсивна имплементация

Имплементация на горния алгоритъм ползвайки рекурсия би била следната:

int
 fastPow(int
 num,
 int power) {
    if
 (power =
=
 0)
        return
 1
;
    int
 half =
 fastPow(num,
 power >
>
 1
);
    int
 squared =
 ((long
 long
)half *
 half) %
 MOD;
    return
 (power &
 1
) ?
 ((long
 long
)squared *
 num) %
 MOD :
 squared;
}

Тук базовият ни случай е когато степента е 0: всяко число на степен 0 дава 1. След това изчисляваме half = num^power/2 викайки рекурсията за power / 2. Накрая, ако степента е била нечетна, връщаме half * half * num, а ако е била четна - half * half.

Итеративна имплементация

Аз лично съм по-голям фен на итеративната имплементация, която е не по-сложна, а е малко по-бърза. За съжаление тя е и малко по-трудна за разбиране, съответно ще я срещате доста по-рядко.

Идеята е следната. Нека, например, имаме N = 3 и P = 13 (тоест искаме да сметнем 3¹³). Тринадесет, представено в двоичен запис, е 1101. Което означава, че искаме да сметнем 3⁸ * 3⁴ * 3¹. Забележете, че тъй като имаме вдигане на степен, умножаваме отделните степени на двойката. Ако имахме умножение, бихме ги събирали.

Сега да видим как допринасят за резултата всеки от битовете в двоичния запис на степента:

0001: 3^0001(2) = 3¹ = 3
0010: 3^0010(2) = 3² = 9 = 3 * 3
0100: 3^0100(2) = 3⁴ = 81 = 9 * 9
1000: 3^1000(2) = 3⁸ = 6561 = 81 * 81

Всеки следващ бит (в ред от най-младшите към най-старшите) допринася с числото от предходния бит на квадрат. Така можем да направим цикъл по битовете на числото (вдигайки го на квадрат на всяка итерация) и добавяйки го в произведението, където имаме 1-ца в двоичния запис на степента.

Тъй като степента P = 13, което е 1101 в двоичен запис, резултатът е: 3¹³ = 3 * 81 * 6561 = 1594323.

В код това изглежда по следния начин:

int
 fastPow(int
 num,
 int power) {
    int
 ret =
 1
;
    for
 (int
 i =
 0
;
 i <
 31
;
 i+
+) {
        if
 (power &
 (1
 <
< i))
            ret =
 ((long
 long
)ret *
 num) %
 MOD;
        num =
 ((long
  long
)num *
 num) %
 MOD;
    }
    return
 ret;
}

Леко модифицирана версия на горния код, която спира след намирането на най-старшия сетнат бит (съответно работи малко по-бързо при малки степени), е следната:

int
 fastPow(int
 num,
 int power) {
    int
 ret =
 1
;
    for
 (;
 power >
 0
;
 power >
>
=
 1) {
        if
 (power &
 1)
            ret =
 ((long
 long
)ret *
 num) %
 MOD;
        num =
 ((long
 long
)num *
 num) %
 MOD;
    }
    return
 ret;
}

Бавно умножение

Забележете, че горният алгоритъм работи и за други операции - например умножение! Да кажем 3 * 13 = 3 * 8 + 3 * 4 + 3 * 1 = 39. В този момент би трябвало да се питате защо по дяволите бихме правили това, при положение, че умножението в C++ е O(1). Защо бихме го направили O(log(M))? Та това е бавно умножение!

Това всъщност е хитър трик, с който можем да се справим с модули, надхвърлящи лимитите на 32-битов int - например MOD = 1,000,000,000,000,000,003. Ако трябва да върнем отговор по този модул, например не бихме могли да умножим две числа A и B. Вярно, че ги пазим по този модул (тоест са по-малки от него), но временната стойност на умножението ще изисква около 120 бита - доста повече от 64-те, които позволяват сегашните компютри.

Вместо това можем да ползваме горния алгоритъм, модифициран за умножение, вместо вдигане на степен. При него най-"тежката" операция, която ще имаме, е събиране. Забележете, че сумата на две числа по този модул все още не надхвърля лимитите на 64-битово цяло число - тоест можем да събираме числа по този модул колкото си искаме. Така с повече на брой стъпки, но такива, които можем да направим, можем да направим и умножение по голям модул!

int
 safeMul(int
 a,
 int b) {
    int
 ret =
 0
;
    for
 (;
 b >
 0
;
 b >
>
=
 1) {
        if
 (b &
 1
) ret =
 (ret +
 a) %
 MOD;
        a =
 (a +
 a) %
 MOD;
    }
    return
 ret;
}

Референции

Exponentiation by squaring (en.wikipedia.org)

За да предложите корекция, селектирайте думата или текста, който искате да бъде променен,
натиснете Enter и изпратете Вашето предложение.

Страницата е посетена 11466 пъти.