Floating point power function without recursion

Question

GLSL doesn't have a pow function for double precision types. Classical implementations of this function make use of recursion. Since recursion is not allowed in GLSL I need a recursion free implementation.

I also need the exp function.

Answer 1

This problem can simply be solved by seperating it into two sub problems.

1. Seperate the whole and fractional part of the exponent
2. Calculate the power with the whole part (which is rather trivial)
3. Calculate the power with the fractional part (using an approximation)
4. Multiply the results
5. Invert the result if the initial exponant was negative

/// This is a simple non recursive implementation for pow, where the exponent is an integer.
double pow(double base, int power) {
    double res = 1.0;// Initialize result

    while (power > 0.0) {
        // If power is odd, multiply x with result
        if (power % 2 == 1) {
            res *= base;
        }

        // n must be even now
        power /= 2;
        base *= base;// Change x to x^2
    }

    return res;
}

// The desired precision. More precision results in slower execution.
const double EPSILON = 0.00000001;

double pow(double base, double power) {
    // We want to ignore negative exponents for now. We invert our result at the if necessary.
    bool negative = power < 0.0LF;
    if (negative) {
        power *= -1.0LF;
    }

    // Seperate the whole and fractional parts.
    double fraction = power - int(power);
    int integer = int(power - fraction);

    // The whole part is easily calculated.
    double intPow = pow(base, integer);

    // The fractional part uses an approximation method.
    double low = 0.0LF;
    double high = 1.0LF;

    double sqr = sqrt(base);
    double acc = sqr;
    double mid = high / 2.0LF;

    while (abs(mid - fraction) > EPSILON) {
        sqr = sqrt(sqr);

        if (mid <= fraction) {
            low = mid;
            acc *= sqr;
        } else {
            high = mid;
            acc *= (1.0LF / sqr);
        }

        mid = (low + high) / 2.0LF;
    }

    // Exponential rules allow us to simply multiply our results.
    double result = intPow * acc;

    // If we started with a negative exponent we invert the result.
    if (negative) {
        return 1.0LF / result;
    }

    return result;
}

const double E = 2.7182818284590452353602874LF;

/// e^x is as simple as that.
double exp(double power) { return pow(E, power); }

Performance

Average iterations for exp(random(0.0, 100.0)):

EPSILON = 0.0001     -> 16.0763
EPSILON = 0.00001    -> 19.4108
EPSILON = 0.000001   -> 22.6639
EPSILON = 0.0000001  -> 26.0477
EPSILON = 0.00000001 -> 29.3929