I was writing a program to calculate sin(x) with x in [0, 1], using Taylor series that include raising some of its terms to a certain power. Since I'm new to C, I decided to write my own pow() function just to practice. But then I noticed that my program doesn't calculate the sine function correctly. Upon an examination, I discovered that raising a number to a power in a loop doesn't produce the same result as simple straightforward multiplication ...v * v * v...

Here's my powr() function and difference in results:

#include <stdio.h>    

double powr(double base, int power)
{
    if (power == 0) {
        base = 1;
    } else if (power > 0) {
        for (int i = 1; i < power; i  ) {
            base *= base;
        }
    } else {
        // power is negative.
    }
    printf("In loop number to a power of 3: %lf\n", base);
    double hard_coded = 0.99 * 0.99 * 0.99;
    printf("Out of loop number to a power of 3: %lf\n", hard_coded); 

    return base;
}

int main(void)
{
    double num = 0.99;
    double num_to_power = powr(num, 3);
    return 0;
}

Output:

In loop number to a power of 3: 0.960596
Out of loop number to a power of 3: 0.970299

0.960596 and 0.970299 may not look that different, but with smaller values an error becomes even bigger. Why is that and how can I solve it?

Also, I just noticed that if I take small values like 0.123 and raise them to a power of 6, for example. I have an output like 0.000000. I understand that it's just too small to show up. But how can I print an a floating-point number that has more than 6 after point numbers? I've tries long double, didn't work either.

CodePudding user response：

Your loop calculation is wrong. The base *= base; line will, on each loop, multiply the 'accumulated' value by the value that has been modified on the previous loop, rather than by the original value (as it should).

You need to store the passed base value and use that saved value as the multiplier on each loop:

double powr(double base, int power)
{
    double saved = base;
    if (power == 0) {
        base = 1;
    }
    else if (power > 0) {
        for (int i = 1; i < power; i  ) {
            base *= saved;
        }
    }
    else {
        // power is negative.
    }

    printf("In loop number to a power of 3: %lf\n", base);
    double hard_coded = 0.99 * 0.99 * 0.99;
    printf("Out of loop number to a power of 3: %lf\n", hard_coded);

    return base;
}

The two output vales will then be equivalent:

In loop number to a power of 3: 0.970299
Out of loop number to a power of 3: 0.970299

CodePudding user response：

First, your loop does not calculate what you need:

for (int i = 1; i < power; i  ) {
    base *= base;
}

In you case you are calculating:

1. first iteration base *= base will give base=0.99*0.99=0.9801
2. second iteration base *= base will give base=0.9801*09801=0.96059601

By accumulating previous result over time you are calculating base^(2^(exp-1)). No surprise then. For exponentiation a simple loop like this will make it:

result = base;
for (int i = 1; i < power; i  ) {
    result *= base;
}

There exists a simple fast exponentiation that use squaring if you want to make it faster (https://en.wikipedia.org/wiki/Exponentiation_by_squaring).

Second, the hard coded version can be computed by the compiler in any other way, and as floating point computation is very sensible to rounding errors, the result may be different.