Using CTFE in D to speed up Sine and Cosine
By Danny Arends in D programmingPosted at: 01 Sep 2012, 21:35, last edited: 04 Sep 2012, 21:35
In 3D applications Sine and Cosine functions are used a lot to compute angles and rotation matrices. Usually these numbers are taken from a math library like std.math or something similar. But this involves many such function calls at each frame to compute the correct value.
Here however I show how to use the Compile Time Function Execution (CTFE) feature of the D 2.0 programming language to easily create Sine and Cosine look-up tables, with values pre-computed at compile time, making them constants. I'm going to start by defining some helper functions. This is because CTFE doesn't allow for unsafe operations at compile time. (see more about safe functions: http://dlang.org/safed.html)
To calculate Sine I use the Cordic algorithm (NOTE: it is also possible to use the std.math library for this (see on reddit)), and just translated the matlab code found on wikipedia: http://en.wikipedia.org/wiki/CORDIC to D 2.0 code. First I define a struct which will hold the sin and cosine of an angle in an array (T[2] d) and provide functions to (1) Negate the array and (2) multiply the array with a constant or (3) multiple by a 2x2 matrix. The struct has a Type parameter (T) making it so that it can hold either a float[] or a double[].
struct Coord(T = float) if (isFloatingPoint!T){
T[2] d = [1.0, 0.0];
T[2] opUnary(string s)() if (s == '-') {
d = [ -d[0], -d[1] ];
return d;
}
pure T[2] mult(in T a, in T precision = 1e-12){
d = [d[0] * a, d[1] * a];
if(d[0] < precision) d[0] = 0;
if(d[1] < precision) d[1] = 0;
return d;
}
pure T[2] mult(in T[2][2] R){
d = [ R[0][0]*d[0] + R[0][1]*d[1], R[1][0]*d[0] + R[1][1]*d[1] ];
return d;
}
}
We also need a two simple helper function to go from Degrees to Radians, and a function I call degreeloop. The degreeloop function shifts angles smaller then 0 and larger (or equal) to 360 inside the range of the array we are going to generate [0 .. 360].
pure U degToRad(U = float, V = int)(in V deg){
return to!U((deg * PI) / 180.0);
}
pure int degreeloop(int deg){
while(deg < 0 || deg >= 360){
(deg < 0) ? deg += 360 : deg -= 360;
}
return deg;
}
The Cordic algortihm in the D2.0 code is using the Coord struct, it has 2 parameters: 'beta' - the angle in radians and 'iter' - the number of iterations the algorithm will do. More iterations will lead to more precise results. It uses the same two global arrays as on the wikipedia page:
immutable Kangles = [ 0.78539816339745, ... , 0.00000000745058 ];
immutable Kvalues = [ 0.70710678118655, ... , 0.60725293500888 ];
I use a so called 'in{ } contract' to assert that the output variable type provided to the function can only be of a FloatingPoint type.
pure T[] cordic(T = float)(T beta, in uint iter = 50)
in{
assert(isFloatingPoint!T, 'T must be a FloatingPoint type');
}
body{
Coord!T v;
if(beta < -(.5*PI) || beta > (.5*PI)){
T newbeta = (beta < 0)? (beta + PI) : (beta - PI);
v.d = cordic!T(newbeta, iter);
return -v;
}
int sigma;
float Kn = Kvalues[min(iter, (Kvalues.length-1))];
float factor;
T R[2][2];
float onedivpowtwo = 1;
float angle = Kangles[0];
for(size_t j = 0; j < iter; j++){
sigma = (beta < 0)? -1 : 1;
factor = sigma * onedivpowtwo;
R = [[1, -factor],[factor, 1]];
v.mult(R);
beta = beta - sigma * angle;
onedivpowtwo = onedivpowtwo / 2.0;
if(j+1 > (Kangles.length-1)){
angle = angle / 2.0;
}else{
angle = Kangles[j+1];
}
}
return v.mult(Kn);
}
Now we have everything in place to generate our Sine and Cosine, now when we call:
cordic(degToRad(45));
We end up with a Coord struct back holding the Sine(45) in d[0] and Cosine(45) in d[1] in floats. To generate all Sines and Cosines from 0 to 360 we use the CTFE feature and put them into an enum or immutable. This is done by creating a simple loop function which is able to fill our immutable look-up table:
pure T[2][] gen_trigonometric(T = float)(in uint iter = 50)
in{
assert(isFloatingPoint!T, 'T must be a FloatingPoint type');
}
body{
T[2][] result = new T[2][](360);
foreach(i; 0 .. 360){
result[i] = cordic!T(degToRad!T(i), iter);
}
return result;
}
Then can we fill the table (with doubles) by calling the gen_trigonometric function.
immutable trigono = gen_trigonometric!double();
And finally you can generate 2 helper functions to get the Sine and Cosine of any angle by using the degreeloop function, and the look-up table, these functions return auto because based on the gen_trigonometric it could return a double or a float:
pure auto sin(int deg){
deg = degreeloop(deg);
return trigono[deg][0]; // 0 = SIN
}
pure auto cos(int deg){
deg = degreeloop(deg);
return trigono[deg][1]; // 1 = COSINE
}
Using these two function will now avoid any calculation (Besides the shifting into the 0 .. 360 range) and use the CTFE computed values to improve performance, You can easily extend this approach to generate all possible 360 Pitch, Roll and Yaw matrices used in Euler geometry.
Find the whole code as gist on github: https://gist.github.com/3526430
Last modified: 04 Sep 2012, 21:35 | Edit