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?

Thyroiditis

It's working exactly as advertised. _pow 'rounds' the fractional part to the nearest 1/n. You could make this work recursively. So after computing _pow(2, 0.125), you calculate _pow(2,0.125-123456) and so on. – Thermobarometer

Ah, now I understand. So doesn't bcmath have exp and log or are there other reasons why a^b = exp(b*log(a)) isn't an option? The recursion Jeffrey suggests would of course work, but its speed may not be satisfactory if you need many 1/k to represent the exponent. Is writing the exponent as a rational number n/d and calculating (a^n)^(1/d) an option, or must too large n and d be expected? Perhaps worth an investigation is approximating the exponent by a rational number with smallish denominator (continued fraction expansion) and doing the rest with recursion. – Priapus

@JeffreySax: Ah, I see... That's a bummer but still doesn't seem to work (codepad.org/eI4ykyQU) or am I missing something? – Thyroiditis

@DanielFischer: Thanks for getting back to me! =) Well, the bcmath API is pretty poor, besides */+- we have sqrt and a crippled pow: php.net/manual/en/ref.bc.php. One problem I see with calculating (a^n)^(1/d) is that 1/d might also be a irrational number. Either way, I asked this mostly because I was curious -- I doubt I'll need to use irrational exponents on such big numbers. =) – Thyroiditis

I think we can safely ignore irrational numbers. We can approximate them arbitrarily well with rational numbers. The problem is that the numerator and denominator of such an approximation may be huge. Can you specify what sort of input you want to treat and what accuracy you want in the result? The fewer digits you need, the smaller numerators and denominators you can get away with in the approximations. – Priapus

@DanielFischer: Actually, I just needed to be able to calculate roots, and your answer solved that problem perfectly. =) This question was mostly motivated because some of my tests (ideone.com/KNztd) weren't passing and I was curious if this could be done in a straightforward way. Since not, and since this was mostly to satisfy my own curiosity I'm voting to close, feel free to post your fist comment as an answer and I'll accept it. =) Thanks for your time BTW. – Thyroiditis

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))

Priapus

Great, detailed answer! Thank you. – Thyroiditis