Uoti Urpala uoti.urpala
Thu May 14 16:33:43 CEST 2009

```On Thu, 2009-05-14 at 15:36 +0200, Reimar D?ffinger wrote:
> On Thu, May 14, 2009 at 12:09:31PM +0300, Siarhei Siamashka wrote:
> > On Thu, May 14, 2009 at 11:29 AM, Reimar D?ffinger wrote:
> > > On Thu, May 14, 2009 at 11:06:21AM +0300, Siarhei Siamashka wrote:
> > >>
> > >> +    alpha = (scale < 0 ? 0.5 * M_PI : 0) + 0.25 * M_PI / n;
> > >> +    scale = -sqrt(fabs(scale));
> > >>      for(i=0;i<n4;i++) {
> > >> -        alpha = 2 * M_PI * (i + 1.0 / 8.0) / n;
> > >> -        s->tcos[i] = -cos(alpha);
> > >> -        s->tsin[i] = -sin(alpha);
> > >> +        s->tcos[i] = cos(alpha) * scale;
> > >> +        s->tsin[i] = sin(alpha) * scale;
> > >> +        alpha += 2.0 * M_PI / n;
> > >
> > > Rounding errors accumulate like that, IMO you really should never do
> > > that unless a) you know the value you add can be represented exactly
> > > or b) you absolutely can't avoid it for performance reasons.
> >
> > Everything is calculated in double precision floating point, with the
> > final results stored in single precision.
>
> Have you actually checked how bad it is?
> Assuming only the least significant bit of "2.0 * M_PI / n" is wrong,
> and given we support up to nbits == 18, that means the 16 least
> significant bits of alpha will be wrong at the end of the loop, even
> if you do not include possible rounding errors in the addition.

This analysis is completely bogus. X*i/Y is no better than starting with
s = 0 and adding X/Y 'i' times _if_ the addition is done with infinite
precision. Yes, the absolute error in X/Y is multiplied by i - but the
leading digit is multiplied by the same amount, meaning the error is
still in the least significant digit of the result given _its_ magnitude
as opposed to the magnitude of X/Y.