Pow implementation for double

Question

I am developing a code that will be used for motion control and I am having a issue with the pow function. I am using VS2010 as IDE.

Here is my issue: I have:

    double p = 100.0000;
    double d = 1000.0000;
    t1 = pow((p/(8.0000*d),1.00/4.000);

When evaluating this last function, I don't get the better approximation as result. I am getting a 7 decimal digits correct, and the consequent digits are all trash. I am guessing that pow function only casts any input variable as float and proceds with calculation.

Am I right?
If so, is there any code I can get "inspired" with to reimplement pow for a better precision?

Edit: Solved.

After all, I was having problems with FPU config bits, caused by Direct3D which was being used by OGRE 3D framework.

If using OGRE, on the config GUI, just set "Floating-point mode=Consistent".

If using raw Direct3D, when calling CreateDevice, make sure to pass "D3DCREATE_FPU_PRESERVE" flag to it.

Original post:

You may be using a libray that is changing the default precision of the FPU to single-precision. Then all floating-point operations, even on doubles, will actually be performed as single-precision operations.

As a test, you can try calling _controlfp( _CW_DEFAULT, 0xfffff ); (you need to include ) before performing the calculation to see if you get the correct result. This will reset the floating-point control word to default values. Note that it will reset other settings as well, which may cause issues.

One common library that changes the floating-point precision is Direct3D 9 (and maybe other versions too): By default, it changes the FPU to single-precision when creating a device. If you use it, specify the flag D3DCREATE_FPU_PRESERVE when creating the device to prevent it from changing the FPU precision.

see the prototype of `pow()` in the standard library comes with your compiler. glibc's `pow()` has a prototype with all `double`. The `float` version is called `powf()`. — Aftnix, Jul 17 '12 at 14:25
C99 has the `long double` type for extra precision, but apparantly [Microsoft decided to get rid that](http://msdn.microsoft.com/en-us/library/9cx8xs15.aspx). MS C++ does have `pow` for `long double` arguments. — Fred Foo, Jul 17 '12 at 14:40
gcc produced this with "50 digits" precision 0.33437015248821100321663379872916266322135925292969 — Aftnix, Jul 17 '12 at 14:44
I did... I am using in deed double pow(double, double)... But, I mean, don't you find it very wierd that the function returns exactly 7 decimal digits precise results, the exact number of decimal digits that float type can bear? — Paulo Sherring, Jul 17 '12 at 14:57
Excel gives (1/80)^(1/4) == 0.334370152488211, in agreement with both what Aftnix got with gcc and what I got using VS2005. Also, your example doesn't paste well: prefix `t1` with double and fix the unbalanced parenthesis. — Prashant Kumar, Jul 17 '12 at 15:14

score 4 · Accepted Answer · answered Jul 17 '12 at 15:45

4

You may be using a libray that is changing the default precision of the FPU to single-precision. Then all floating-point operations, even on doubles, will actually be performed as single-precision operations.

As a test, you can try calling _controlfp( _CW_DEFAULT, 0xfffff ); (you need to include <float.h>) before performing the calculation to see if you get the correct result. This will reset the floating-point control word to default values. Note that it will reset other settings as well, which may cause issues.

One common library that changes the floating-point precision is Direct3D 9 (and maybe other versions too): By default, it changes the FPU to single-precision when creating a device. If you use it, specify the flag D3DCREATE_FPU_PRESERVE when creating the device to prevent it from changing the FPU precision.

answered Jul 17 '12 at 15:45

interjay

107,303
21
270
254

Interjay: Thank you sooooo much! I am using D3D and finally got rid of this terrible problem! – Paulo Sherring Jul 18 '12 at 00:22
Other thing: I am using OGRE 3D for rendering. I thought of reading FPU status, changing it, then, afterr calculations, change it back. But this might not be thread safe, since cpu may be preempted before changing it back. What is the way to go: trying to do it through OGRE or trying some critical section thingy? Thank you so much, again and again! – Paulo Sherring Jul 18 '12 at 00:36
Ok, just made it through configuration GUI of direct3d rendering system. Thanks!!! – Paulo Sherring Jul 18 '12 at 02:52

score 0 · Answer 2 · answered Jul 17 '12 at 15:31

0

How did you determine you're only getting 7 digits of precision? Are you printing t1 and specifying the correct output format? On my machine, with VS2010, the following code:

int main()
{
    double p = 100.0000;
    double d = 1000.0000;
    double t1 = pow(p/(8.0000*d),1.00/4.000);

    printf("t1=%.15f\n", t1); // C
    std::cout << "t1=" << std::setprecision(15) << t1 << '\n'; // C++
}

Produces this output:

t1=0.334370152488211
t1=0.334370152488211

answered Jul 17 '12 at 15:31

Blastfurnace

18,411
56
55
70

I could identify it by analysing the result. I was getting 7 digits right, the other digits all diverged from the expected results. I solved this by using interjay's hint. – Paulo Sherring Jul 18 '12 at 00:32

Pow implementation for double

2 Answers2

Linked