7
\$\begingroup\$

I have been developing a real time application in C and I need to calculate sine and cosine of a given angle in radians. In respect to the fact that the application is a real time application the functions from the standard library are not suitable for my purposes due to the fact that there is not fixed calculation time. Based on that I have decided to implement those functions in my own.

From the execution time point of view it seems to me that the best approach for real time application is to use the look-up table based method of calculation with improvement of precion based on linear interpolation. I have exploited the symmetry of the sine and cosine functions and I have defined the look-up table covering the quarter of the period.

Look-up table

#define TABLE_SIZE 256
#define PI         3.14
#define STEP_SIZE  (PI/(2*(TABLE_SIZE-1))) // N points divide PI/2 into N-1 segments

// look-up table containing values of sine function over one quarter of period
// angel step size is pi/(2*256) i.e. pi/512 approximately 0.35°
double sinLUT[TABLE_SIZE] = 
{
0.000000,  0.006160,  0.012320,  0.018479,  0.024637,  0.030795,  0.036951,  0.043107,  
0.049260,  0.055411,  0.061561,  0.067708,  0.073853,  0.079994,  0.086133,  0.092268,  
0.098400,  0.104528,  0.110653,  0.116773,  0.122888,  0.128999,  0.135105,  0.141206,  
0.147302,  0.153392,  0.159476,  0.165554,  0.171626,  0.177691,  0.183750,  0.189801,  
0.195845,  0.201882,  0.207912,  0.213933,  0.219946,  0.225951,  0.231948,  0.237935,  
0.243914,  0.249883,  0.255843,  0.261793,  0.267733,  0.273663,  0.279583,  0.285492,  
0.291390,  0.297277,  0.303153,  0.309017,  0.314870,  0.320710,  0.326539,  0.332355,  
0.338158,  0.343949,  0.349727,  0.355491,  0.361242,  0.366979,  0.372702,  0.378411,  
0.384106,  0.389786,  0.395451,  0.401102,  0.406737,  0.412356,  0.417960,  0.423549,  
0.429121,  0.434676,  0.440216,  0.445738,  0.451244,  0.456733,  0.462204,  0.467658,  
0.473094,  0.478512,  0.483911,  0.489293,  0.494656,  0.500000,  0.505325,  0.510631,  
0.515918,  0.521185,  0.526432,  0.531659,  0.536867,  0.542053,  0.547220,  0.552365,  
0.557489,  0.562593,  0.567675,  0.572735,  0.577774,  0.582791,  0.587785,  0.592758,  
0.597707,  0.602635,  0.607539,  0.612420,  0.617278,  0.622113,  0.626924,  0.631711,  
0.636474,  0.641213,  0.645928,  0.650618,  0.655284,  0.659925,  0.664540,  0.669131,  
0.673696,  0.678235,  0.682749,  0.687237,  0.691698,  0.696134,  0.700543,  0.704926,  
0.709281,  0.713610,  0.717912,  0.722186,  0.726434,  0.730653,  0.734845,  0.739009,  
0.743145,  0.747253,  0.751332,  0.755383,  0.759405,  0.763398,  0.767363,  0.771298,  
0.775204,  0.779081,  0.782928,  0.786745,  0.790532,  0.794290,  0.798017,  0.801714,  
0.805381,  0.809017,  0.812622,  0.816197,  0.819740,  0.823253,  0.826734,  0.830184,  
0.833602,  0.836989,  0.840344,  0.843667,  0.846958,  0.850217,  0.853444,  0.856638,  
0.859800,  0.862929,  0.866025,  0.869089,  0.872120,  0.875117,  0.878081,  0.881012,  
0.883910,  0.886774,  0.889604,  0.892401,  0.895163,  0.897892,  0.900587,  0.903247,  
0.905873,  0.908465,  0.911023,  0.913545,  0.916034,  0.918487,  0.920906,  0.923289,  
0.925638,  0.927951,  0.930229,  0.932472,  0.934680,  0.936852,  0.938988,  0.941089,  
0.943154,  0.945184,  0.947177,  0.949135,  0.951057,  0.952942,  0.954791,  0.956604,  
0.958381,  0.960122,  0.961826,  0.963493,  0.965124,  0.966718,  0.968276,  0.969797,  
0.971281,  0.972728,  0.974139,  0.975512,  0.976848,  0.978148,  0.979410,  0.980635,  
0.981823,  0.982973,  0.984086,  0.985162,  0.986201,  0.987202,  0.988165,  0.989092,  
0.989980,  0.990831,  0.991645,  0.992421,  0.993159,  0.993859,  0.994522,  0.995147,  
0.995734,  0.996284,  0.996795,  0.997269,  0.997705,  0.998103,  0.998464,  0.998786,  
0.999070,  0.999317,  0.999526,  0.999696,  0.999829,  0.999924,  0.999981,  1.000000
};

