So in the implementation of cos_table_*_LERP, you did technically 2 step range reduction:
1. Reduce x = x mod 2*pi
2. Reduce index = floor(x / 10^-n), and i - index = 10^n * (x mod 10^-n)
With limited input range and required precision as in the tests, you can combine these 2 range reduction steps:
1. Choose the reduced range as power of 2 instead of power of 10 for cheaper modulus operation, let say 2^-N = 2^-7.
2. Avoid the division in modd(x, CONST_2PI) by multiplying by 2^N / pi.
3. Avoid the round trip double -> int -> double by using the floor function / instruction.
Here is the updated version of cos_table_*_LERP which should have higher throughput and lower latency:
double cos_table_128_LERP(double x) {
x = fabs(x);
double prod = x * TWO_TO_SEVEN_OVER_PI;
double id = floor(prod);
double x = prod - id; /* after this step, 0 <= x < 1 */
int i = ((int)id) & 0xFF; /* i = id mod 2^8 */
return lerp(x, COS_TABLE[i], COS_TABLE[i + 1]);
}
You can also optimize lerp a bit more with the formula:
We do employ this range reduction strategy in a more accurate way for trig functions in LLVM libc:
- with FMA instructions: https://github.com/llvm/llvm-project/blob/main/libc/src/math/generic/range_reduction_fma.h
- without FMA instructions: https://github.com/llvm/llvm-project/blob/main/libc/src/math/generic/range_reduction.h
Your intuition is actually proved mathematically for two-player games of perfect information (even for infinite games), in which chess is included, by Wolfgang Schmidt (Theorem 7 in https://www.jstor.org/stable/1994619). Roughly speaking, if a two-player games of perfect information has a winning strategy (which might need all previous moves of a play to make decisions), then it has a positional winning strategy, which only need the current state/position to make decisions. The `only` requirement needed for that theorem is the Axiom of Choice / Well-ordering Principle.
Maybe correct rounding (and hence consistency) will be made into the standard requirements (at least for single and double precisions) in the near future?
I don't know much about numpy.signal package, but there are other python packages that have good implementation for generating polynomial approximation. Among them are pythonsollya (a Python wrapper for Sollya), and SageMath.
Other than that, some of the overview references about Remez algorithm (should be a bit closer to your taste) that I like are:
- J.M. Muller et. al., "Handbook of Floating-Point Arithmetic", section 10.3.2 (and follow references in that section).
I think the function inputs are in double-double format, so the assumption is that the magnitude of y is significantly smaller than the magnitude of x, ideally |y| < 2^(-52) |x|. So that 1 is a very good approximation for cos(y), since |cos(y) - 1| < |y|^2 < 2^(-104) |x|^2. Similarly, if we assume |x| is also small, |sin(x) - x| < |x|^3 is also a very good approximation (same with sin(y) ~ y).
Some of the more recent implementations of cosine for single precision that are correctly rounded to all rounding modes according to IEEE 754 standard:
1. Reduce x = x mod 2*pi
2. Reduce index = floor(x / 10^-n), and i - index = 10^n * (x mod 10^-n)
With limited input range and required precision as in the tests, you can combine these 2 range reduction steps:
1. Choose the reduced range as power of 2 instead of power of 10 for cheaper modulus operation, let say 2^-N = 2^-7.
2. Avoid the division in modd(x, CONST_2PI) by multiplying by 2^N / pi.
3. Avoid the round trip double -> int -> double by using the floor function / instruction.
Here is the updated version of cos_table_*_LERP which should have higher throughput and lower latency:
You can also optimize lerp a bit more with the formula:
We do employ this range reduction strategy in a more accurate way for trig functions in LLVM libc: