how does the pow function work?

Question

I am looking to write my own power function to work with NSDecimalNumbers and exponents that are not whole numbers. I first tried to use a combination of newtons method and the built in integer power method, but due to newtons method i am getting overflow errors when I have exponents with more than 2 decimals. So I thought maybe the float value pow function might serve as a good model for my own function. So I was wondering if anyone knows where I can fond some sort of documentation on the inner workings of the pow function?

Edit:

@wombat57, those links look like they could be what I am looking for however I have no idea to read them. The algorithm you suggest is in fact what I am using. the overflow comes from newtons method due to very large exponents. Because I am getting exponents in decimal form I have to convert it to a fraction first. the only way of ding this in code, as far as I know, multiplying the decimal by ten until you have a whole number, and using that as the numerator. Doing this you get exponents of 100+ for numbers with 3 or more decimals. this causes an overflow error.

Answer 1

EDIT 1: Here are links to the actual source

http://opensource.apple.com/source/Libm/Libm-2026/Source/Intel/expf_logf_powf.c http://opensource.apple.com/source/Libm/Libm-315/Source/ARM/powf.c

I got the links from this question, which has a bunch of relevant discussion

self made pow() c++

This page describes an algorithm: Link. x^(1/n) = the nth root of x, and x^mn = (x^m)^n. Thus, x^(m/n) = (the nth root of x)^m. Arbitrary roots can be calculated with Newton's method. Integer powers can be calculated with Exponentiation by squaring. For irrational exponents, you can use increasingly accurate rational approximations until you get the desired number of significant digits.

EDIT 2:

Newton's method involves raising your current guess to the power of the root that you're trying to find. If that power is large, and the guess is even a little too high, this can result in overflow. One solution here is to identify this case. If overflow ever occurs, this means that the guess was too high. You can solve the problem by (whenever a guess results in overflow), setting the current guess to a value between the last guess that did not overflow and the current guess (you may have to do this several times). That is, whenever Newton's method overflows, do a binary search down toward the last guess that did not overflow. Here's some python that implements all of this:

def nroot(n, b, sig_figs = 10):
    g1 = 1.0
    g2 = 1.0
    while True:
        done = False
        while not done:  
            try:
                g3 = g2 - ((g2**b) - n) / (b * (g2**(b-1)))
                done = True
            except OverflowError:
                g2 = (g1 + g2) / 2.0 

        if abs(g2 - g3) < 1.0 / (10**sig_figs):
            return g3
        g1 = g2
        g2 = g3

def npowbysqr(n, p):
    if p == 0:
        return 1.0
    if p % 2 == 0:
        v = npowbysqr(n, p/2)
        return v*v 
    else:
        return n*npowbysqr(n, p-1)

def npow(n, p):
    return npowbysqr(nroot(n, 1000000), int(p*1000000))

print npow(5, 4.3467)
print 5**4.3467   
         

I should add that there are probably much better solutions. This does seem to work, however

Answer 2

I happened to need something like this a while ago. Thankfully, Dave DeLong had been tinkering with this in his DDMathParser, so I built off of that. He yanked his implementation from his code in this commit, but I took that and modified it. This is my version of his NSDecimal power function:

extern NSDecimal DDDecimalPower(NSDecimal d, NSDecimal power) {
    NSDecimal r = DDDecimalOne();
    NSDecimal zero = DDDecimalZero();
    NSComparisonResult compareToZero = NSDecimalCompare(&zero, &power);
    if (compareToZero == NSOrderedSame) {
        return r;
    }
    if (DDDecimalIsInteger(power))
    {
        if (compareToZero == NSOrderedAscending)
        {
            // we can only use the NSDecimal function for positive integers
            NSUInteger p = DDUIntegerFromDecimal(power);
            NSDecimalPower(&r, &d, p, NSRoundBankers);
        }
        else
        {
            // For negative integers, we can take the inverse of the positive root
            NSUInteger p = DDUIntegerFromDecimal(power);
            p = -p;
            NSDecimalPower(&r, &d, p, NSRoundBankers);
            r = DDDecimalInverse(r);
        }
    } else {
        // Check whether this is the inverse of an integer        
        NSDecimal inversePower = DDDecimalInverse(power);
        NSDecimalRound(&inversePower, &inversePower, 34, NSRoundBankers); // Round to 34 digits to deal with cases like 1/3

        if (DDDecimalIsInteger(inversePower))
        {
            r = DDDecimalNthRoot(d, inversePower);
        }
        else
        {
            double base = DDDoubleFromDecimal(d);
            double p = DDDoubleFromDecimal(power);
            double result = pow(base, p);
            r = DDDecimalFromDouble(result);
        }        
    }
    return r;
}