Sine function calculation

double sine(double arg, double lut[TABLE_SIZE]){
double retval;
double rem;
uint8_t index;

if(arg >= 0 && arg <= PI/2){
    // first quadrant
    index = arg/STEP_SIZE;
    rem = arg - index*STEP_SIZE;
    if(rem > 0){
        // sine value for given argument isn't directly in the lut
        if(index == (TABLE_SIZE-1)){
            // last point in the lut so the interval for the interpolation
            // is the last interval bounded with the index-1 and index
            index -= 1;
        }
        retval = (lut[index+1] - lut[index])/STEP_SIZE*rem + lut[index];
    }else{
        // sine value for given argument is directly in the lut
        retval = lut[index];
    }
}else if(arg > PI/2 && arg <= PI){
    // second quadrant
    index = (PI - arg)/STEP_SIZE;
    rem = (PI - arg) - index*STEP_SIZE;
    if(rem > 0){
        if(index == (TABLE_SIZE-1)){
            index -= 1;
        }
        retval = (lut[index+1] - lut[index])/STEP_SIZE*rem + lut[index];
    }else{
        retval = lut[index];
    }
    
}else if(arg > PI && arg <= 3*PI/2){
    // third quadrant
    index = (arg - PI)/STEP_SIZE;
    rem = (arg - PI) - index*STEP_SIZE;
    if(rem > 0){
        if(index == (TABLE_SIZE-1)){
            index -= 1;
        }
        retval = (-lut[index+1] + lut[index])/STEP_SIZE*rem - lut[index];
    }else{
        retval = -lut[index];
    }
    
}else{
    // fourth quadrant
    index = (2*PI - arg)/STEP_SIZE;
    rem = (2*PI - arg) - index*STEP_SIZE;
    if(rem > 0){
        if(index == (TABLE_SIZE-1)){
            index -= 1;
        }
        retval = (-lut[index+1] + lut[index])/STEP_SIZE*rem - lut[index];
    }else{
        retval = -lut[index];
    }
}
return retval;
}

Cosine function

double cosine(double arg, double lut[TABLE_SIZE]){
 double temp;

 temp = (arg + PI/2);
 if(temp > 2*PI){
    temp -= 2*PI;
 }
 return sine(temp, lut);
}
\$\endgroup\$
9
  • 8
    \$\begingroup\$ What is the target processor? \$\endgroup\$
    – Edward
    Apr 14, 2021 at 9:52
  • 6
    \$\begingroup\$ math.h functions are very likely using table look-up too, so what makes you think they are non-deterministic? \$\endgroup\$
    – Lundin
    Apr 14, 2021 at 11:24
  • 8
    \$\begingroup\$ I'm curious what the benchmark result showed. How fast is your code compared to the standard functions? \$\endgroup\$ Apr 14, 2021 at 16:33
  • 14
    \$\begingroup\$ #define PI 3.14 I LOLed. \$\endgroup\$ Apr 14, 2021 at 20:15
  • 4
    \$\begingroup\$ FYI this is how people programmed 3d graphics on the x86 processors before floating point coprocessors became common. So there's code like this all over the web to compare. \$\endgroup\$
    – don bright
    Apr 15, 2021 at 5:07

