# C++ : Will an (affine transformation...no) better pow() function improve accuracy?

This topic is 4903 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

This thread previously considered affine transformations to improve accuracy. It seems the answer is no, and the best solution appears to be to use rational/bigint data types (Thanks Fruny/ZQJ). However, what about an improved pow() function instead...this seems to be the most typical point for truncation 'error' accumulation. Any ideas? ==========================Previous thread below========= ======================= I've tested this equation on several platforms and in several languages: f = 333.75 * y^6 + x^2 * ( 11 * x^2 * y^2 - y^6 - 121 * y^4 - 2 ) + 5.5 * y^8 + x/(2 * y) References: Stolfi and Figueiredo, "Self Validated Numerical Methods and Applications", May 1997, P 3, and Rump, "Algorithms for verified inclusions: theory and practice", Perspectives in Computing, 1988. It was originally crafted by Rump to illustrate the danger of double precision in modern computers. The correct answer is: f = -0.827396... This can be achieved in Lisp by rational type. A simple lisp code is:
(defun f (a b x y)
(+ (* a y y y y y y)
(* x x (- (* 11 x x y y) (* y y y y y y) (* 121 y y y y) 2))
(* b y y y y y y y y)
(/ x (* 2 y))))

;; Values entered in rational form.
(defun rump-stolfi ()
(f 33375/100 55/10 77617 33096))

;; Take float as last step.
(float (rump-stolfi))


Which will yield the correct result in 32 bit CLisp. However, this same equation in C++ on a 32 bit machine, as shown below:
#include <cstdlib>
#include <cstdio>
#include <iostream>

double f(double, double, double, double);

int main (int argc, const char * argv[])
{
std::cout << "\n" << f(333.75,5.5,77617,33096) << std::endl;
return EXIT_SUCCESS;
}

double f(double a, double b, double x, double y)
{
// f = 333.75 * y^6
//     + x^2 * ( 11 * x^2 * y^2 - y^6 - 121 * y^4 - 2 )
//     + 5.5 * y^8
//     + x / ( 2 * y )
return ( a * pow(y,6))
+ ( pow(x,2) * ( 11 * pow(x,2) * pow(y,2) - pow(y,6) - 121 * pow(y,4) - 2 ))
+ ( b * pow(y,8))
+ ( x / ( 2 * y ));
}


yields -1.18059e+21. This answer has also been noted within MSExcel and Visual Fortran, even though CLisp on the same machine yields the correct answer. My question is, if we make an affine transformation of the numbers, will this lead to less round-off error, and are there alternatives to making this come closer to the correct answer in C++? --random [Edited by - random_thinker on August 21, 2005 6:31:20 AM]

##### Share on other sites
You can take the same approach as you used in Lisp - you could create a Rational class. Or you could pick up a multiprecision library such as GMP.

Also, keep in mind that the result of double operations are 80 bits, but when the result is written back, it is truncated to 64 bits. Which means that a given equation may not produce the same result depending on whether the compiler writes back intermediate results or does everything on the FP stack.

##### Share on other sites
Wouldn't a rational class require a big integer class as well? OItherwise wouldn't there be overruns?

--random

##### Share on other sites
If you are going to take the difference of large numbers then that's what's going to happen. It's well known in science that subtracting two large quantities produces large uncertainty in the result. Since the large terms in this equation are of order 10e36 that means the fractional error in your calculation is of order 10e-15. The C constant DBL_EPSILON is 2.2e-16 which means the error in the outcome is not significantly greater than the accuracy of doubles. In answer to you're question I've no idea if an affine transformation would help, but if it gets you anywhere near the right result it's probably more blind luck than anything else. You might do better by reordering terms to group identical powers together (that is, largest to smallest).

##### Share on other sites
One reason I suggested this is that I have seen a number of advantages in the field of computational molecular dynamics regarding the affine transformation of space co-ordinate systems. For example, the square root function can be reasonably estimated using a 3rd order equation in affine space, which could speed particle position computation on some platforms.

My worry about a rational class and the associated big integer types is speed; it is just quite slow. That is why I was wondering if affine transformation would help.

I've attempted to do this without success so far on the above equation but maybe I'm doing something wrong.

--random

Just adding; in molecular dynamics the affine transformation is one that transforms the co-ordinate system to a unit (0 to 1) basis. This permits many benfits when carrying out positional computation, including square root estimation.

[Edited by - random_thinker on August 20, 2005 11:00:51 AM]

##### Share on other sites
It has occurred to me that perhaps a better pow() function that allows fractional exponnents and perhaps encapsulates rational/bigint types might be a solution to error accumulation. Any ideas?

--random

• ### What is your GameDev Story?

In 2019 we are celebrating 20 years of GameDev.net! Share your GameDev Story with us.

• 15
• 9
• 11
• 9
• 9