Writing your own exponential power function with d

Question

So, I want to write a function in code using some sort of algorithm to calculate any number to any power, including decimals. I use JavaScript and it already has an inbuilt pow function:

Math.pow(2, 0.413) // 2^0.413 = 1.331451613236371, took under 1 second.

Now I want to write my own like this:

function pow(x, y) {
    // Algorithm
}

This is a function that calculates the square root of any number (x^0.5), and it's very accurate with only 10 loops:

function sqrt(x, p) { // p = precision (accuracy)
    var a = 1;
    var b = x;

    while (p--) {
        a = (a + b) / 2
        b = x / a
    }

    return a
}

Is there any simple formula to calculate any exponential?

If there isn't a simple one, is there a hard one?

If the solution is slow, how can JavaScript's pow estimate under a single second?

Answer 1

Those are some really nice examples, here is a simpler one too.

function exponential(a,b){
    var c = 1;
    for(var i=1; i<=b; i++){
        c = c * a;
    }
    return c;
}

now call the function:

exponential(2,4);

Edit: It only works on integer, but it's simple and quick.

Answer 2

I checked this post, but it worked only for whole numbers (1,2,3... not 0.1, 0.3...)

Recursive power function: Why does this work if there's no initial return value?

Then,

I got this from here: Algorithm for pow(float, float)

function power(x,n) {
    if(n === 0) return 1;
    if(n === -1) return 1/x;
    if(n === 1) return x;
    return Math.exp(n*Math.log(x))
}

console.log(power(2,3.5));

I added some basic checks (n===0)... To fasten things up in case.

Flexo sums it up:

The general algorithm tends to be computing the float power as the combination of the integer power and the remaining root. The integer power is fairly straightforward, the root can be computed using either Newton - Raphson method or Taylor series. IIRC numerical recipes in C has some text on this. There are other (potentially better) methods for doing this too, but this would make a reasonable starting point for what is a surprisingly complex problem to implement. Note also that some implementations use lookup tables and a number of tricks to reduce the computation required.

http://mathworld.wolfram.com/NewtonsMethod.html

http://mathworld.wolfram.com/TaylorSeries.html

https://en.wikipedia.org/wiki/Logarithm#Power_series

https://rads.stackoverflow.com/amzn/click/0521431085

Answer 3

Heres a nice algorithm for positive integer powers, it starts by dealing with some simple cases and then uses a loop testing the binary bits of the exponent. For example to find 3^11 11 in binary is 1011 so the stages in the loop are

• bitMask = 1011, evenPower = 3, result = 3
• bitMask = 101, evenPower = 3*3 = 9, result = 3*9 = 27
• bitMask = 10, evenPower = 9*9 = 81, result = 27
• bitMask = 1, evenPower = 81*81 = 6561, result = 27*6561 = 177147

That is the evenPower squares at each loop, and the result gets multiplied by the evenPower if the bottom bit is 1. The code has been lifted from Patricia Shanahan’s method http://mindprod.com/jgloss/power.html which in turn has its roots in Kunth and can be traced back to 200 BC in india.

/**
 * A fast routine for computing integer powers.
 * Code adapted from {@link <a href="http://mindprod.com/jgloss/power.html">efficient power</a>} by Patricia Shanahan [e-mail]
 * Almost identical to the method Knuth gives on page 462 of The Art of Computer Programming Volume 2 Seminumerical Algorithms.
 * @param l number to be taken to a power.
 * @param n power to take x to. 0 <= n <= Integer.MAX_VALUE
 * Negative numbers will be treated as unsigned positives.
 * @return x to the power n
 * 
 */
public static final double power(double l,int n)
{
    assert n>=0;

    double x=l;
    switch(n){
    case 0: x = 1.0; break;
    case 1: break;
    case 2: x *= x; break;
    case 3: x *= x*x; break;
    case 4: x *= x; x *= x; break;
    case 5: { double y = x*x; x *= y*y; } break;
    case 6: { double y = x*x; x = y*y*y; } break;
    case 7: { double y = x*x; x *= y*y*y; } break;
    case 8: x *= x; x *= x; x *= x; break;
    default:
    {
        int bitMask = n;
        double evenPower = x;
        double result;
        if ( (bitMask & 1) != 0 )
            result = x;
        else
            result = 1;
        bitMask >>>= 1;
        while ( bitMask != 0 ) {
            evenPower *= evenPower;
            if ( (bitMask & 1) != 0 )
                result *= evenPower;
            bitMask >>>= 1;
        } // end while
        x = result;
    }
    }
    return x;
}

