-1

Here is a part in a Physics engine.
The simplified function centerOfMass calculates 1D-center-of-mass of two rigid bodies (demo) :-

#include <iostream>
#include <iomanip> 
float centerOfMass(float pos1,float m1, float pos2,float m2){
    return (pos1*m1+pos2*m2)/(m1+m2);
}
int main(){
    float a=5.55709743f;
    float b= centerOfMass(a,50,0,0);
    std::cout << std::setprecision(9) << a << '\n';  //5.55709743
    std::cout << std::setprecision(9) << b << '\n';  //5.55709696
}

I need b to be precisely = 5.55709743.

The tiny difference can, sometimes (my real case = 5%), introduces a nasty Physics divergence.
There are some ways to solve it e.g. heavily do some conditional checking.
However, it is very error-prone for me.

Question: How to solve the calculation error while keep the code clean, fast, and still easily to be maintained?

By the way, if it can't be done elegantly, I would probably need to improve the caller to be more resistant against such numerical error.

Edit

(clarify duplicate question)

Yes, the cause is the precision error from the storage/computing format (mentioned in Is floating point math broken?).

However, this question asks about how to neutralize its symptom in a very specific case.

cppBeginner
  • 1,114
  • 9
  • 27
  • 2
    Possible duplicate of [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) – Code-Apprentice Jan 21 '18 at 04:14
  • Note that floating point numbers are an inexact approximation to real numbers. Operations with floating point numbers have many well-known limitations. – Code-Apprentice Jan 21 '18 at 04:16
  • @Code-Apprentice Hmm, it is similar. In my case, I know the cause and I want a workaround, but that question focus on the numerical cause. – cppBeginner Jan 21 '18 at 04:16
  • *How to solve the calculation error* what error? – CroCo Jan 21 '18 at 04:20
  • @CroCo : `5.55709743 != 5.55709696` , i.e. a should == b. – cppBeginner Jan 21 '18 at 04:21
  • 2
    you should not rely on `==` to compare floating-point numbers. In your case, switching from float to double solves your particular case as stated by @Aziz, but you should seriously address the limitations of floating-point representation, otherwise, you will be surprised with the results. – CroCo Jan 21 '18 at 04:39
  • @CroCo Agree, Aziz's can solve the specific case. .... In my case, even I never compare it directly, the tiny error can indirectly cause a very small angular velocity on a rigid body that has zero inertia. It is forbidden in my program. (it is so complex) .... I just dream for a beautiful way that I don't have to check. Now, I come to believe such thing doesn't exist. ... Thank. – cppBeginner Jan 21 '18 at 04:44
  • 1
    @cppBeginner -- The real issue is that you took a "math book" calculation, and without modification, wrote the computer program equivalent, expecting the same results. Floating point isn't a math book -- in many cases doing a direct translation will not work for a binary computing machine due to round-off error. Many formulas had to be rewritten to account for this issue, so in a nutshell, you should do the research so that you can rearrange the calculations in some way to reduce the error. – PaulMcKenzie Jan 21 '18 at 06:52
  • @PaulMcKenzie I believe what you meant : there is a solution, and it should not be very hard to find? As far as I know, there are some techniques e.g. change order of operation, apply distributivity, do some checking if-else. None of them seems to achieve my objective : beautiful, precise and simple. I googled "weight average precision/numerical error" , found nothing useful. May you provide some clues, please? I probably miss something that is very obvious. – cppBeginner Jan 21 '18 at 07:18
  • If you think that you really need accuracy in the seventh digit for your "physics engine", then there's something wrong with that "engine". For example, a bad ODE integration scheme will give wrong results that depend very sensitively on minute variations of the initial conditions. In classical physics such small differences only have a measurable impact in the case of chaotic systems. – RHertel Jan 21 '18 at 08:29
  • @RHertel Good to hear. I fixed other parts instead, and it works now. :) – cppBeginner Jan 21 '18 at 10:43

2 Answers2

2

You are trying to get 9 decimal digits of precision , but the datatype float has a precision of about 7 decimal digits.

Use double instead. (demo)

Aziz
  • 20,065
  • 8
  • 63
  • 69
1

Use double, not float. IEEE 754 double has about 16 decimal places of precision.

#include <iostream>
#include <iomanip> 
double centerOfMass(double pos1, double m1, double pos2, double m2) {
    return (pos1*m1 + pos2 * m2) / (m1 + m2);
}
int main() {
    double a = 5.55709743;
    double b = centerOfMass(a, 50, 0, 0);
    std::cout << std::setprecision(16) << a << '\n';  //5.55709743
    std::cout << std::setprecision(16) << b << '\n';  //5.55709743
    std::cout << std::setprecision(16) << (b - a) << '\n';  // 0
}

For the example given, centerOfMass(a, 50, 0, 0), the following will give exact results for all values of a, but of course the example does not look realistic.

double centerOfMass(double pos1, double m1, double pos2, double m2) {
    double divisor = m1 + m2;
    return pos1*(m1/divisor) + pos2*(m2/ divisor);
}
Jive Dadson
  • 16,680
  • 9
  • 52
  • 65
  • Will it still have the same issue - but in a different numerical degree? – cppBeginner Jan 21 '18 at 08:10
  • @cppBeginner - Yes, always. If your a's are "close" to 1.0, the difference between a and b will be of the order 1e-15 (assuming IEEE 754 double). Taking an example off the top of my head, a=6.7820743978696196 will yield b=6.7820743978696187 – Jive Dadson Jan 21 '18 at 08:13
  • Thank for the number! Can I get rid ALL of the issue? I mean - not just alleviate it. PaulMcKenzie (comment near OP's) somehow mentioned that it is possible to re-arrange the formula, but I don't know how. – cppBeginner Jan 21 '18 at 08:14
  • @cppBeginner Can't be done without an infinite amount of memory, and who can afford that? – Jive Dadson Jan 21 '18 at 08:15
  • @cppBeginner - Rearranging the formula is not likely to help much if any. – Jive Dadson Jan 21 '18 at 08:20