Re: fast power function

In article <Xns9A2CA71CBAD72toelavabitcom@194.125.133.14>,
toe@lavabit.com says...


It already has such overloads, though there's no guarantee that they
work any differently than for floating point exponents.


This isn't a particularly good algorithm unless exp is _quite_ small
(low double digits at most). For an N-bit exponent, this takes a maximum
of 2^N-1 multiplications. Assuming values of exp were randomly
distributed, it would take about 2^(N-1) multiplications on average.

You can do a lot better than that, taking roughly 2N multiplications in
the worst case, and something like 1.5N multiplications in the typical
case (actually, it should be somewhat less than that, but I haven't
bothered to figure up the exact number -- in any case, it's linear on
the number of bits instead of exponential.

Here's an implementation with some test code and instrumentation:

#include <vector>
#include <numeric>

unsigned mults;

template <class T>
T mult(T a, T b) {
    mults++;
    return a * b;
}

#define bits(object) (sizeof(object)*CHAR_BIT)

template <class T>
T my_pow(T X, unsigned Y) {
    std::vector<T> sq;

    sq.reserve(bits(Y));
    
    T current_square(X);
    mults = 0;

    while (Y) {
        if (Y&1)
            sq.push_back(current_square);
        Y>>=1;
        current_square *= current_square;
        mults++;
    }
    T result = std::accumulate(sq.begin(), sq.end(), T(1), mult<T>);

    return result;
}

#ifdef TEST

#include <iostream>
#include <cmath>

int main() {
    unsigned max_mults = 0;
    unsigned total_mults = 0;
    unsigned numbers = 0;

    for (double i=0; i<10; i++)
        for (int j=0; j<300; j++) {
            std::cout << i << "^" << j << "=" << my_pow(i, j);
            std::cout << "\t" << mults << " multiplications\n";
            if (std::pow(i,j) != my_pow(i,j))
                std::cout << "^Error!\n";
            if (mults > max_mults)
                max_mults = mults;
            numbers ++;
            total_mults += mults;
        }
    std::cout << "Maximum Multiplications: " << max_mults;
    std::cout << "\nAverage Multipliations: " << (double)total_mults /
numbers;
    return 0;
}

#endif

When I run this, I get a worst-case of 16 and an average of 11.23
multiplications per call. In every case, the result matches the standard
library exactly -- which gives a strong indication that the standard
library is implementing the code the same way. I'd expect to see at
least a few minor differences if the standard library were doing the job
by taking the logarithm and multiplying.

--
    Later,
    Jerry.

The universe is a figment of its own imagination.