Floating Point Power Calculation (PHP / BCMath)

Question

I am writing a wrapper for the bcmath extension and error # 10116 regarding bcpow() particularly annoying - it distinguishes $right_operand ( $exp ) from integer (native PHP, not arbitrary length), so when you try to calculate the square root (or any other root above than 1 ), the number you always end with 1 instead of the correct result.

I started looking for algorithms that would allow me to calculate the nth root from a number, and I found this answer that looks pretty solid, I actually expanded the formula with WolframAlpha, and I was able to improve its speed by about 5%, keeping with the accuracy of the results.

Here is a clean PHP implementation that mimics 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 works fine , except when 1 / fmod($exp, 1) does not give an integer . For example, if $exp is 0.123456 , its inverse will be 8.10005 , and the result of pow() and _pow() will be slightly 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 algorithm used to find the root of n ^th from the (positive) number a is Newton's algorithm for finding zero

 f(x) = x^n - a. 

This includes only forces with natural numbers as indicators, so it’s easy to implement them.

Calculation of power with the exponent 0 < y < 1 , where y does not have the form 1/n with integer n , is more complicated. Performing an analog, solution

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

will again include power calculation with a non-integer metric, the very problem we are trying to solve.

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

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

but for most rational ones 0 < y < 1 numerator and denominator are quite large, and the intermediate x^n will be huge, so the calculation will use a lot of memory and take a (relatively) long time. (For an approximate figure of 0.123456 = 1929/15625 this is not so bad, but 0.1234567 will be taxable.)

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

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

with integers a < b < c < ... < q and multiply / divide the individual x^(1/k) . (Each rational 0 < y < 1 has such representations, and the shortest such representations, as a rule, do not contain many terms, for example.

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

using only additions in the decomposition leads to longer representations with large denominators, for example

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

which will require more work.)

Some improvement is possible by mixing approaches, first find a close rational approximation to y with a small denominator by continuing the fractional fraction of y - for the approximate indicator 1929/15625 = [0;8,9,1,192] and using the first four partial fractions gives an approximation of 10/81 = 0.123456790123... [note that 10/81 = 1/8 - 1/648 , the partial sums of the shortest decomposition into pure fractions are convergent), and then we decompose the remainder into pure fractional parts.

However, generally speaking, this approach leads to the calculation of the roots of n ^th for large n , which is also slow and memory intensive if the desired accuracy of the final result is high.

In general, it may be easier and faster to implement exp and log and use

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