5 Answers 5

24
\$\begingroup\$

Your implementation is going to be slow, and the excuse "I need it to take a fixed amount of time" does not justify this. Using plain tables smells like cargo culting as well.

So I'm not tackling what you did wrong in your code, but what you did wrong in even thinking about your implementation. First, google how to implement cosine on a micro controller, or how to approximate cosine in general. Absolutely do not try to google "How to implement X with Y method" because you are only going to get that method (especially not at first). This will expand what you thought possible, for example, one much faster way that is much simpler to implement is using Chebyshev approximation. You would generate the polynomials before hand and compile them in depending on how accurate you want your application to be. You would also take advantage of symmetry to only approximate for every 45 degrees. Chebyshev approximation takes an easily quantifiable fixed amount of time.

Additionally, you may not want to use PI internally. Use base pi or 2.0 * pi. What I mean is, when you generate you Chebyshev polynomials, do it so that your input is normalized 1.0 -> -1.0 or some other normalized range. This can help with precision issues (as pi is irrational and can't be completely represented anyway), and allows you to use functions like cos2pi and return more accurate results if the user can make their angle easily in terms of normalized angle.

Another method that is much faster is actually using a real numerical approximation specific to the function ala https://stackoverflow.com/a/28050328/. No branches see? And there are further refinements of methods like this that can give even more accurate answers at the cost of speed. In general these methods outclass the Chebyshev formulation in both speed and precision, but the Chebyshev formulation doesn't require an analytical or advanced application specific numerical formulation to work, it just needs the input and output values, so if you have a more complicated function you only care about for a specific range of values, it can often be used to approximate the whole thing.

Because of the way you approached your problem, you ignored two massively better solutions to your problem, deterministic numerical computation and more advanced deterministic polynomial based approximations.

\$\endgroup\$
5
  • 1
    \$\begingroup\$ Very good points. Also, because Chebyshev approximations of even functions such as cosine only use even order polynomials, one can use Horner’s rule to reduce the number of multiplications to very few. \$\endgroup\$
    – Edward
    Apr 15, 2021 at 0:33
  • 1
    \$\begingroup\$ For a given precision, you can lop off one or two terms, with suitable adjustment of the Chebyshev coefficients. \$\endgroup\$
    – Rick James
    Apr 15, 2021 at 20:52
  • \$\begingroup\$ stackoverflow.com/a/28050328 useful for range -pi to pi. OP apparently wants 0...2pi. So direct usage does not meet the goal without changes (and maybe an if) but does serve to demo a simple approximation. \$\endgroup\$ Apr 17, 2021 at 22:27
  • \$\begingroup\$ @chux-ReinstateMonica This is simply not a problem... trig functions repeat... \$\endgroup\$
    – Krupip
    Apr 18, 2021 at 14:06
  • \$\begingroup\$ Yes trig functions repeat, but the linked code does not. \$\endgroup\$ Apr 18, 2021 at 14:50
4
\$\begingroup\$

One remark is that you should get rid of all the needless branching and code repetition. It's bad for performance and code maintenance both.

Given an angle you should be able to:

  • Take it's absolute value.
  • Divide by PI/2.
  • Convert to unsigned integer, truncating decimals.

Then you'll either have an index from 0 to 3 or you started with an angle larger than 2*PI (and might have to deal with that first, if that's a valid use case?).

Pseudo code:

double angle = ...;
angle /= PI/2.0;
unsigned long q = (unsigned long)angle;
// optionally some error handling/sanity checking here

quadrant(q, angle);

Where the quadrant function can use the variable q as an index to various other tables containing mulitpliers etc necessary for your calculation. It should boil down to the 5-6 lines that you currently have inside each if-else if.

