Calculating Floating Point Powers (PHP/BCMath)

Question

I'm writing a wrapper for the bcmath extension, and bug #10116 regarding bcpow() is particularly annoying -- it casts the $right_operand ($exp) to an (native PHP, not arbitrary length) integer, so when you try to calculate the square root (or any other root higher than 1) of a number you always end up with 1 instead of the correct result.

I started searching for algorithms that would allow me to calculate the nth root of a number and I found this answer which looks pretty solid, I actually expanded the formula using WolframAlpha and I was able to improve it's speed by about 5% while keeping the accuracy of the results.

Here is a pure PHP implementation mimicking my BCMath implementation and its limitations:

function _pow($n, $exp)
{
    $result = pow($n, intval($exp)); // bcmath casts $exp to (int)

    if (fmod($exp, 1) > 0) // does $exp have a fracional part higher than 0?
    {
        $exp = 1 / fmod($exp, 1); // convert the modulo into a root (2.5 -> 1 / 0.5 = 2)

        $x = 1;
        $y = (($n * _pow($x, 1 - $exp)) / $exp) - ($x / $exp) + $x;

        do
        {
            $x = $y;
            $y = (($n * _pow($x, 1 - $exp)) / $exp) - ($x / $exp) + $x;
        } while ($x > $y);

        return $result * $x; // 4^2.5 = 4^2 * 4^0.5 = 16 * 2 = 32
    }

    return $result;
}

The above seems to work great except when 1 / fmod($exp, 1) doesn't yield an integer. For example, if $exp is 0.123456, its inverse will be 8.10005 and the outcome of pow() and _pow() will be a bit different (demo):

• pow(2, 0.123456) = 1.0893412745953
• _pow(2, 0.123456) = 1.0905077326653
• _pow(2, 1 / 8) = _pow(2, 0.125) = 1.0905077326653

How can I achieve the same level of accuracy using "manual" exponential calculations?

Answer 1

The employed algorithm to find the n^th root of a (positive) number a is the Newton algorithm for finding the zero of

f(x) = x^n - a.

That involves only powers with natural numbers as exponents, hence is straightforward to implement.

Calculating a power with an exponent 0 < y < 1 where y is not of the form 1/n with an integer n is more complicated. Doing the analogue, solving

x^(1/y) - a == 0

would again involve calculating a power with non-integral exponent, the very problem we're trying to solve.

If y = n/d is rational with small denominator d, the problem is easily solved by calculating

x^(n/d) = (x^n)^(1/d),

but for most rational 0 < y < 1, numerator and denominator are rather large, and the intermediate x^n would be huge, so the computation would use a lot of memory and take a (relatively) long time. (For the example exponent of 0.123456 = 1929/15625, it's not too bad, but 0.1234567 would be rather taxing.)

One way to calculate the power for general rational 0 < y < 1 is to write

y = 1/a ± 1/b ± 1/c ± ... ± 1/q

with integers a < b < c < ... < q and to multiply/divide the individual x^(1/k). (Every rational 0 < y < 1 has such representations, and the shortest such representations generally don't involve many terms, e.g.

1929/15625 = 1/8 - 1/648 - 1/1265625;

using only additions in the decomposition leads to longer representations with larger denominators, e.g.

1929/15625 = 1/9 + 1/82 + 1/6678 + 1/46501020 + 1/2210396922562500,

so that would involve more work.)

Some improvement is possible by mixing the approaches, first find a close rational approximation to y with small denominator via the continued fraction expansion of y - for the example exponent 1929/15625 = [0;8,9,1,192] and using the first four partial quotients yields the approximation 10/81 = 0.123456790123... [note that 10/81 = 1/8 - 1/648, the partial sums of the shortest decomposition into pure fractions are convergents] - and then decompose the remainder into pure fractions.

However, in general that approach leads to calculating n^th roots for large n, which also is slow and memory-intensive if the desired accuracy of the final result is high.

All in all, it is probably simpler and faster to implement exp and log and use

x^y = exp(y*log(x))