Algorytmy

Szybkie Potęgowanie

  • 2010-10-18 22:14
  • 2 komentarze
  • 10577 odsłon
  • Oceń ten tekst jako pierwszy
Spis treści

     1 Opis
     2 Uwagi
     3 Przykładowe implementacje
          3.1 Java
          3.2 C


Opis


Algorytm szybkiego potęgowania jest bardzo dobrym rozwiązaniem jeżeli musimy obliczać całkowite potęgi dużych liczb. Praktycznie używa się go w implementacjach  szyfrowania RSA gdzie wykorzystywany jest przy obliczaniu reszty z dzielenia potęgi przez liczbę.

Opiera się o właściwości potęgowania:

x<sup>{n} = x * x</sup>{n-1}

Wykonanie działania x^n metodą naiwną pociąga za sobą konieczność wykonania n mnożeń. W przypadku dużych wartości nzłożoność obliczeniowa zaczyna rosnąć wykładniczo w stosunku do liczby cyfr n. Korzystając z wymienionej wyżej własności można zauważyć, że do obliczenia x^n wystarczy wykonanie 1+ n/2 mnożeń. Całość można zatem zapisać w postaci rekurencyjnej:

potęga(x, n)
  jeżeli n = 0 
    zwróć 1;
  jeżeli n nieparzysta
    zwróć x * potęga (x, n-1);
  jeżeli n parzysta
    a = potęga (x, n/2); 
    zwróć a * a


Rekurencję tą można zastąpić wersją iteracyjną opartą o przesunięcia bitowe.
potęga(x, n)
   a = 1;
   dopóki n > 0
     jeżeli n & 1
       a = a * x
     x = x^2
     n = n >>1  
  zwróć a;


Wersja ta jest szybsza od rekurencyjnej, ale jej zrozumienie wymaga znajomości arytmetyki bitowej. Jej zaletą jest mniejsze zużycie pamięci z uwagi na wykonywanie obliczeń "w miejscu" to znaczy bez sięgania po dodatkowe ramki stosu (co jest bolączką wszystkich algorytmów rekurencyjnych).

Uwagi


Algorytm posiada pewne ograniczenia. Najpoważniejszym jest konieczność użycia liczby całkowitej jako potęgi. Można obejść to ograniczenie rozkładając liczbę zmiennoprzecinkową na część całkowitą i ułamkową, podniesieniu x najpierw do potęgi całkowitej za pomocą algorytmu szybkiego potęgowania, a następnie pomnożenie wyniku przez x podniesione do potęgi równej części ułamkowej za pomocą innego algorytmu. Należy jednak pamiętać, że tego typu działanie jest narażone na błędy wynikające z niedokładności zapisu liczby zmiennoprzecinkowej.
Algorytmu nie opłaca się stosować do obliczeń w przypadku liczb znanych w momencie kompilacji. W takim przypadku kompilatory, co do zasady, wyliczają wartość wyrażenia i wstawiają wynik do kodu wynikowego.

Przykładowe implementacje


Java


public class QuickExponentiating {
 
        public static double pow(double base, double power) {
                if (power == 0)
                        return 1;
                if (power % 2 == 1)
                        return pow(base, power - 1) * base;
                else {
                        double a = pow(base, power / 2);
                        return a * a;
                }
        }
 
        public static double bit_pow(double base, Long power) {
                double res = 1.0;
                while (power > 0) {
                        int i = power.byteValue() & (byte) 1;
                        if (i != 0)
                                res *= base;
                        base *= base;
                        power >>= 1;
                }
                return res;
        }
 
        public static double bit_pow(BigDecimal base, BigInteger power) {
                BigDecimal res = BigDecimal.ONE;
                while (power.compareTo(BigInteger.ZERO) > 0) {
                        BigInteger i = power.and(BigInteger.ONE);
                        if (!i.equals(BigInteger.ZERO))
                                res = res.multiply(base);
                        base = base.multiply(base);
                        power = power.shiftRight(1);
                }
                return res.doubleValue();
        }
}


Uwaga, kod rekurencyjny oparty o klasy BigDecimal i BigInteger może powodować StackOverflowError. Wersja iteracyjna jest pozbawiona tego mankamentu.

C


Wersja ricewinda z komentarzy.
double bin_pow(double base, unsigned int exp)
{
        double res = 1.0;
        while (exp > 0)
        {
                if (exp & 1)
                        res *= base;
                base *= base;
                exp >>= 1;
        }
        return res;
}
 
double rek_pow(double base, unsigned int exp)
{
        if (exp == 0)
                return 1;
        if (exp % 2 == 1)
                return rek_pow(base, exp - 1) * base;
        else
        {
                double a = rek_pow(base, exp / 2);
                return a * a;
        }
}

2 komentarze

Koziołek 2010-10-18 16:32

@rincewind, przysiądę i wieczorem dopiszę co trzeba. Wersję w C twojego autorstwa też wrzucę a co :)

rincewind 2010-10-18 16:25

Rekurencja powoduje dość spory narzut czasowy, a ten algorytm można ładnie przełożyć na wersję iteracyjną (ze względu na rekurencję ogonową). Poniżej kod w C++:

double bin_pow(double base, unsigned int exp)
{
    double res = 1.0;
    while (exp > 0)
    {
        if (exp & 1)
            res *= base;
        base *= base;
        exp >>= 1;
    }
    return res;
}


Jest około 3x szybszy na podstawie tego "benchmarka" (kompilacja pod GCC 4.5 z flagą -O2):

#include <iostream>
#include <ctime>
 
double bin_pow(double base, unsigned int exp)
{
    double res = 1.0;
    while (exp > 0)
    {
        if (exp & 1)
            res *= base;
        base *= base;
        exp >>= 1;
    }
    return res;
}
 
double rek_pow(double base, unsigned int exp)
{
    if (exp == 0)
        return 1;
    if (exp % 2 == 1)
        return rek_pow(base, exp - 1) * base;
    else
    {
        double a = rek_pow(base, exp / 2);
        return a * a;
    }
}
 
int main()
{
    double x;
    int exp;
    double res;
 
    unsigned int start = clock();
    for (int i = 0; i < 1000; ++i)
        for (x = 2.0; x < 1000.0; x += 1.0)
            for (exp = 2; exp < 100; ++exp)
                res = bin_pow(x, exp);
    std::cout << "[bin_pow] Ticks: " << (clock() - start) << std::endl;
 
    start = clock();
    for (int i = 0; i < 1000; ++i)
        for (x = 2.0; x < 1000.0; x += 1.0)
            for (exp = 2; exp < 100; ++exp)
                res = rek_pow(x, exp);
    std::cout << "[rek_pow] Ticks: " << (clock() - start) << std::endl;
}


Wynik na mojej maszynie:

> pow
[bin_pow] Ticks: 4049
[rek_pow] Ticks: 12491


Po zmianie na wersję iteracyjną odpada też niebezpieczeństwo Stack Overflow, bo wszystko wykonuje się w miejscu.