Geekly Articles each Day

Most likely, you know the following ratios from school:

$$ sin( alpha+ beta)= sin alpha times cos beta+ cos alpha times sin beta cos( alpha+ beta)= cos alpha times cos beta− sin alpha times sin beta

When you first became acquainted with this formula in childhood, most likely, your first feeling was pain due to the fact that this formula must be remembered. This is very bad, because in fact you do not need to remember this formula - it is displayed when you rotate the triangle on paper. In fact, I do the same when I write down this formula. This interpretation will be apparent by the middle of this article. But now, to leave all the fun for later and to postpone the moment when you say “Eureka!”, Let's think about why we should even think about this formula.

The trigonometric functions `sin()`

and `cos()`

probably the most popular in computer graphics, since they are the basis for describing any round shape in a parametric way. Among the places of their possible application are the generation of circles or volumetric objects of rotation, when calculating the Fourier transform, procedural generation of waves on the water plane, generators for a software sound synthesizer, and so on. In all these cases, `sin()`

and `cos()`

are called inside the loop, as here:

` ``for(int n=0; n < num; n++) { const float t = 2.0f*PI*(float)n/(float)num; const float s = sinf(t); const float c = cosf(t); // do something with s and c ... }`

We begin to rewrite the cycle incrementally (see the code below), so it is easier for us to imagine that at iteration `n`

this cycle with phase `t`

, the next iteration, `n+1`

, will consider `sin()`

and `cos()`

for `t+f`

. In other words, `sin(t)`

and `cos(t)`

are counted for us and we need to count `sin(t+f)`

and `cos(t+f)`

:

` ``const float f = 2.0f*PI/(float)num; const float t = 0.0f; for( int n=0; n < num; n++ ) { const float s = sinf(t); const float c = cosf(t); // do something with s and c ... t += f; }`

It doesn't matter how we got `t`

and what its range of values is (in the example above - $$[0;2 pi] ) The only thing that bothers us is that there is a loop that constantly calls `sin()`

and `cos()`

with a parameter that increases in constant steps (in this case, $$ frac2 pi textnum ) This article is about how to optimize this code for speed in such a way that the same calculations can be performed without using the `sin()`

or `cos()`

functions (in the inner loop), or even the faster `sincos()`

function.

But if you look at the first formula in the article, we will see that if $$f= alpha and $$t= beta we can rewrite it as

` ``sin(t+f) = sin(f)*cos(t) + cos(f)*sin(t) cos(t+f) = cos(f)*cos(t) - sin(f)*sin(t)`

or in other words

` ``new_s = sin(f)*old_c + cos(f)*old_s new_c = cos(f)*old_c - sin(f)*old_s`

Since `f`

is a constant, `sin(f)`

and `cos(f)`

also. Call them `a`

and `b`

respectively:

` ``new_s = b*old_c + a*old_s new_c = a*old_c - b*old_s`

This expression can be directly used in our code, and then we get a loop in which expensive (and indeed no) trigonometric functions are called!

` ``const float f = 2.0f*PI/(float)num; const float a = cosf(f); const float b = sinf(f); float s = 0.0f; float c = 1.0f; for( int n=0; n < num; n++ ) { // do something with s and c ... const float ns = b*c + a*s; const float nc = a*c - b*s; c = nc; s = ns; }`

To date, we have blindly played with math, not understanding what is really happening. Let's rewrite the inner loop like this:

$$sn+1=sna+cnbcn+1=cna−snb

Some of you may have noticed that this is a formula for rotating an object in two-dimensional space. If you still do not understand this, perhaps the matrix form will help you.

$$\ left (\ begin {matrix} s_ {n + 1} \\ c_ {n + 1} \ end {matrix} \ right) = \ left (\ begin {matrix} a & b \\ -b & a \ end {matrix} \ right) \ cdot \ left (\ begin {matrix} s_ {n} \\ c_ {n} \ end {matrix} \ right)

\ left (\ begin {matrix} s_ {n + 1} \\ c_ {n + 1} \ end {matrix} \ right) = \ left (\ begin {matrix} a & b \\ -b & a \ end {matrix} \ right) \ cdot \ left (\ begin {matrix} s_ {n} \\ c_ {n} \ end {matrix} \ right)

