Estimation Off

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.57.6 = 1057.65509515..., but my function gives 89321.72336....

Here it is:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
/*   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 ex 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;
    else if(exp == 1)   return *this;
    else if(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;
}


Appreciate any help.
Last edited on
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 e6.963809562

Using an online compiler and directly summing the power series for ex 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 ex your answer would be too small since all the terms in the power series for ex 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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#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;
}


This yields 1057.66
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.57.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....

Here is the new code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
Impl_Float_& Impl_Float_::operator^=(const Impl_Float_& exp){
    if(exp == 0)        return *this = 1;
    else if(exp == 1)   return *this;
    else if(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.

I am really braindead over this problem.
Last edited on
Update:

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?
1.) Mathematically, a negative base raised to a non-integer power is undefined. Note the conditions of x and y in the power function in <cmath>:

http://www.cplusplus.com/reference/cmath/pow/

double pow (double base, double exponent);

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

http://www.amazon.com/Elementary-Functions-Prentice-Hall-computational-mathematics/dp/0138220646/ref=sr_1_1?ie=UTF8&qid=1378258040&sr=8-1&keywords=software+manual+for+the+elementary+functions

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.

Sorry I couldn't be of more help.

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 = 11/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.57.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.1283749182736468128200837.99022110201983

And I have searched how to approximate ln(x). I even tried to find methods that would work better than eb*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.
Last edited on
(-1)7.6 = 1
If this was true, the base 7.6 logarithm of 1 would be -1, however log_k(1) = 0 for all positive k.

(-1)^(38/5) is approx. 0.31 - 0.95 i (if we assume that x^y = |x| |y| e^(arg(x) arg(y) i)).
Last edited on
Oh, thanks for the clarification. Didn't think about it that way.
Durr. I used the wrong operation in my argument. You are indeed correct.
Topic archived. No new replies allowed.