\$\endgroup\$
1
4
\$\begingroup\$
  • #define PI 3.14... could use a few more digits!
  • sine should work for numbers greater than 2.0 * PI.
  • sine should work for negative numbers.
  • same for cosine.
  • if(temp > 2*PI) { temp -= 2*PI; } is ineffective for numbers greater than 4.0 * PI.

if(rem > 0){
    // sine value for given argument isn't directly in the lut
    if(index == (TABLE_SIZE-1)){
        // last point in the lut so the interval for the interpolation
        // is the last interval bounded with the index-1 and index
        index -= 1;
    }
    retval = (lut[index+1] - lut[index])/STEP_SIZE*rem + lut[index];
}else{
    // sine value for given argument is directly in the lut
    retval = lut[index];
}

Is it actually possible to have a remainder if index == TABLE_SIZE-1?

Maybe it would be better to always do the interpolation? We could change the table to double sinLUT[TABLE_SIZE + 1] = ..., and add an extra 1.0 value to the end. Then we always have an index+1 available (note: we'd have to change the index type from uint8_t).


There's a lot of duplicate code here that could be eliminated. We could do that by checking the quadrant, and flipping the argument direction / return value sign as necessary.

I'd be inclined to do something like this:

#include <assert.h>
#include <math.h>
#include <stdbool.h>
#include <stdint.h>

#define TABLE_SIZE 256
#define PI         3.1415926535897932384626433832795
#define STEP_SIZE  (1.0 / (TABLE_SIZE - 1))

// look-up table containing values of sine function over one quarter of period
// angle step size is pi/(2*256) i.e. pi/512 approximately 0.35°
double sinLUT[TABLE_SIZE + 1] = 
{
    0.000000,  0.006160,  0.012320,  0.018479,  0.024637,  0.030795,  0.036951,  0.043107, 
    0.049260,  0.055411,  0.061561,  0.067708,  0.073853,  0.079994,  0.086133,  0.092268, 
    0.098400,  0.104528,  0.110653,  0.116773,  0.122888,  0.128999,  0.135105,  0.141206, 
    0.147302,  0.153392,  0.159476,  0.165554,  0.171626,  0.177691,  0.183750,  0.189801, 
    0.195845,  0.201882,  0.207912,  0.213933,  0.219946,  0.225951,  0.231948,  0.237935, 
    0.243914,  0.249883,  0.255843,  0.261793,  0.267733,  0.273663,  0.279583,  0.285492, 
    0.291390,  0.297277,  0.303153,  0.309017,  0.314870,  0.320710,  0.326539,  0.332355, 
    0.338158,  0.343949,  0.349727,  0.355491,  0.361242,  0.366979,  0.372702,  0.378411, 
    0.384106,  0.389786,  0.395451,  0.401102,  0.406737,  0.412356,  0.417960,  0.423549, 
    0.429121,  0.434676,  0.440216,  0.445738,  0.451244,  0.456733,  0.462204,  0.467658, 
    0.473094,  0.478512,  0.483911,  0.489293,  0.494656,  0.500000,  0.505325,  0.510631, 
    0.515918,  0.521185,  0.526432,  0.531659,  0.536867,  0.542053,  0.547220,  0.552365, 
    0.557489,  0.562593,  0.567675,  0.572735,  0.577774,  0.582791,  0.587785,  0.592758, 
    0.597707,  0.602635,  0.607539,  0.612420,  0.617278,  0.622113,  0.626924,  0.631711, 
    0.636474,  0.641213,  0.645928,  0.650618,  0.655284,  0.659925,  0.664540,  0.669131, 
    0.673696,  0.678235,  0.682749,  0.687237,  0.691698,  0.696134,  0.700543,  0.704926, 
    0.709281,  0.713610,  0.717912,  0.722186,  0.726434,  0.730653,  0.734845,  0.739009, 
    0.743145,  0.747253,  0.751332,  0.755383,  0.759405,  0.763398,  0.767363,  0.771298, 
    0.775204,  0.779081,  0.782928,  0.786745,  0.790532,  0.794290,  0.798017,  0.801714, 
    0.805381,  0.809017,  0.812622,  0.816197,  0.819740,  0.823253,  0.826734,  0.830184, 
    0.833602,  0.836989,  0.840344,  0.843667,  0.846958,  0.850217,  0.853444,  0.856638, 
    0.859800,  0.862929,  0.866025,  0.869089,  0.872120,  0.875117,  0.878081,  0.881012, 
    0.883910,  0.886774,  0.889604,  0.892401,  0.895163,  0.897892,  0.900587,  0.903247, 
    0.905873,  0.908465,  0.911023,  0.913545,  0.916034,  0.918487,  0.920906,  0.923289, 
    0.925638,  0.927951,  0.930229,  0.932472,  0.934680,  0.936852,  0.938988,  0.941089, 
    0.943154,  0.945184,  0.947177,  0.949135,  0.951057,  0.952942,  0.954791,  0.956604, 
    0.958381,  0.960122,  0.961826,  0.963493,  0.965124,  0.966718,  0.968276,  0.969797, 
    0.971281,  0.972728,  0.974139,  0.975512,  0.976848,  0.978148,  0.979410,  0.980635, 
    0.981823,  0.982973,  0.984086,  0.985162,  0.986201,  0.987202,  0.988165,  0.989092, 
    0.989980,  0.990831,  0.991645,  0.992421,  0.993159,  0.993859,  0.994522,  0.995147, 
    0.995734,  0.996284,  0.996795,  0.997269,  0.997705,  0.998103,  0.998464,  0.998786, 
    0.999070,  0.999317,  0.999526,  0.999696,  0.999829,  0.999924,  0.999981,  1.000000, 1.000000
};

double fmodn(double x, double y)
{
    return x - y * floor(x/y);
}

double mix(double x, double y, double a) // (lerp)
{
    return x * a + y * (1.0 - a);
}

double sine(double arg)
{
    arg = fmodn(arg, 2.0 * PI); // range 0.0 to 2.0 * PI

    const int quadrant = (int)(arg / (0.5 * PI)); // range 0 to 3
    const bool flip = (quadrant > 2); // -retval
    const bool invert = (quadrant % 2); // 1.0 - arg

    arg = fmod(arg, 0.5 * PI) / (0.5 * PI); // range 0.0 to 1.0
    arg = invert ? 1.0 - arg : arg;

    const int index = arg / STEP_SIZE;
    const double rem = arg - index * STEP_SIZE;
    double retval = mix(sinLUT[index], sinLUT[index + 1], rem / STEP_SIZE);
    
    retval = flip ? -retval : retval;

    return retval;
}

#include <stdio.h>

int main()
{
    for (double i = -2.5 * PI; i < 2.5 * PI; i += PI * 0.05)
        printf("%f %f\n", sin(i), sine(i));
}

This version maps arg to a range 0.0 to 1.0 before calculating the index in the table (and changes STEP_SIZE to match). It's a bit slower to do that though, and not really necessary.


The standard library may not guarantee a "fixed" run-time, but that doesn't mean a specific platform implementation won't suit your needs fine.

I'm not too familiar with this stuff; I'm sure there are much faster and more accurate algorithms out there...

\$\endgroup\$
2
  • \$\begingroup\$ For constant time, you need to skip the test for temp > 2*pi. Instead always do mod. \$\endgroup\$
    – Rick James
    Apr 15, 2021 at 20:53
  • \$\begingroup\$ @RickJames "Instead always do mod" --> If you are thinking of fmod(), that function is not specified to be constant time. If otherwise, how are you suggesting to "do mod"? \$\endgroup\$ Apr 17, 2021 at 18:13
2
\$\begingroup\$

In my testing (using C++ std::sin(), which should be the same function) the Standard Library sin() is about 9 times faster than the float version of LUT sine shown in @user673679's answer.

If approximation suites you, as suggested already, You might find some help in Faster Math Functions (part 2) by Robin Green (part 1 also available).

Here are couple approximation examples: https://www.desmos.com/calculator/dm5wdeqjy0

Can you take advantage of SSE/AVX?

EDIT: variation data (std::sin and Chebychev approximation found in linked Desmos sheet):

Approx. rdtsc/val.... datasize: 524288
std::sin         : 11.0862
ChebySineApprox    : 1.74816

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9639
ChebySineApprox    : 1.74816

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9639
ChebySineApprox    : 1.74799

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9639
ChebySineApprox    : 1.74799

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9634
ChebySineApprox    : 1.74786

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9607
ChebySineApprox    : 1.74786

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9606
ChebySineApprox    : 1.74786

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9603
ChebySineApprox    : 1.74786

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9603
ChebySineApprox    : 1.74786

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9603
ChebySineApprox    : 1.74786

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9603
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9603
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9603
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9603
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74775

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9599
ChebySineApprox    : 1.74771

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.7477

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9593
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743

Approx. rdtsc/val.... datasize: 524288
std::sin         : 10.9519
ChebySineApprox    : 1.74743
\$\endgroup\$
3
  • 4
    \$\begingroup\$ Your testing was probably done on hardware with strong floating-point support, and might not be valid on the target platform. However, the question is strangely silent on what the target actually is, or any benchmarking done, and there's not much you can do about that (it's already been asked in comments)... \$\endgroup\$ Apr 15, 2021 at 6:38
  • 1
    \$\begingroup\$ While you were testing, how much variation of timing of sin() was there? \$\endgroup\$
    – Rick James
    Apr 15, 2021 at 20:54
  • 1
    \$\begingroup\$ Added data from my test which may tell something related to variation. \$\endgroup\$
    – Juha P
    Apr 16, 2021 at 11:05
