A..2 Potenzen mit ganzzahligen Exponenten

Mit der Funktion pow() aus der Standardbibliothek können Potenzen aus zwei Gleitkommazahlen gebildet werden.

Dies ist aber in vielen Fällen unnötiger Rechenaufwand, da Potenzen mit ganzzahligen Exponenten schneller berechnet werden können.

In [DEK II] und [Amm] sind recht schnelle Verfahren dazu angegeben.

Beispielhaft wird hierfür folgende Funktion angegeben:

/* popowi.c 30. 7.95 kw
 * Berechnung von Potenzen mit ganzzahligen Exponenten nach:
 *  Donald E. Knuth: The Art of Computer Programming
 *  Vol. II (Seminumerical Algorithms) Second Edition
 *  Abschnitt 4.6.3, S. 441, Algorithmus A.
 * Je nach Option kann man eine Funktion fuer eine ganzzahlige Basis
 * (#define TYP LONG) oder eine gebrochene Basis (#define TYP DOUBLE)
 * erzeugen.
 * Der Funktionsname ist dementsprechend lpowi() oder dpowi().
 *
 * Da in dem genannten Algorithmus weder n noch x nach Start der
 * Rechnung benoetigt werden, kann man auf die lokalen Variablen N und
 * Z (aus Knuth) verzichten. Daher folgende Variablenentsprechung:
 *      Knuth:      | powi.c:
 *
 *      n, N        | n
 *      Y           | y
 *      x, Z        | x
 *
 *
 * Zur leichteren Orientierung sind als Kommentar die entsprechenden
 * Fragmente des Originaltextes eingefuegt.
 */

/* LONG oder DOUBLE:
 */
#ifndef TYP
#define TYP         DOUBLE
#endif

#define LONG        1
#define DOUBLE      2

#if     TYP==DOUBLE
#define BASIS_T     double      /* Typ                              */
#define POW         dpowi       /* Funktionsname                    */
#elif   TYP==LONG
#define BASIS_T     long        /* Typ                              */
#define POW         lpowi       /* Funktionsname                    */
#else
error "TYP nicht LONG oder DOUBLE!"
#endif

/* TEST oder NTEST fuer oder gegen Testprogramm:
 */
#define TEST

BASIS_T POW( BASIS_T x, int n )
{
  BASIS_T     y = 1;          /* A1. [Initialize.]                */
  do
  {
    if( n & 0x0001 )            /* A2. [Halve N.] "wether N was
                                 *                even or odd"
                                 */
    {
      n >>= 1;            /* A2. [Halve N.] "N <- N/2"
                           */
      /* A3. [Multiply Y by Z.] "Set Y <- Z times Y"
       */
      y *= x;

      /* A4. [N = 0?] "If n=0, the algorithm
       *              terminates with Y as the answer."
       */
      if( !n ) return y;
    }
    else
    {
      n >>= 1;            /* A2. [Halve N.] "N <- N/2"
                           *                "if N was even,
                           *                skip to section
                           *                A5."
                           */
    }
    /* A5. [Square Z.]         "Set Z <- Z times Z"             */
    x *= x;
  }
  while( 1 );               /* "and return to step A2."         */
}

#ifdef TEST

#if     TYP==LONG
#define WRFMT           "%10ld"
#elif   TYP==DOUBLE
#define WRFMT           "%10.3f"
#else
error "TYP nicht LONG oder DOUBLE!"
#endif

#include <time.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define IMAX                10000000

void testpow( BASIS_T b, int exp )
{   BASIS_T     ret;
 double      dret;

 long        i;
 clock_t     start,
   ende;

 printf( "Testprogramm fuer Potenzieren mit "
         "ganzzahligen Exponenten.\n"
         );

 /* Test mit Originalroutine:                                    */
 start = clock();
 for( i=0; i<IMAX; i++ )
   dret = pow( (double)b, (double)exp );
 ende = clock();
 printf( "Original: pow(%10.3f,%7.1f) = %10.3f; in %.9f sec\n",
         (double)b,
         (double)exp,
         dret,
         (double)( ende - start )/((double)CLOCKS_PER_SEC*IMAX)
         );

 /* Test mit Algorithmus aus Knuth:                              */
 start = clock();
 for( i=0; i<IMAX; i++ )
   ret = POW( b, exp );
 ende = clock();
 printf( "Knuth:    pow(" WRFMT ",%7d) = " WRFMT "; in %.9f sec\n",
         b,
         exp,
         ret,
         (double)( ende - start )/((double)CLOCKS_PER_SEC*IMAX)
         );
 printf( "\n" );
}

int main()
{       testpow( 2, 3 );
 testpow( 2, 10 );
 return EXIT_SUCCESS;
}
#endif

In diesem Quelltext hat man mit #define TYP DOUBLE oder #define TYP LONG die Auswahl, eine von zwei verschiedenen Funktionen kompilieren zu lassen: Entweder die Funktion long lpowi( long x, int n ) oder double dpowi( double x, int n ).

Weiterhin hat man mit #define TEST die Möglichkeit, ein Testprogramm erzeugen zu lassen. #define NTEST unterbindet dies; dann wird nur eine Funktion übersetzt.

Auf einem bestimmten Rechner (Laptop Pentium 3, 1.1GHz, Linux, gcc 2.95.3) werden von der long-Version folgende Ausgaben für $2^3$ und $2^{10}$ erzeugt:

Testprogramm fuer Potenzieren mit ganzzahligen Exponenten.
Original: pow(     2.000,    3.0) =      8.000; in 0.000000172 sec
Knuth:    pow(         2,      3) =          8; in 0.000000011 sec

Testprogramm fuer Potenzieren mit ganzzahligen Exponenten.
Original: pow(     2.000,   10.0) =   1024.000; in 0.000000176 sec
Knuth:    pow(         2,     10) =       1024; in 0.000000018 sec

Dieser Algorithmus ist nicht so universell wie die Funktion pow() aus der Standardbibliothek, aber dafür in ihrem Einsatzgebiet wesentlich schneller.

Die Moral dieses Abschnitts ist damit, daß man die Rechenleistung von Programmen teilweise erheblich steigern kann, wenn man brauchbare Algorithmen für einen bestimmten Einsatzzweck entwickelt.

In der double-Version ist die Beschleunigung ebenfalls kräftig:

Testprogramm fuer Potenzieren mit ganzzahligen Exponenten.
Original: pow(     2.000,    3.0) =      8.000; in 0.000000164 sec
Knuth:    pow(     2.000,      3) =      8.000; in 0.000000014 sec

Testprogramm fuer Potenzieren mit ganzzahligen Exponenten.
Original: pow(     2.000,   10.0) =   1024.000; in 0.000000189 sec
Knuth:    pow(     2.000,     10) =   1024.000; in 0.000000023 sec

Auf anderen Rechnern kann der Vergleich aber ganz anders aussehen; solche Optimierungen sind in aller Regel nicht universell übertragbar.

In [Amm] sind teilweise noch etwas schnellere Algorithmen angegeben.