I've been implementing my own floating-point class, but I got stuck on writing an efficient, quick, and accurate exponentiation function.
The problems are that the function takes up a lot of runtime for two calls (like 1 - 2 seconds) and the answer is way off. For example, 2.5^{7.6} = 1057.65509515..., but my function gives 89321.72336....
/* Some notes to answer some questions I know you'll ask.
- I use other operators in this function, and those have already been tested.
However, if you want to see their source codes, I'll be happy to provide
them. It would not be the first time I mistook the root of the problem.
- Most of you are probably going to frown upon me overloading the XOR operator
for this, but I reason that while XOR makes no sense for FP, it is a valid
exponentiation symbol in mathematics.
- My algorithm uses the Maclaurin series for e^{x} and for ln(1+x).
- I use another overloaded XOR operator, but that one will call
Impl_Float_::operator^=(const Integer&), which has also been tested and
uses the exponentiation by squaring algorithm.
- m_precision is a size_t variable holding the maximum number of digits to
the right of the decimal point.
*/
Impl_Float_& Impl_Float_::operator^=(const Impl_Float_& exp){
if(exp == 0) return *this = 1;
elseif(exp == 1) return *this;
elseif(exp == -1) return *this = 1 / *this;
/*
a^b [not XOR]
= e^(ln(a^b))
= e^(b*ln(a))
ln(a) = SUM[n:1->INFINITY]( ( (-1^(n+1)) * (a-1)^n )/n )
e^C = SUM[i:0->INFINITY]( C^i/i! )
*/
Impl_Float_ nlogged(0);
for(size_t i(1); i < m_precision; ++i){
nlogged
+= (i % 2 == 1 ? 1 : -1)
* (*this - 1) ^ i
/ i
;
}
nlogged *= exp;
*this = 0;
for(size_t i(0); i < m_precision; ++i)
*this += (nlogged ^ i) / Precision::Factorial(i);
return *this;
}
Without really digging too deeply into your code, your b*ln(a) would be:
2.5 * ln( 7.6 ) = 6.963809562 so you would be calculating e^{6.963809562}
Using an online compiler and directly summing the power series for e^{x} adding the first 40 terms yields 1057.66 which is fairly accurate. The fact that your result is bigger than the correct value probably means that your calculation for ln(x) is not accurate enough or has an error. If you were not taking enough terms in calculating e^{x} your answer would be too small since all the terms in the power series for e^{x} are positive.
Try watching in the debugger to see how accurate the b*ln(a) value is and if it is significantly off, either increase the precision or look for an error in the implementation of the ln function.
#include <iostream>
#include <cmath>
double fact( int n )
{
double product = 1;
for ( int i = 1; i < n + 1; ++ i)
{
product *= i;
}
return product;
}
int main()
{
double sum = 0;
for ( int i = 0; i < 40; ++i )
{
sum += pow( 6.963809562, i ) / fact( i );
}
std::cout << sum << std::endl;
return 0;
}
Edit:
Disregard this post. I completely forgot how to find the Taylor series._{I am so sorry my dear 'ol Calculus teacher...}
You're right. My estimation of the natural log was horrid. My nlogged variable absolutely exploded and returned a huge number. To counteract this, I switched the Maclaurin series to Taylor. The center is at the Integer cast of the base, which should ensure that the accuracy is focused within a (base - 1, base + 1) range. And indeed, it did help, but now the function is underestimating.
For 2.5^{7.6}, ln(2.5) became 0.559615... and 7.6*ln(2.5) became 4.2530799...., both of which are way off. The function returned 70.32166906....
Impl_Float_& Impl_Float_::operator^=(const Impl_Float_& exp){
if(exp == 0) return *this = 1;
elseif(exp == 1) return *this;
elseif(exp == -1) return *this = 1 / *this;
/*
a^b [not XOR]
= e^(ln(a^b))
= e^(b*ln(a))
ln(a) = SUM[n:1->INFINITY](
((-1^(n+1))
* (a-1)^n
* (a - int(a))^n
/ n
)
e^C = SUM[i:1->INFINITY]( C^i/i! )
*/
Impl_Float_ nlogged(0);
for(size_t i(1); i < m_precision; ++i){
nlogged
+= (i % 2 == 1 ? 1 : -1)
* ( (*this - 1) ^ i )
* ( (*this - this->integer()) ^ i ) //Center the series near the base
/ i
;
}
nlogged *= exp;
*this = 0;
for(size_t i(0); i < m_precision; ++i)
*this += (nlogged ^ i) / Precision::Factorial(i);
return *this;
}
Also, the number of terms is based on the maximum number of digits to the right of the decimal--m_precision. And by default, m_precision is assigned 100, so there should be plenty of terms to make the series as accurate as possible. I even messed with the number terms--doubling it, halving it--but it only affected the last few digits out of 100.
Okay, so my function actually is accurate (though still slow) for a range of (0, 1.8]. I can relate this to how series estimations blow up outside of a (x-1, x+1) range.
So what is an efficient way to increase accuracy of my series for large bases? It's already using 100 terms (by default).
And is there a way to center the radius of convergence (for the Taylor expansion of ln(x)) near the value of the base?
Also, how should I handle negative bases since the exponents are not whole? Would taking the exponent, converting it to a fraction, and seeing if the denominator is odd be sufficient?
The result of raising base to the power exponent.
If x is finite negative and y is finite but not an integer value, it causes a domain error.
2.) Implementing either e^{x} or ln(x) to achieve efficient and accurate results is a non-trivial problem. I took a look at the book The Standard C Library by P. J. Plauger, 1992 which contains code and comments with an implementation of the old C library "math.h". This code is fairly complicated with scaling of the x parameter, approximations using the ratio of certain polynomials and numerous error checks (underflow, overflow etc.). His comment about the power function pow() is concise:
The function pow, which raises x to the y power is easily the most complex of all the math functions.
Trying to approximate the power series is simply not going to work as you have found.
Here is a reference that he uses and borrows from extensively:
Software Manual for the Elementary Functions by William J. Cody, Jr and William Waite, 1980
Since your primary goal is to overload operator ^ (we wont comment on that) :) why not bypass this complicated problem and just call pow() from <cmath>? Or google "implementing ln(x)" and see what you come up with.
1) I could argue that if I took the exponent--such as 7.6--and found the denominator for the most simplified fraction--5 in this case--I can tell if the answer will be real or complex. In fact, since I am dealing with base10, all I'd have to do is check if the right most, non-zero digit (if the exponent is a FP) is even or odd. Then, to see if -1 will be squared away, I could check if the numerator is even
So for (-2.5)^{7.6} I separate it to (-1)^{7.6}*(2.5)^{7.6}
-->7.6 -->76
-->Even? Yes
-->76/2 = 38
-->Still even? Yes
-->(-1)^{7.6} = 1
Which makes sense because (-1)^{7.6} = (-1)^{38/5} = ( (-1)^{38} )^5^{-1} = 1^{1/5} = 1.
The test would also show that (-1)^{0.6} = -1.
Or would this be too inefficient?
2) That really is a shame that pow(double, double) requires a complex algorithm... I always thought the standard functions wouldn't be too complex.
I wish I could just call <cmath> functions, but this type will deal with numbers too large to store in the primitive types. 2.5^{7.6} was just a simple test to see if my algorithm was working.
For example, my type would have to be able to handle: 198237469.1283749182736468^{128200837.99022110201983}
And I have searched how to approximate ln(x). I even tried to find methods that would work better than e^{b*ln(a)}. I haven't really found a good alternative... Though there are plenty of algorithms for integer exponents.
Thank you for the insight though. I guess I'll have to leave my exponentiation function unfinished for the time being.