Reimar Döffinger Reimar.Doeffinger
Thu May 14 15:43:32 CEST 2009

```On Thu, May 14, 2009 at 11:00:48AM +0200, Reimar D?ffinger wrote:
> On Thu, May 14, 2009 at 10:29:37AM +0200, 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.
>
> Note that I think you can get a) with
>
> alpha = 0.125 / n;
> if (scale < 0)
>     alpha += 0.25;
> for ....
> cos(2*M_PI*alpha) * scale;
> sin(2*M_PI*alpha) * scale;
> alpha += 1.0 / n;

It could be considered at bit obfuscated, but even better might be:
// sampling_pos from 0 to 8n corresponds to a full period of sin/cos
int sampling_pos = 1;
if (scale < 0)
sampling_pos += 8*n/4; // shift by 90 degrees
scale = -sqrt(FFABS(scale));
for ...
// convert from 0 .. 8n to 0 .. 2pi
alpha = sampling_pos * (2 * M_PI / (8 * n));
cos(alpha) * scale;
sin(alpha) * scale;
sampling_pos += 8;

```