For a real exponent you will basically need ways of finding exp and log. You can use Taylor series which are the simplest to get but there are much better method. We have

exp(x) = 1 + x + x^2/2 + x^3/6 + x^4/24 + x^5/120 + x^6/6! + ...

ln(1+x) = x - x^2 /2 + x^3 /3 - x^4 / 4 + x^5 / 5 - x^6/6 + ... |x|<1

To find x^y note ln(x^y) = y*ln(x). Now we need to get the argument in the right range so we can use our power series. Let x = m * 2^ex, the mantissa and exponent chosen so 1/sqrt(2)<= m < sqrt(2) and ln(m*2^ex) = ln(m) + ex*ln(2). Let h = m-1 and find ln(1+h).

Using java and floats as there is an easy way to find the internals of the IEEE representation (I've used float as there are fewer bits to cope with)

int b = Float.floatToIntBits(x);
int sign = (b & 0x80000000) == 0 ? 1 : -1;
int mattissa = b & 0x007fffff;
int ex = ((b & 0x7f800000) >> 23 ) - 127;

in javascript it might be easiest to us Number.toExponential and parse the results.

Next construct a number z in the desired range 1/sqrt(2) < z < sqrt(2)

int bits = mattissa | 0x3f800000;
float z = Float.intBitsToFloat(bits);
if(z>root2) { 
    z = z/2;
    ++ex;
}

Use this function to find the log of 1+x using a taylor series

static float ln1px(float x) {
    float x_2 = x*x; // powers of x
    float x_3 = x_2 * x;
    float x_4 = x_3 * x;
    float x_5 = x_4 * x;
    float x_6 = x_5 * x; 
    float res = x - x_2 /2 + x_3 /3 - x_4 / 4 + x_5 / 5 - x_6/6;
    return res;
}

this seems to be good to three significant figures, often much better when x is close to 0.

The log of our number x can then be found

float w = z - 1;
float ln_z = ln1px(w);
float ln_x = ln_z + ln2 * ex;
System.out.println("ln "+ln_x+"\t"+Math.log(x));

Now to the actual power if we write y = n + a where n is an integer and a is fractional. So x^y=x^(n+a) = x^n * x^a. use the first algorithm in this answer to find the x^n. Writing x=m*2^ex then ln((m*2^ex)^a) = yln(m) + yex*ln(2) and

x^a=exp(ln((m*2^ex)^a)) = exp(a * ln(m)) * exp(a * ln(2))^ex

the two exponential terms have fairly small values so the taylor series should be good.

We need one function for the taylor series of the exponential function

static float exp(float x) {
    float x_2 = x*x; // powers of x
    float x_3 = x_2 * x;
    float x_4 = x_3 * x;
    float x_5 = x_4 * x;
    float x_6 = x_5 * x; 
    float res = 1+ x + x_2 /2 + x_3 /6 + x_4 / 24 + x_5 / 120 + x_6/ 720;
    return res;
}

finally we can put the pieces together

// Get integer and fractional parts of y
int n = (int) Math.floor(y);
float a = y-n;

float x_n = power(x,n);         // x^n
float a_ln_m = a * ln_z;        // ln(m^a) = a ln(m)
float a_ln_2 = a * ln2;         // ln(2^a) = a ln(2)
float m_a = exp(a_ln_m);        // m^a = exp(a ln(m))
float _2_a = exp(a_ln_2);       // 2^a = exp(a ln(2))
float _2_a_ex = power(_2_a,ex); // (2^ex)^a = 2^(a*ex) = (2^a)^ex 
float x_a = m_a * _2_a_ex;      // x^a = m^a * 2^(a*ex)

float x_y = x_n * x_a;          // x^y = x^n * x^a

System.out.println("x^y "+x_y+"\t"+Math.pow(x,y));

That should be the complete program, you need some smarts to cope with negative arguments etc.

Note this is not particularly accurate as I've only used a few terms of the taylor series. Other SO questions have more detailed answers How can I write a power function myself?