# Twiddle table computing

I'm trying to optimize this twiddle computing function:

void twiddle(int N)
{
int i;
for (i=0;i<N;i++)
{
twiddle_table[i].re = (float)   cos((float)i * 2.0 * PI /(float)N);
twiddle_table[i].im = (float) - sin((float)i * 2.0 * PI /(float)N);
}
}


where N = 4096 size of twiddle table and could be bigger!

And then, I did the following:

void twiddle(int N)
{
int i;
float Tconst;
Tconst = 2.0 * PI /(float)N;
for (i=0;i<N;i++)
{
twiddle_table[i].re = (float)   cos((float)i * Tconst);
twiddle_table[i].im = (float) - sin((float)i * Tconst);
}
}


But, I'm getting a performance of 340,000 cycles for the for loop, which I think is bad.

Any hints that could enhance the performance of this function?

• Your question is quite unclear. What does it mean to "twiddle" something? Where is your twiddle_table defined? You didn't show all of your code. Also, what does 340 000 cycles for the for loop mean, when we don't even know the value of N or the size of the twiddle_table? Nov 22, 2016 at 23:03
• While the OP certainly should have explained himself, Google quickly told me that "twiddle table" is an established jargon term related to the implementation of FFTs (fast Fourier transforms). See e.g. en.wikipedia.org/wiki/Twiddle_factor Nov 22, 2016 at 23:56

The optimization you made manually is most likely already noticed and implemented by the compiler.

Instead, you should notice that $\cos\frac{2\pi k}{N} - i \sin\frac{2\pi k}{N} = e^{-\frac{2\pi i}{N}k}$, and replace multiple invocations of sin and cos with multiplications. Naively:

    complex base = { cos(2.0 * PI / N), -sin(2.0 * PI / N) };
twiddle_table[0] = { 1.0, 0.0 };
for (i = 1; i < N; ++i) {
twiddle_table[i] = twiddle_table[i - 1] * base;
}


In the production code you need to pay attention to the accumulation of numerical errors, and once in a while fall back to direct computation of sin and cos (or cexp).

Of course, if your compiler doesn't support complex type, you's have to implement complex multiplication manually.

Another optimization comes from the inherent symmetry of the table, e.g. $cos\frac{2\pi k}{N} = \cos\frac{2\pi (N - k)}{N}$, so you only need to compute half of the table; there are more symmetries you may exploit.

PS: That said, I strongly recommend to precompute it once for largest possible N. Notice that it can be reused for any power of 2.

• "In the production code you need to pay attention to the accumulation of numerical errors," - this is VERY important as multiplication quickly grows the error! Nov 23, 2016 at 8:00

A standard trick to reduce the error accumulation in vnp's code is to notice that the real part of base will be close to 1; we can use the multiple angle formula to compute base-1 (with base as in vnp's post) more accurately as

base_1 = -2*t*t + I * sin(a)


where

a = 2.0*PI/N
t = sin(a/2.0)


(I is the square root of -1). Then update your table with

tab[i] = tab[i-1] + tab[i-1]*base_1.


Though this more computation, it will have less rounding error.

Unless hardware rules it out you should consider computing base as doubles, and accumulating the complex exponentials in doubles; you can use one variable W say with

W += W*base_1


then stored W in tab[i] -- hence converting it to float. Again more computation but less accumulation of error.

this line:

Tconst = 2.0 * PI /(float)N;


will be performed in double due to the 2.0 to perform it in float, change the 2.0 to 2.0f

extract the expression:

(float)i*Tconst


to assign a float variable, inside the top of the for() loop and use that variable in the actual calculations

please show the definition of: twiddle_table[]

the parameter 'N' must NEVER be >= the number of entries in the twiddle)table[]. Otherwise, undefined behavior will occur as the data is being saved beyond the end of the array, which can lead to a seg fault event

Since the data and results are float, the function: cos() should not be used. Suggest: cosf().

the function: sin() should not be used. Suggest: sinf().

Note: then there would then be no need to cast the result of the call to cosf() nor sinf() to float