PHP: How to raise number to (tiny) fractional exponent?

Question

I'm doing a calculation in PHP using bcmath, and need to raise e by a fractional exponent. Unfortunately, bcpow() only accepts integer exponents. The exponent is typically higher precision than a float will allow, so normal arithmetic functions won't cut it.

For example:

$e = exp(1);
$pow = "0.000000000000000000108420217248550443400745280086994171142578125";
$result = bcpow($e, $pow);

Result is "1" with the error, "bc math warning: non-zero scale in exponent".

Is there another function I can use instead of bcpow()?

Answer 1

Your best bet is probably to use the Taylor series expansion. As you noted, PHP's bcpow is limited to raising to integer exponentiation.

So what you can do is roll your own bc factorial function and use the wiki page to implement a Taylor series expansion of the exponential function.

function bcfac($num) { 
    if ($num==0) return 1;
    $result = '1';
    for ( ; $num > 0; $num--)
        $result = bcmul($result,$num);
    return $result;
}
$mysum = '0';
for ($i=0; $i<300; $i++) {
    $mysum = bcadd($mysum, bcdiv(bcpow($pow,$i), bcfac($i)) );
}
print $mysum;

Obviously, the $i<300 is an approximation for infinity... You can change it to suit your performance needs.

With $i=20, I got

1.00000000000000000010842021724855044340662275184110560868263421994092888869270293594926619547803962155136242752708629105688492780863293090291376157887898519458498571566021915144483905034693109606778068801680332504212458366799913406541920812216634834265692913062346724688397654924947370526356787052264726969653983148004800229537555582281617497990286595977830803702329470381960270717424849203303593850108090101578510305396615293917807977774686848422213799049363135722460179809890014584148659937665374616

This is comforting since that small of an exponent should yield something really close to 1.0.

Answer 2

Old question, but people might still be interested nonetheless.

So Kevin got the right idea with the Taylor-polynomial, but when you derive your algorithm from it directly, you can get into trouble, mainly your code gets slow for long input-strings when using large cut-off values for $i.

Here is why: At every step, by which I mean with each new $i, the code calls bcfac($i). Everytime bcfac is called it performs $i-1 calculations. And $i goes all the way up to 299... that's almost 45000 operations! Not your quick'n'easy floating point operations, but slow BC-string-operations - if you set bcscale(100) your bcmul has to handle up to 10000 pairs of chars!

Also bcpow slows down with increasing $i, too. Not as much as bcfac, because it propably uses something akin to the square-and-multiply method, but it still adds something.

Overall the time required grows quadraticly with the number of polynomial terms computed.

So... what to do?

Here's a tip:

Whenever you handle polynomials, especially Taylor-polynomials, use the Horner method.

It converts this: exp(x) = x^0/0! + x^1/1! + x^2/2! + x^3/3! + ...

...into that: exp(x) = ((( ... )*x/3+1 )*x/2+1 )*x/1+1

And suddenly you don't need any powers or factorials at all!

function bc_exp($number) {
    $result = 1;
    for ($i=299; $i>0; $i--)
        $result = bcadd(bcmul(bcdiv($result, $i), $number), 1);
    return $result;
}

This needs only 3 bc-operations for each step, no matter what $i is. With a starting value of $i=299 (to calculate exp with the same precision as kevin's code does) we now only need 897 bc-operations, compared to more than 45000. Even using 30 as cut-off instead of 300, we now only need 87 bc-operations while the other code still needs 822 for the factorials alone.

Horner's Method saving the day again!

Some other thoughts:

1) Kevin's code would propably crash with input="0", depending on how bcmath handles errors, because the code trys bcpow(0,0) at the first step ($i=0).

2) Larger exponents require longer polynomials and therefore more iterations, e.g. bc_exp(300) will give a wrong answer, even with $i=299, whyle something like bc_exp(3) will work fine and dandy. Each term adds x^n/n! to the result, so this term has to get small before the polynomial can start to converge. Now compare two consecutive terms:

 ( x^(n+1)/(n+1)! ) / ( x^n/n! ) = x/n

Each summand is larger than the one before by a factor of x/n (which we used via the Horner method), so in order for x^(n+1)/(n+1)! to get small x/n has to get small as well, which is only the case when n>x.

Inconclusio: As long as the number of iterations is smaller than the input value, the result will diverge. Only when you add steps until your number of iterations gets larger than the input, the algorithm starts to slowly converge.

In order to reach results that can satisfie someone who is willing to use bcmath, your $i needs to be significantly larger then your $number. And that's a huge proplem when you try to calculate stuff like e^346674567801

A solution is to divide the input into its integer part and its fraction part. Than use bcpow on the integer part and bc_exp on the fraction part, which now converges from the get-go since the fraction part is smaller than 1. In the end multiply the results.

e^x = e^(intpart+fracpart) = e^intpart * e^fracpart = bcpow(e,intpart) * bc_exp(fracpart)

You could even implement it directly into the code above:

function bc_exp2($number) {
    $parts = explode (".", $number);
    $fracpart = "0.".$parts[1];
    $result = 1;
    for ($i=299; $i>0; $i--)
        $result = bcadd(bcmul(bcdiv($result, $i), $fracpart), 1);
    $result = bcmul(bcpow(exp(1), $parts[0]), $result);
    return $result;
}

Note that exp(1) gives you a floating-point number which propably won't satisfy your needs as a bcmath user. You might want to use a value for e that is more accurate, in accordance with your bcscale setting.

3) Talking about numbers of iterations: 300 will be overkill in most situations while in some others it might not even be enough. An algorithm that takes your bcscale and $number and calculates the number of required iterations would be nice. Alraedy got some ideas involving log(n!), but nothing concrete yet.

4) To use this method with an arbitrary base you can use a^x = e^(x*ln(a)). You might want to divide x into its intpart and fracpart before using bc_exp (instead of doing that within bc_exp2) to avoid unneccessary function calls.

function bc_pow2($base,$exponent) {
    $parts = explode (".", $exponent);
    if ($parts[1] == 0){
        $result = bcpow($base,$parts[0]);
    else $result = bcmul(bc_exp(bcmul(bc_ln($base), "0.".$parts[1]), bcpow($base,$parts[0]);
    return result;
}

Now we only need to program bc_ln. We can use the same strategy as above:

Take the Taylor-polynomial of the natural logarithm function. (since ln(0) isn't defined, take 1 as developement point instead) Use Horner's method to drasticly improve performance. Turn the result into a loop of bc-operations. Also make use of ln(x) = -ln(1/x) when handling x > 1, to guarantee convergence.