Бързо степенуване
Накратко: Бързото степенуване смята aⁿ с O(log n) умножения, като повдига на квадрат и гледа двоичните цифри на n; с остатък след всяко умножение смята и степен по модул.
Проблем
Ключът на HTTPS връзката
Когато браузърът отваря HTTPS страница, той и сървърът се договарят за общ таен ключ, без да го пращат по мрежата. В класическия вариант (Diffie–Hellman) всеки смята степен по модул: gᵃ mod p, където a е тайно число, дълго десетки цифри. Същата сметка стои в RSA подписите и в hash на текст, който умножава по степени на основата.
Наивният начин е цикъл, който умножава n пъти:
def power_loop(a, n, m): result = 1 for _ in range(n): result = result * a % m return resultПри показател от 77 цифри (например 2²⁵⁵ − 19) това са около 6 · 10⁷⁶ умножения. При милиард умножения в секунда цикълът ще върви 10⁶⁷ секунди, а Вселената е на около 4 · 10¹⁷. Сървърът обаче прави хиляди такива връзки в секунда.
Идеята
Квадратът на половината
3¹³ не иска 13 умножения. 3¹³ = 3⁶ · 3⁶ · 3, а 3⁶ = 3³ · 3³. Всяка стъпка разполовява показателя:
- aⁿ = (aⁿᐟ²)², ако n е четно;
- aⁿ = (aⁿᐟ²)² · a, ако n е нечетно (n / 2 е целочислено деление);
- a⁰ = 1: това е дъното.
def power(a, n, m): if n == 0: return 1 half = power(a, n // 2, m) square = half * half % m if n % 2 == 1: return square * a % m return square
print(power(3, 13, 1_000_000_007))print(power(5, 1_000_000_000, 1_000_000_007))1594323142848001Остатъкът по m след всяко умножение пази числата малки: (x · y) mod m е същото като ((x mod m) · (y mod m)) mod m. Тук m = 1 000 000 007, просто число, честа граница в hash-ове и задачи по програмиране.
Колко струва? Всяко извикване дели n на 2, точно като двоичното търсене, затова дълбочината е около log₂ n: при n = 10⁹ са 30 рамки. Извикването не е опашково (след него има повдигане на квадрат), но от урока „Стекът от извиквания отвътре“ знаеш, че 30 рамки са нищо за стека. Работата е най-много две умножения на ниво: O(log n) вместо O(n).
Идеята
Без рекурсия: двоичните цифри на n
Същото може и с цикъл. 13 в двоично е 1101: 13 = 8 + 4 + 1, затова 3¹³ = 3⁸ · 3⁴ · 3¹. Цикълът чете цифрите на n отдясно наляво (n % 2 е последната, а n /= 2 я маха) и пази a, повдигнато на квадрат: 3¹, 3², 3⁴, 3⁸… При цифра 1 текущото a влиза в резултата.
def power(a, n, m): result = 1 a %= m while n > 0: if n % 2 == 1: result = result * a % m a = a * a % m n //= 2 return result
print(power(3, 13, 1_000_000_007))print(power(5, 1_000_000_000, 1_000_000_007))1594323142848001Умноженията са колкото цифрите на n (по едно повдигане на квадрат на цифра) плюс единиците в него (по едно умножение на резултата). Това е версията, която ползват библиотеките: няма рекурсия и не зависи от размера на стека.
Виж го
Квадратите влизат в резултата
def power(a, n, m): result = 1 a %= m while n > 0: if n % 2 == 1: result = result * a % m a = a * a % m n //= 2 return result
print(power(3, 13, 1_000_000_007)) # 1594323Стъпка по стъпка, с думи
- result започва от 1. n = 13, в двоично 1101: всяка двоична цифра решава дали a влиза в резултата.
- n = 13 е нечетно: result = 1 · 3 = 3. После a става 3² = 9, а n става 6.
- n = 6 е четно: result остава 3. Само a се повдига на квадрат: 9 · 9 = 81 = 3⁴, а n = 3.
- n = 3: result = 3 · 81 = 243 = 3⁵. a = 6561 = 3⁸, n = 1.
- n = 1: result = 243 · 6561 = 1594323 = 3¹³, защото 13 = 8 + 4 + 1. Цикълът пак повдига a на квадрат: 3¹⁶ е умножение, което никой не ползва.
- n = 0 и цикълът спира: 7 умножения вместо 13. При n = 10⁹ са 43 вместо милиард.
Виж го
Броят умножения
Двете версии с брояч, на показатели от 10 до 10 000 000, а бързата и на милиард:
multiplications = 0
def power_loop(a, n, m): global multiplications result = 1 for _ in range(n): multiplications += 1 result = result * a % m return result
def power(a, n, m): global multiplications result = 1 a %= m while n > 0: if n % 2 == 1: multiplications += 1 result = result * a % m multiplications += 1 a = a * a % m n //= 2 return result
m = 1_000_000_007for n in [10, 1_000, 100_000, 10_000_000]: multiplications = 0 power_loop(5, n, m) slow = multiplications multiplications = 0 power(5, n, m) print(f"n = {n}: цикъл {slow}, бързо {multiplications} умножения")multiplications = 0power(5, 1_000_000_000, m)print(f"n = 1000000000: бързо {multiplications} умножения")n = 10: цикъл 10, бързо 6 умноженияn = 1000: цикъл 1000, бързо 16 умноженияn = 100000: цикъл 100000, бързо 23 умноженияn = 10000000: цикъл 10000000, бързо 32 умноженияn = 1000000000: бързо 43 умноженияПусни ги и тук. Bench брои по същото правило (едно умножение по модула е една операция), затова числата са тези от изхода горе. Милиард с наивния цикъл не го пускаме: Bench го оценява по кривата.
- Цикъл: n умножения
O(n) - Повдигане на квадрат
O(log n)
| Алгоритъм | n = 10 | n = 1000 | n = 100 000 | n = 10 000 000 | n = 1 000 000 000 |
|---|---|---|---|---|---|
| Цикъл: n умножения | 10 | 1000 | 100 000 | 10 000 000 | ≈ 1 000 000 000 |
| Повдигане на квадрат | 6 | 16 | 23 | 32 | 43 |
def power(a, n, m): result = 1 a %= m while n > 0: if n % 2 == 1: result = result * a % m a = a * a % m n //= 2 return resultНаивният цикъл расте като n, бързото като броя цифри на n: десет пъти по-голямо n добавя около 3,3 цифри, тоест 3 до 7 умножения.
Пробвай
15 или 16
Кой показател иска повече умножения?
multiplications = 0
def power(a, n, m): global multiplications result = 1 a %= m while n > 0: if n % 2 == 1: multiplications += 1 result = result * a % m multiplications += 1 a = a * a % m n //= 2 return result
for n in [15, 16]: multiplications = 0 power(3, n, 1_000_000_007) print(f"{n}: {multiplications} умножения")15 = 1111: четири цифри и четири единици, 8 умножения. 16 = 10000: пет цифри и една единица, 6 умножения. По-голям показател не значи повече работа: важат цифрите и единиците.
Идеята
Когато произведението не се побира
Остатъкът държи числата под m, но произведението на две такива числа е до m². Ако m² не се побира в типа, сметката греши, и то тихо. Ето същата степен с по-голям модул:
m = 10_000_000_019print("power:", power(5, 1_000_000_000, m))print("pow: ", pow(5, 1_000_000_000, m))power: 475867028pow: 475867028int расте колкото трябва, затова същият код дава верния отговор и при голям модул. Цената е, че всяко умножение на голямо число е по-бавно. Вграденото pow(a, n, m) прави същото повдигане на квадрат, но в C.Идеята
Вграденото и показател от 77 цифри
В истинския код степента по модул не се пише на ръка:
secret = 2**255 - 19p = 2**127 - 1print("Цифри в показателя: " + str(len(str(secret))))print("Умножения при повдигане на квадрат: " + str(secret.bit_length() + secret.bit_count()))print("5^secret mod p =", pow(5, secret, p))Цифри в показателя: 77Умножения при повдигане на квадрат: 5085^secret mod p = 97285997395055909576908201857106207893pow(a, n, m) с три аргумента прави повдигане на квадрат и взима остатъка по пътя. a ** n % m е вярно, но първо строи цялото aⁿ: при показател от 77 цифри никога няма да свърши. 508 умножения на числа от 39 цифри вместо 6 · 10⁷⁶: затова ключът се договаря за милисекунди. Библиотеките броят по-малко и от това, защото гледат по няколко цифри наведнъж, но идеята е същата.
Пробвай
Намери грешката
Тази версия връща верния резултат. Само че при n = 1 000 000 прави 2 048 575 умножения вместо 27. Кой ред е виновен?
Избери реда с грешката.
Половината се смята два пъти. Всяко ниво удвоява извикванията: 2, 4, 8… и на дъното са около n. Дълбочината пак е log₂ n, затова стекът не забелязва нищо, а работата се връща на O(n). Поправката е половината в променлива и после half * half.
Пробвай
Попълни
Кога текущото a влиза в резултата?
def power(a, n, m):
result = 1
a %= m
while n > 0:
if :
result = result * a % m
a = a * a % m
n //= 2
return resultКогато последната двоична цифра на n е 1, тоест n е нечетно.
Упражнение
Лихва до цент
Реши задачата в редактора си.
Колко става депозитът. Влагаш 1000 € при 1% лихва на месец, която се добавя към сумата. Прочети броя месеци и отпечатай сумата, закръглена до цент. Сметни я точно, с цели числа: 1000 € са 100 000 цента, а след n месеца сумата е 100 000 · 101ⁿ / 100ⁿ цента. Напиши свое бързо степенуване с int, без модул, и закръгли нагоре от половин цент.
120След 120 месеца: 3300,39 €При 120 месеца умноженията са малко и при наивния цикъл. Тук числата растат (101¹²⁰ има 241 цифри), затова всяко следващо умножение е по-скъпо, а пестенето на умножения важи още повече.
Още задачи
- Дневна лихва за 30 години. Същото с 0,01% на ден за 10 950 дни. Преброй умноженията на двете версии.
- Последните цифри. Отпечатай последните 6 цифри на 2¹⁰⁰⁰⁰⁰⁰⁰⁰⁰ (модул 1 000 000), без да строиш цялото число.
По-точно
Още за степента
По-точно, за любопитните
Числата на Фибоначи за O(log n). Същото повдигане на квадрат работи с матрица 2 × 2 вместо число: n-тото число на Фибоначи излиза от n-тата степен на [[1, 1], [1, 0]], с около 2 log₂ n умножения на матрици.
Елиптичните криви. Днес TLS договаря ключа най-често с елиптична крива (X25519). Там вместо умножение има „събиране на точки“, а вместо повдигане на квадрат „удвояване“, но алгоритъмът е същият: double-and-add по двоичните цифри на тайното число.
Постоянно време. Цикълът от урока умножава резултата само при цифра 1. Някой, който мери времето, може да познае колко единици има тайното число. Криптографските библиотеки затова правят еднаква работа за всяка цифра (Montgomery ladder). Затова криптография не се пише на ръка, а се ползва библиотека.
- Python:
pow(a, -1, m)(от 3.8) дава обратното на a по модул m.a ** nс цели числа е точно, но расте без граница.
Тест