In fact, `sin(t)`

and `cos(t)`

can be grouped into a vector of length 1 and drawn as a point on the screen. Call this vector `x`

. Then, $$x = \ {\ sin \ beta, \ cos \ beta \}x = \ {\ sin \ beta, \ cos \ beta \} . So the vector form of expression is $$xn+1=Rxn where $$R = \ left (\ begin {matrix} a & b \\ - b & a \ end {matrix} \ right)R = \ left (\ begin {matrix} a & b \\ - b & a \ end {matrix} \ right) . We see that our loop performs a small two-dimensional rotation each iteration so that `x`

rotates in a circle during the execution of the cycle. Each iteration `x`

rotates by $$ frac2 pi textnum radian.

So basically

$$ sin( alpha+ beta)= sin alpha times cos beta+ cos alpha times sin beta cos( alpha+ beta)= cos alpha times cos beta− sin alpha times sin beta

this is the point motion formula $$x = \ {\ cos \ alpha, \ sin \ alpha \}x = \ {\ cos \ alpha, \ sin \ alpha \} circumference in increments of $$ beta radian. To do this, we will build one of two rotation axes, for example, $$u = \ {\ cos \ beta, \ sin \ beta \}u = \ {\ cos \ beta, \ sin \ beta \} . The first component of the rotation is the projection. $$x on $$u . Because $$x and $$u normalized (have a length of 1), the projection is their scalar product. Consequently, $$s=x cdotu= sin alpha cdot cos beta+ cos alpha cdot sin beta , and of course the second component is the anti-projection, which can be found by projecting onto the perpendicular axis, $$v . We can create this vector by expanding the coordinates $$u and change the sign to the opposite at the first coordinate: $$c=x cdotv= cos alpha cdot cos beta+ sin alpha cdot sin beta

Usually you should be able to perform these tiny rotations over and over again. Indeed, $$| R | = \ left | \ begin {matrix} a & b \\ - b & a \ end {matrix} \ right | = a ^ 2 + b ^ 2 = \ sin ^ 2 \ alpha + \ cos ^ 2 \ alpha = 1| R | = \ left | \ begin {matrix} a & b \\ - b & a \ end {matrix} \ right | = a ^ 2 + b ^ 2 = \ sin ^ 2 \ alpha + \ cos ^ 2 \ alpha = 1 , which means that the matrix $$R does not increase or decrease the space to which it is applied, which means that $$x will move in a perfect circle. However, because computers are not accurate, $$x will move in a spiral and eventually coincides with the center of the circle of rotation. I did not have such problems, but I think that they can occur with very large `num`

, i.e. small turning angles.

In Kindercrasher , the 4096-byte demo from 2008 (screenshot on KDPV), a group of spheres pulsates to the music. For this, I counted the Fourier transform of sound. There are algorithms that do this in real time, for example, FFT . However, my code should fit in a few kilobytes, and I decided to go the other way. Instead of implementing FFT, I wrote DFT by its simple definition. You can check it out on wikipedia.

$$Xk= sumN−1n=0xne− frac2 piiNkn quadk=0,1, ldots,N−1

My function also takes a 16-bit stereo audio buffer, `x`

, and returns the first 128 frequencies of the sound spectrum of sound `y`

. See how the inner loop is organized, the one that runs 4096 times: not a single call to the `sin()`

or `cos()`

functions, although in other implementations these calls will be.

` ``#include <math.h> void iqDFT12( float *y, const short *x ) { for( int i=0; i<128; i++ ) { const float wi = (float)i*(2.0f*3.1415927f/4096.0f); const float sii = sinf( wi ); const float coi = cosf( wi ); float co = 1.0f; float si = 0.0f; float acco = 0.0f; float acsi = 0.0f; for( int j=0; j<4096; j++ ) { const float f = (float)(x[2*j+0]+x[2*j+1]); const float oco = co; acco += co*f; co = co*coi - si*sii; acsi += si*f; si = si*coi + oco*sii; } y[i] = sqrtf(acco*acco+acsi*acsi)*(1.0f/32767.0f); } }`

Source: https://habr.com/ru/post/462725/