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?

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?

• javascript
• algorithm
• math

• 1 "If the solution is slow, how can JavaScript's pow estimate under a single second?" - it's written in C, that's how. like this. – The Paramagnetic Croissant Commented Jun 1, 2014 at 17:32
• Course it's written in C, but what is C written in? – user2039981 Commented Jun 1, 2014 at 17:34
• 1 I'm not sure if you're joking or not, but c isn't written in anything. It is a language. You might be asking what the c piler is written in, and the answer is, of course, it could be anything, but in practice, paradoxically, it is usually mostly c. Native functions in Javascript aren't fast because of the language they're written in, they're fast because they're native code, often highly optimized. They could be written in any piled language, or assembly, or whatever... – gwcoffey Commented Jun 1, 2014 at 17:37
• Yes whatever, but in the bottom of everything, there is an algorithm, and that's what I want. – user2039981 Commented Jun 1, 2014 at 17:39
• @Murplyx you can read the source code of the exp() and log() functions, they are right next to the source file I've linked to. – The Paramagnetic Croissant Commented Jun 1, 2014 at 17:42

3 Answers 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./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 puting integer powers.
 * Code adapted from {@link <a href="http://mindprod./jgloss/power.html">efficient power</a>} by Patricia Shanahan [email protected]
 * 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 plete 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?

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.

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 puting the float power as the bination of the integer power and the remaining root. The integer power is fairly straightforward, the root can be puted 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 plex problem to implement. Note also that some implementations use lookup tables and a number of tricks to reduce the putation required.

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

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

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

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