Бързо степенуване
Накратко: Бързото степенуване смята aⁿ с O(log n) умножения, като повдига на квадрат и гледа двоичните цифри на n; с остатък след всяко умножение смята и степен по модул.
Проблем
Ключът на HTTPS връзката
Когато браузърът отваря HTTPS страница, той и сървърът се договарят за общ таен ключ, без да го пращат по мрежата. В класическия вариант (Diffie–Hellman) всеки смята степен по модул: gᵃ mod p, където a е тайно число, дълго десетки цифри. Същата сметка стои в RSA подписите и в hash на текст, който умножава по степени на основата.
Наивният начин е цикъл, който умножава n пъти:
static long powerLoop(long a, long n, long m) { long result = 1; for (long i = 0; i < n; i++) { 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: това е дъното.
static long power(long a, long n, long m) { if (n == 0) { return 1; } long half = power(a, n / 2, m); long square = half * half % m; if (n % 2 == 1) { return square * a % m; } return square;}
System.out.println(power(3, 13, 1_000_000_007));System.out.println(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 влиза в резултата.
static long power(long a, long n, long m) { long result = 1; a %= m; while (n > 0) { if (n % 2 == 1) { result = result * a % m; } a = a * a % m; n /= 2; } return result;}
System.out.println(power(3, 13, 1_000_000_007));System.out.println(power(5, 1_000_000_000, 1_000_000_007));1594323142848001Умноженията са колкото цифрите на n (по едно повдигане на квадрат на цифра) плюс единиците в него (по едно умножение на резултата). Това е версията, която ползват библиотеките: няма рекурсия и не зависи от размера на стека.
Виж го
Квадратите влизат в резултата
static long power(long a, long n, long m) { long result = 1; a %= m; while (n > 0) { if (n % 2 == 1) { result = result * a % m; } a = a * a % m; n /= 2; } return result;}
System.out.println(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, а бързата и на милиард:
static long multiplications = 0;
static long powerLoop(long a, long n, long m) { long result = 1; for (long i = 0; i < n; i++) { multiplications++; result = result * a % m; } return result;}
static long power(long a, long n, long m) { long result = 1; a %= m; while (n > 0) { if (n % 2 == 1) { multiplications++; result = result * a % m; } multiplications++; a = a * a % m; n /= 2; } return result;}
long m = 1_000_000_007;for (long n : new long[] {10, 1_000, 100_000, 10_000_000}) { multiplications = 0; powerLoop(5, n, m); long slow = multiplications; multiplications = 0; power(5, n, m); System.out.println("n = " + n + ": цикъл " + slow + ", бързо " + multiplications + " умножения");}multiplications = 0;power(5, 1_000_000_000, m);System.out.println("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 |
static long power(long a, long n, long m) { long 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
Кой показател иска повече умножения?
static long multiplications = 0;
static long power(long a, long n, long m) { long result = 1; a %= m; while (n > 0) { if (n % 2 == 1) { multiplications++; result = result * a % m; } multiplications++; a = a * a % m; n /= 2; } return result;}
for (long n : new long[] {15, 16}) { multiplications = 0; power(3, n, 1_000_000_007); System.out.println(n + ": " + multiplications + " умножения");}15 = 1111: четири цифри и четири единици, 8 умножения. 16 = 10000: пет цифри и една единица, 6 умножения. По-голям показател не значи повече работа: важат цифрите и единиците.
Идеята
Когато произведението не се побира
Остатъкът държи числата под m, но произведението на две такива числа е до m². Ако m² не се побира в типа, сметката греши, и то тихо. Ето същата степен с по-голям модул:
long m = 10_000_000_019L;System.out.println("long: " + power(5, 1_000_000_000, m));BigInteger exact = BigInteger.valueOf(5).modPow(BigInteger.valueOf(1_000_000_000), BigInteger.valueOf(m));System.out.println("BigInteger: " + exact);long: -7139217113BigInteger: 475867028long стига до 9,2 · 10¹⁸. При m = 10 000 000 019 произведението е до 10²⁰, преминава границата и се обръща в отрицателно число, без грешка. Безопасно е, докато m е под около 3 · 10⁹ (m² под 2⁶³); 1 000 000 007 е под границата. За по-голям модул има BigInteger.modPow, а Math.multiplyExact поне хвърля ArithmeticException вместо да мълчи.Идеята
Вграденото и показател от 77 цифри
В истинския код степента по модул не се пише на ръка:
BigInteger secret = BigInteger.TWO.pow(255).subtract(BigInteger.valueOf(19));BigInteger p = BigInteger.TWO.pow(127).subtract(BigInteger.ONE);System.out.println("Цифри в показателя: " + secret.toString().length());System.out.println("Умножения при повдигане на квадрат: " + (secret.bitLength() + secret.bitCount()));System.out.println("5^secret mod p = " + BigInteger.valueOf(5).modPow(secret, p));Цифри в показателя: 77Умножения при повдигане на квадрат: 5085^secret mod p = 97285997395055909576908201857106207893BigInteger.modPow прави повдигане на квадрат и взима остатъка по пътя. 508 умножения на числа от 39 цифри вместо 6 · 10⁷⁶: затова ключът се договаря за милисекунди. Библиотеките броят по-малко и от това, защото гледат по няколко цифри наведнъж, но идеята е същата.
Пробвай
Намери грешката
Тази версия връща верния резултат. Само че при n = 1 000 000 прави 2 048 575 умножения вместо 27. Кой ред е виновен?
Избери реда с грешката.
Половината се смята два пъти. Всяко ниво удвоява извикванията: 2, 4, 8… и на дъното са около n. Дълбочината пак е log₂ n, затова стекът не забелязва нищо, а работата се връща на O(n). Поправката е половината в променлива и после half * half.
Пробвай
Попълни
Кога текущото a влиза в резултата?
static long power(long a, long n, long m) {
long 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ⁿ цента. Напиши свое бързо степенуване с BigInteger, без модул, и закръгли нагоре от половин цент.
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). Затова криптография не се пише на ръка, а се ползва библиотека.
- Java:
BigInteger.pow(int)е степен без модул.modPowприема и отрицателен показател, ако основата има обратно по модула (иначеArithmeticException).
Тест