Szybkie Potęgowanie

Koziołek

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 xn 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</sup>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 komentarzy

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

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.