2

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.

Jhonny007
  • 1,698
  • 1
  • 13
  • 33
  • While your answer is quite good, you can always simulate recursion using a stack. It's usually pretty simple too. – President James K. Polk Dec 10 '19 at 15:05
  • see [Power by squaring for negative exponents](https://stackoverflow.com/a/30962495/2521214) but I would first try `pow(x,y) = exp2(y*log2(x))` or any other `exp,log` base IIRC they are implemented in GLSL – Spektre Dec 11 '19 at 08:17
  • While they are implemented for GLSL they are only available for single precision floats not for double precision floats. Otherwise I could easily use the default `pow` function provided by GLSL. – Jhonny007 Dec 24 '19 at 15:48

1 Answers1

4

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
Jhonny007
  • 1,698
  • 1
  • 13
  • 33