C, PHP, VB, .NET

Степенуването е основна математическа операция и е очаквано да се използва често. Интуитивната реализация на алгоритъм за изчисление на [mathi]x^{n}[/mathi] би била „умножаваме [mathi]x[/mathi] по себе си [mathi]n[/mathi] на брой пъти“. Тоест нещо подобно на следното:

static BigInteger naivePow(BigInteger x, long n){
   BigInteger result = BigInteger.ONE;
   for (long i = 0; i < n; i++) {
      result = result.multiply(x);
   }
   return result;
}

С това решение при достигане до прекалено големи числа ще започне да се отнема сериозно време. Например на компютъра, на който пиша в момента тази статия (Intel Core-i5 9400) изчислението на [mathi]569890561^{121391}[/mathi] отнема около 8 секунди. Не може ли да се случи по-добре?

Основният проблем с бързината идва от там, че извършваме [mathi]n[/mathi] на брой умножения. Вместо това бихме могли да се възползваме от следното наблюдение:

• [mathi]x^{n}=x^{n/2}.x^{n/2}[/mathi] ако [mathi]n[/mathi] е четно;
• [mathi]x^{n}=x^{n/2}.x^{n/2}.x[/mathi] ако [mathi]n[/mathi] е нечетно.

Ето един бърз пример как това може да получи приложение. Представете си, че трябва да изчислим [mathi]5^8[/mathi]. С представената интуитивна реализация ще пресметнем [mathi]5.5.5.5.5.5.5.5=5^{8}[/mathi], т.е. 8 умножения. Вместо това можем да подходим със следната последователност:

1. Изчисляваме [mathi]5.5=5^{2}[/mathi];
2. Изчисляваме [mathi]5^{2}.5^{2}=5^{4}[/mathi] (като [mathi]5^{2}[/mathi] е вече пресметнато в предишната стъпка число);
3. Изчисляваме [mathi]5^{4}.5^{4}=5^{8}[/mathi]

Оказва се, че го направихме с едва 3 умножения. Остава да алгоритмизираме този процес. Най-популярната реализация се нарича „степенуване чрез повдигане на квадрат и умножение“ (square and multiply). При него се възползваме от спецификите на двоичната бройна система и най-вече, че операцията „повдигане на квадрат“ е изключително лека – просто се добавя 0 в края на числото. Ще го демонстрирам директно с пример. Нека да речем, че трябва да изчислим [mathi]5^{45}[/mathi]. Представяме числото 45 в бинарен вид:

[math]101101_{2}[/math]

Нека номерираме битовете с индекси отляво надясно с 1, 2, 3, 4, 5 и 6. Тогава за всеки от тях ще имаме следното:

1. [mathi]5^{1_{2}}=5^{1}[/mathi];
2. [mathi]5^{1_{2}}.5^{1_{2}}=5^{10_{2}}=5^{2}[/mathi];
3. [mathi]5^{10_{2}}.5^{10_{2}}.5^{1_{2}}=5^{101_{2}}=5^{5}[/mathi];
4. [mathi]5^{101_{2}}.5^{101_{2}}.5^{1_{2}}=5^{1011_{2}}=5^{11}[/mathi];
5. [mathi]5^{1011_{2}}.5^{1011_{2}}=5^{10110_{2}}=5^{22}[/mathi];
6. [mathi]5^{10110_{2}}.5^{10110_{2}}.5^{1_{2}}=5^{101101_{2}}=5^{45}[/mathi];

Вижда се как на всяка стъпка след първата се извършва повдигане на квадрат (когато конкретния бит е 0) или повдигане на квадрат и умножение по числото на 1-ва степен (когато конкретния бит е 1). Ето и как се реализира алгоритъма рекурсивно:

static BigInteger fastPowRecursive(BigInteger x, long n){
   if(n==0){
      return BigInteger.ONE;
   }
   BigInteger half = fastPowRecursive(x, n/2);
   BigInteger squared = half.multiply(half);
   if((n % 2)==0){
      return squared;
   } else {
      return squared.multiply(x);
   }
}

Ще забележите веднага, че тази реализация ще даде резултат изключително по-бързо, отколкото предишното решение. Итеративната реализация освен това ще използва и по-оптимално паметта:

static BigInteger fastPowIterative(BigInteger x, long n){
   BigInteger result = BigInteger.ONE;
   while(n > 0){
      if((n % 2) == 1){
         result = result.multiply(x);
      }
      x = x.multiply(x);
      n = n / 2;
   }
   return result;
}

Бихте могли допълнително да подобрите бързодействието ако използвате директно побитови операции.

П.П. Този алгоритъм освен това получава изключително често приложение в криптографията, където често се налагат изчисления от типа [mathi]x^{n} mod(m)[/mathi] (например при алгоритъма RSA). Стандартно може да се помисли първо да се изчисли степенуването [mathi]x^{n}=A[/mathi], а след да се намери остатъка от деление [mathi]A mod(m)[/mathi]. По-оптимален вариант би бил да се извършва делението по модул на всяка стъпка от алгоритъма – това няма да промени крайния резултат (проверете). Вярно е, че така броят на извършените операции се увеличава, но за сметка на това вече не се налага да се работи с много големи числа, т.е. заетата памет значително намалява.

П.П.2. Единственият недостатък на така показания алгоритъм е, че при единичен бит се извършва една операция повече (square + multiply), отколкото при нулев (само square). Това прави алгоритъмът податлив на т.нар. „side channel attacks“. Едно възможно решение е алгоритъмът на Монтгомери. Негова примерна реализация е следната:

static BigInteger montgomeryPow(BigInteger x, long n) {
   BigInteger x1 = x;
   BigInteger x2 = x.multiply(x);
   // взима броя битове на числото (премахвайки нулевите в началото)
   int n_bits_len = Long.SIZE - Long.numberOfLeadingZeros(n);
   // обхождаме числото бит по бит, пропускайки най-старшия
   for(int i = n_bits_len-2; i >= 0; i--){
      // (n & (1 << i)) >> i връща i-ти бит на числото n
      if (((n & (1 << i)) >> i) == 0) {
         x2 = x1.multiply(x2);
         x1 = x1.multiply(x1);
      } else {
         x1 = x1.multiply(x2);
         x2 = x2.multiply(x2);
      }
   }
   return x1;
}