I need to calculate the sin and cos of a lot of large (1E+6 or so) numbers. Do I have any other options besides just using __sincosf() and sincosf()? e.g. is it more efficient to subtract an integer multiple of 2pi first?
may be of interest:
The need to calculate trigonometric functions of large arguments is definitely a red flag, numerically speaking.
Are these single-precision operands? If so, when applying accurate trigonometric functions to operands of that magnitude, your results will have an effective precision of about 4 bits (about twenty leading mantissa bits will cancel during argument reduction), meaning you are processing mostly noise. If you switch to __sincosf(), which doesn’t even use accurate argument reduction, the output will be nothing but noise.
Why are these inputs so large? Could you use sincospi() instead? Note that this helps processing large arguments faster and accuracy could be slightly better, but the issue of subtractive cancellation during argument reduction will not go away, meaning effective precision in the results will be low.
Note that the basic problem with large arguments is primarily the limited granularity of the floating-point representation. If my mental arithmetic is correct, arguments of magnitude 1e6 represented in IEEE 754 single precision resolve angles to about the closest multiple of 25 degrees.
Yes, they are single-precision operands. The operand represents the phase of an electromagnetic wave. I have multiple waves with different phase shifts and need to predict whether constructive of destructive interference will occur.
If you are right about the precision of my phase only being about 25 degrees then I definitely need to rethink how I’m doing things.
Thanks for the tips.
You can easily find out what the minimal separation of arguments of large magnitude is, e.g.
printf ("smallest representable difference near %15.8e is %15.8e\n", x, x - nextafterf(x, 0.0f));
That gives you the granularity in radians. [Later:] Seems like I was off by an order of magnitude in my mental arithmetic. The above code tells me that near x = 1e6, the resolution of the IEEE-754 single-precision floating-point format is 6.25e-2 = 1/16 (equivalent to 3.6 degrees).
Note that manually removing multiples of 2π runs the risk of phase errors that grows as time increases, since 2π can’t be represented accurately in a floating-point number. In addition, you incur rounding error every time you remove a multiple of 2π, and those rounding error accumulate.
Another thing to consider, in addition to floating-point quantization, is inherent error in whatever data your are processing. The frequency accuracy for radiation sources may be 1 ppm at best, and your sensors may only have 16 bit resolution (or, more commonly, 12 bits). In which case some low-order bits of such data stored in single-precision operands would be noise, and one would not want to compute with noise.
I wonder whether mathematical treatment of this issue is possible. Given two sine waves of different phase, it should be possible to determine function that indicates the effective result of interference. Something akin to the beat frequency that results from the interference of two sine waves with different frequency. Now, if we are unlucky, such a function (if it exists) could involve subtractive cancellation, so that numerically, we are looking at the same problem again, just in a different disguise. In this case, we might be interested in the equation of the envelope function?