Jump to content
  • Advertisement
Sign in to follow this  
random_thinker

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

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

If you intended to correct an error in the post then please contact us.

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 this post


Link to post
Share on other sites
Advertisement
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 this post


Link to post
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 this post


Link to post
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 this post


Link to post
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

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!