Re: complexity of trigonometric functions
On 6 Sep., 14:33, Vladimir Jovic wrote:
SG wrote:
In case you want to compute a sequence of sine values a la
f[k] = sin(a+b*k)
you can make use of a simple linear recurrence:
f[0] = sin(a)
f[1] = sin(a+b)
f[k] = 2 cos(b) f[k-1] - f[k-2]
where 2*cos(b) is obviously a constant. And since cos(x)=sin(pi/2+x)
you can use this to compute cosine values as well, of course. Only two
FLOPs per value is a very low number of operations. The only drawback
is rounding errors. Every one in a while one should start new or
somehow "normalize" the values so that the amplitude doesn't decrease
or increase due to rounding errors.
That is exactly what I am doing. Thanks for the tip.
btw is 1000 iterations (k=1000) going to produce big error?
I'm not sure. The rounding error in the precomputed constant for
2*cos(b) obviously affects the frequency to some degree. And a
rounding error 'e' in one of the samples produces an infinite response
(a sinusoid with amplitude e/sin(b)). Assuming the rounding errors are
not correlated and ignoring the frequency error, the expected absolute
error for the k-th sample should be proportional to sqrt(k)/sin(b).
If the frequency dependency (1/sin(b)) is bugging you and you need
both, sin AND cos an alternative which is almost as fast would be to
use the following complex-valued 1st order linear recurrence:
c[k] = exp(i*(a+b*k))
= exp(i* b ) * c[k-1]
where -1=i*i, real(c[k])=cos(a+b*k) and imag(c[k])=sin(a+b*k). With
this recurrence you need 4 multiplications and two additions for a
pair of sin/cos values. That's twice as many multiplications compared
to above. You still have the sqrt(k)-scaling but the frequency
dependency is gone now (yay!).
Since you're concerned about speed, let me share my recent experiences
with GCC's std::complex<double> implementation: Without a compiler
option like -ffast-math operator* maps to a function called __muldc3
(according to the profiler) which is terribly slow compared to 4 MULs
and 2 ADDs for some reason. I also tested operator* for std::complex<>
using a recent Microsoft compiler. According to my measurements
operator* was only half as fast compared to 4 MULs and 2 ADDs on an
old Pentium 4 HT. I guess there's some special case handling involved
with respect to NaN or Inf that is slowing down the multiplication.
Cheers!
SG