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 * _…

Accepted Answer

The employed algorithm to find the nth 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 nth 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))
```

Calculating Floating Point Powers (PHP/BCMath)

Solution

Problem

Related problems