On Jan 5, 10:16*pm, Rune Allnor <(EMail Removed)> wrote:
> On 6 Jan, 04:11, bei <(EMail Removed)> wrote:
>
>
>
> > On Jan 5, 7:09*pm, bei <(EMail Removed)> wrote:
>
> > > Hello,
>
> > > I was bothered by the following problem for a while and could not
> > > understand how it happened.
>
> > > They are a lot of floating point calculation in the code.
> > > Basically, i am doing an interpolation from a random position to the
> > > closest grid point position. It is a 1d problem.
> > > The domain is from (L, L) and N is the number of grid point at a half
> > > domain.
> > > The grid size is dv, where dv = L/N; The grid is cell centered (0.5
> > > in the code below).
> > > L, dv are declared to be double type and N to be integer type.
>
> > > Below is a part of the code that bothered me (where v is the position,
> > > j is the closest grid point index number)
>
> > > *int j0 = (int)(floor((v/dv+(N0.5))));
> > > *double a = 2.5*pow(((v/dv+(N0.5))j0), 2.0);
> > > *double b = 1.5*pow(((v/dv+(N0.5))j0), 3.0);
> > > *double q = 1.02.5*pow(((v/dv+(N0.5))j0), 2.0)+1.5*pow(((v/dv+
> > > (N0.5))j0), 3.0);
> > > *double q0 = 1.0 +a + b; (q0 is equal to q in the language of math.)
>
> > > the output is
> > > v = 2.818975726362143e+00;
> > > dv = 4.394531250000000e03;
> > > N = 1146;
> > > j0 = 1786;
> > > a = 2.369681598890769e+00;
> > > b = 1.384255444361780e+00;
>
> > > nothing wrong until here, below is the problem
> > > q = 1.02.5*pow(((v/dv+(N0.5))j0), 2.0)+1.5*pow(((v/dv+(N0.5))j0),
> > > 3.0)
> > > * *= 1.369681598890769e+00; (wrong answer)
> > > q0 = 1.0 +a + b
> > > * * = 1.457384547101137e02; (correct answer)
>
> > > It looked like if a formula involvs a lot of floating point
> > > calculation, it is safer to save the temporary data into a variable.
>
> > > Why this is happened and how could I fix this problem? Any suggestion
> > > is appreciated. Thank you!
>
> > btw, i am using g++ 4.3.3 with O3 optimization turned on. The machine
> > is 64bits.
>
> The first place to look is data types. You have j0 as an int,
> and all the other factors as doubles. That's the kind of thing
> that almost always causes surprises and trouble. Declare j0 as
> double first and try again. Repeat for all integers in the
> computations.
>
> Next have a look at the pow function. Since you use powers of
> small integers, do something like
>
> pbase = (((v/dv+(N0.5))j0);
> p2 = pbase*pbase;
> p3 = p2*pbase;
>
> instead of using the pow() function. This would almost
> certainly be significantly faster than pow().
>
> As for why the problem occurs, I would *guess* that the
> optimizer recognizes that pow() is called twice in the same
> statement, with the same mantissa and different exponents.
> If this is a template implementation (which might well be
> the case), the optimizer might try to reduce the number of
> repeated computations in the two calls, messing things up
> in the process.
>
> Rune
Rune, very helpful suggestion. Thank you.