It tries to identify common cases and use more precise calculations for those. It does fall back on pow() for things that don't fit in these cases, though.

The rest of the NSDecimal functions I use can be found here and here.

Answer 3

I have come up with a function that suits my needs and will hopefully suit the needs of many others. the following method is fully annotated and works for any power function that has a real value. This method also only uses NSDecimalNumbers meaning you will not loose any precision due to float rounding error. This method takes two arguments one for the base and one for the power, and both are NSDecimalNumbers. So here it is:

//these are constants that will be used
NSDecimalNumber *ten = [NSDecimalNumber decimalNumberWithString:@"10"];
NSDecimalNumber *one = NSDecimalNumber.one;

//these will together hold the power in fractional form
NSDecimalNumber *numerator = power, *denominator = one;

//this will hold the final answer and all previous guesses the first guess is set to be the base
NSDecimalNumber *powAns = base;

//this will hold the change in your guess, also serves as an idea of how large the error is
NSDecimalNumber *error = one;

//part1 holds f(x) and part2 holds f'(x)
NSDecimalNumber *part1, *part2;

//if the base is < 0 and the power is not whole, answer is not real
if ([base doubleValue] < 0 && [[power stringValue] rangeOfString:@"."].location != NSNotFound)
    return NSDecimalNumber.notANumber;

//converts power to a fractional value
while ([[numerator stringValue] rangeOfString:@"."].location != NSNotFound) {

    numerator = [numerator decimalNumberByMultiplyingBy:ten];
    denominator = [denominator decimalNumberByMultiplyingBy:ten];
}

//conditions here are the precision you wish to get
while ([error compare:[NSDecimalNumber decimalNumberWithString:@"1e-20"]] == NSOrderedDescending ||
       [error compare:[NSDecimalNumber decimalNumberWithString:@"-1e-20"]] == NSOrderedAscending) {

    //if this catches an overflow error it is set to be a very large number otherwise the value cannot be a number, however no other error should be returned.
    @try {
        part1 = [powAns decimalNumberByRaisingToPower:[denominator intValue]];
    }
    @catch (NSException *exception) {

        if ([exception.name isEqual: NSDecimalNumberOverflowException])
            part1 = [NSDecimalNumber decimalNumberWithString:@"10e127"];

        else
            return NSDecimalNumber.notANumber;
    }

    part1 = [part1 decimalNumberBySubtracting:base];

    //if this catches an overflow error it is set to be a very large number otherwise the value cannot be a number, however no other error should be returned.
    @try {
        part2 = [powAns decimalNumberByRaisingToPower:[denominator intValue]-1];
        part2 = [part2 decimalNumberByMultiplyingBy:denominator];
    }
    @catch (NSException *exception) {

        if ([exception.name isEqual: NSDecimalNumberOverflowException])
            part2 = [NSDecimalNumber decimalNumberWithString:@"10e127"];

        else
            return NSDecimalNumber.notANumber;
    }

    //error is the change in the estimated value or y - f(x)/f'(x)
    error = [part1 decimalNumberByDividingBy:part2];
    powAns = [powAns decimalNumberBySubtracting: error];
}

//if the numerator value is negative it must be made positive and the answer is then inverted
if ([numerator intValue] < 0) {
    powAns = [powAns decimalNumberByRaisingToPower:abs([numerator intValue])];
    powAns = [one decimalNumberByDividingBy:powAns];
}
else
    powAns = [powAns decimalNumberByRaisingToPower:[numerator intValue]];

return powAns;

If anyone has any questions about my code I am happy to answer them.