CODEWITHSUNDEEP
c++precision
G. Sophys
G. Sophys

Reputation: 3

Precision concerning double variable in C++

My problem is when I calculate one variable like this:

UXZ=(DRDe*sY*(Ux2 + Uy2) + cY*DRDt*Tau*(Ux2 + Uy2) + DRDx*Tau*Ux*Sqrt(1 + Ux2 + Uy2) + DRDy*Tau*Uy*Sqrt(1 + Ux2 + Uy2))/(Tau*Sqrt(1 - 1/(1 + Ux2 + Uy2)));

This value 1/(1 + Ux2 + Uy2) is set equal to (strictly) 1. Then, we have an error because the denominator is equal to 0.

However, the terms Ux2 and Uy2 are equal to (each) 1E-32. Consequently, the value is not strictly equal to 1.

I don't know how to set the precision on the calculation. Every variable in this equation are double, maybe I have to put them in long?

Edit : This calculation is implemented in a event generator, and Ux2 and Uy2 will change/increase during different times/iterations. We start with a very small value of Ux2 and Uy2, very close to 0 but not strictly 0. If we keep (1 + 1E-32 + 1E-32)=1. that fully break my code.

The important code for my question is this :

UXZ=XX/(YY*Sqrt(1 - 1/(1 + Ux2 + Uy2)));

I cannot change the value of Ux2 and Uy2 because that will change the physics problem under consideration.

Upvotes: 0

Views: 80

Answers (1)

TripeHound
TripeHound

Reputation: 2960

Your "problem" term is:

SQRT( 1 - 1 / ( 1 + Ux2 + Uy2 ) )

Since SQRT(a) is identical to SQRT(b*b * a) / b (i.e. multiply inside the square-root by b-squared, and divide the result by b), then – using (1 + Ux2 + Uy2) for b – you get:

SQRT( (1 + Ux2 + Uy2)*(1 + Ux2 + Uy2) - (1 + Ux2 + Uy2) ) / (1 + Ux2 + Uy2)

which should evaluate to the same value, without triggering a divide-by-zero.

It would probably also be useful to bring out some of the common expressions (e.g. 1 + Ux2 + Uy2) into variables.

To understand why you need to do this, you need to be looking into the area of Numerical Analysis (Wikipedia), in particular, Generation and propagation of errors, which tackle the problem of representing arbitrary real numbers with the limited-precision schemes used in most computer languages.

I'm far from an expert in this, but at it's core, you have to recognise that while 1 - 1 / ( 1 + x ) will always be "valid", assuming x is at least a little-bit positive, in most computer implementations, 1 + x will degenerate to exactly 1 for a small-enough x. Part of Numerical Analysis is recognising this, and reorganising how you perform a calculation to make better use of limited precision.

Upvotes: 1

Related Questions