1
\$\begingroup\$

A key missing specification is the precision needed for the result as that steers algorithm design.

Going forward that arg may be wide ranging and table look-up precision is good enough: I'd recommend simplifications:

  1. Pass the table in via global unless there might be more than 1 table.

  2. Perform as much as possible with integer math.

  3. With such low precision requirements, consider a float table.

  4. Table size should be a power-of-2 (which OP is doing so far) to simplify "modding" the angle.

  5. To achieve some level of "fixed calculation time", make the table 2*pi wide.

  6. I suspect a simple double arg to integer and table look-up will provide the fastest solution, with somewhat fixed time - although at a cost of space and precision.

Simplification:

#define SINE_TABLE_SIZE 256
#define MY_PI 3.1415926535897932384626433832795
#define RAD2INDEX (SINE_TABLE_SIZE / (2 * MY_PI))
#define INDEX2RAD ((2 * MY_PI) / SINE_TABLE_SIZE)

static const float sinLUT[SINE_TABLE_SIZE + 1] = { 
  // pre-compute these and code to 9+ significant digits.
  sinf(0 * INDEX2RAD),
  sinf(1 * INDEX2RAD),
  //  ...
  sinf(SINE_TABLE_SIZE * INDEX2RAD)
};

// Nice simple functions, both use same table

float my_sine(double arg) {
  long index = lround(arg * RAD2INDEX);
  unsigned i = (unsigned) (index & (SINE_TABLE_SIZE - 1));  // "mod"
  return sinLUT[i];
}

float my_cosine(double arg) {
  long index = lround(arg * RAD2INDEX);  // Only floating point math
  index += SINE_TABLE_SIZE/4;
  unsigned i = (unsigned) (index & (SINE_TABLE_SIZE - 1));
  return sinLUT[i];
}

Note: When arg < 0, answer relies on 2's complement long - something very common these days.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.