Sine function calculated via look-up table and linear interpolation

Question

I have been attempting to implement a function calculating values of the sine function. I know that there are several similar threads regarding this topic but my goal was to try to implement such function in my own way as an excercise. I suppose that the argument will be in range \$\left<-2\pi\, +2\pi\right>\$.

In the time being I have below given code in C++ language which seems to be working (I have compared the outputs of that code with the sine values calculated by the Excel).

#define PI         3.14
#define TABLE_SIZE 256
#define STEP_SIZE (2*PI/255)

double lut[TABLE_SIZE] = {
    0.0000, 0.0245, 0.0490, 0.0735, 0.0980, 0.1223, 0.1467, 0.1709, 
    0.1950, 0.2190, 0.2429, 0.2666, 0.2901, 0.3135, 0.3367, 0.3597, 
    0.3825, 0.4050, 0.4274, 0.4494, 0.4712, 0.4927, 0.5139, 0.5348, 
    0.5553, 0.5756, 0.5954, 0.6150, 0.6341, 0.6529, 0.6713, 0.6893, 
    0.7068, 0.7240, 0.7407, 0.7569, 0.7727, 0.7881, 0.8029, 0.8173, 
    0.8312, 0.8446, 0.8575, 0.8698, 0.8817, 0.8930, 0.9037, 0.9140, 
    0.9237, 0.9328, 0.9413, 0.9493, 0.9568, 0.9636, 0.9699, 0.9756, 
    0.9806, 0.9852, 0.9891, 0.9924, 0.9951, 0.9972, 0.9988, 0.9997, 
    1.0000, 0.9997, 0.9988, 0.9974, 0.9953, 0.9926, 0.9893, 0.9854, 
    0.9810, 0.9759, 0.9703, 0.9640, 0.9572, 0.9498, 0.9419, 0.9333, 
    0.9243, 0.9146, 0.9044, 0.8937, 0.8824, 0.8706, 0.8583, 0.8454, 
    0.8321, 0.8182, 0.8039, 0.7890, 0.7737, 0.7580, 0.7417, 0.7251, 
    0.7080, 0.6904, 0.6725, 0.6541, 0.6354, 0.6162, 0.5967, 0.5769, 
    0.5566, 0.5361, 0.5152, 0.4941, 0.4726, 0.4508, 0.4288, 0.4065, 
    0.3840, 0.3612, 0.3382, 0.3150, 0.2917, 0.2681, 0.2444, 0.2205, 
    0.1966, 0.1724, 0.1482, 0.1239, 0.0996, 0.0751, 0.0506, 0.0261, 
    0.0016, -0.0229, -0.0475, -0.0719, -0.0964, -0.1208, -0.1451, -0.1693, 
    -0.1934, -0.2174, -0.2413, -0.2650, -0.2886, -0.3120, -0.3352, -0.3582, 
    -0.3810, -0.4036, -0.4259, -0.4480, -0.4698, -0.4913, -0.5125, -0.5334, 
    -0.5540, -0.5743, -0.5942, -0.6137, -0.6329, -0.6517, -0.6701, -0.6881, 
    -0.7057, -0.7229, -0.7396, -0.7559, -0.7717, -0.7871, -0.8020, -0.8164, 
    -0.8303, -0.8437, -0.8566, -0.8690, -0.8809, -0.8923, -0.9031, -0.9133, 
    -0.9230, -0.9322, -0.9408, -0.9488, -0.9563, -0.9632, -0.9695, -0.9752, 
    -0.9803, -0.9849, -0.9888, -0.9922, -0.9950, -0.9971, -0.9987, -0.9996, 
    -1.0000, -0.9998, -0.9989, -0.9975, -0.9954, -0.9928, -0.9895, -0.9857, 
    -0.9813, -0.9762, -0.9706, -0.9644, -0.9577, -0.9503, -0.9424, -0.9339, 
    -0.9249, -0.9153, -0.9051, -0.8944, -0.8832, -0.8714, -0.8591, -0.8463, 
    -0.8330, -0.8191, -0.8048, -0.7900, -0.7747, -0.7590, -0.7428, -0.7262, 
    -0.7091, -0.6916, -0.6736, -0.6553, -0.6366, -0.6175, -0.5980, -0.5782, 
    -0.5580, -0.5374, -0.5166, -0.4954, -0.4740, -0.4522, -0.4302, -0.4080, 
    -0.3854, -0.3627, -0.3397, -0.3166, -0.2932, -0.2696, -0.2459, -0.2221, 
    -0.1981, -0.1740, -0.1498, -0.1255, -0.1011, -0.0767, -0.0522, -0.0277
};

double sine(double x, double lut[TABLE_SIZE])
{
  bool negateTableValue = false;
  if (x < 0) {
    // sin(-x) = -sin(x)
    x = -x;
    negateTableValue = true;
  }
  
  uint8_t index_01 = x/STEP_SIZE;
  uint8_t index_02 = (index_01 + 1);
      
  double aux = (lut[index_02] - lut[index_01])/STEP_SIZE*(x - index_01*STEP_SIZE) + lut[index_01];
  
  if (negateTableValue) {
    return -aux;
  } else {
    return aux;
  } 
}

The sine values calculation is based on the look-up table containing the pre-computed values of the sine function covering the whole period \$\left<0, 2\pi\right>\$ with 256 values. I have decided to use the linear interpolation method for improving the precision.

I have one doubt regarding the linear interpolation. Namely I have been using table with 256 entries but most of the solutions exploiting the linear interpolation use look-up table with one additional entry. I would say that it isn't necessary in my case because the index variables are uint8_t type i.e. they can store values from range 0-255. But I would like to know other ones opinion. Thank you in advance.

Please do not edit the question, especially the code after an answer has been posted. Everyone needs to be able to see what the reviewer was referring to. What to do after the question has been answered. — pacmaninbw, Jul 14, 2021 at 23:08
@pacmaninbw but the edit you rolled back was not editing the code in the question? — the default., Jul 15, 2021 at 13:23
@thedefault. Any and all answer invalidation (AI) is not allowed. You should note pacmaninbw only emphasized "the code", as "the code" is normally how OPs AI. But pointed out any AI in any part of the question is not allowed; "Please do not edit the question". — Peilonrayz, Jul 16, 2021 at 22:31
" my goal was to try to implement such function in my own way as an excercise" --> Why this approach? Speed? code size? Other approaches are better for both of those. Perhaps it is only an exercise? — chux - Reinstate Monica, Jul 20, 2021 at 19:11

G. Sliepen · Accepted Answer · 2021-07-14 20:22:08Z

You only need to store values for \$\left[0, \frac{1}{2}\pi\right)\$

You can exploit the structure of the sine function, and only store values in the table for inputs between 0 and \$\frac{1}{2}\pi\$. Values between \$\pi\$ and \$2\pi\$ are just the negative of the values between \$0\$ and \$\pi\$, and values left and right of \$\frac{1}{2}\pi\$ are mirrored. The upshot is that keeping the table size the same, you get four times the resolution.

Use a proper value for \$\pi\$

Don't guess or use a low-precision estimate of \$\pi\$. There are several ways to get a good value for it. A simple one is:

static const double PI = std::atan(1) * 4;

In C++20, you can use std::numbers::pi.

The step size is slightly wrong

And not just because of the poor approximation of \$\pi\$, but also because the values in your table are taken at intervals of \$\pi/256\$, so STEP_SIZE should match that.

Prefer `static constexpr` over `#define`

Since C++11, you can use static constexpr to declare constants, and you can even use integers declared as such to specify the size of arrays.

Note that since std::atan() is not constexpr itself, you can't use it to initialize a constexpr variable, but you can make PI static const instead, and most compilers will optimize it anyway. Of course, being able to use std::numbers:pi is even better.

Avoid hard-coding the table

You don't need to write all the values by hand, you can let the compiler generate the table for you at compile time, like so:

#include <array>
#include <cmath>

static constexpr size_t TABLE_SIZE = 256;
static const double PI = std::atan(1) * 4;

static std::array<double, TABLE_SIZE> make_lut() {
    std::array<double, TABLE_SIZE> lut{};

    for(size_t i = 0; i < TABLE_SIZE; ++i) {
        lut[i] = std::sin(i * PI / 2 / TABLE_SIZE);
    }

    return lut;
}

std::array<double, TABLE_SIZE> lut = make_lut();

Interpolation

If you can use C++20, then you can use std::lerp() to interpolate between two numbers. However, your implementation is fine since it interpolates between two numbers between -1 and 1, and not much can go wrong there.

Using an additional entry

The reason for an additional entry is that you only have to calculate index_01, and then you read two consecutive values from the table. This is slightly more efficient: no need for you or the compiler to handle wrapping of index_02, and reading memory locations sequentially might be faster than reading two locations that are far apart.

Your code is fine, I don't think it will make a significant difference to add the additional entry. However, as soon as you want to change the size of the table to something that is not a power of two, it will pay off (modulo operations are basically divisions, and those are expensive).

BTW, you could always use const integers to specify the size of an array, back to Cfront 1.x in the 80's. — JDługosz, Jul 14, 2021 at 21:19
@G. Sliepen thank you very much for your review. I have just added some notes regarding the step size selection to my question. May I ask you for your opinion? — L3sek, Jul 14, 2021 at 21:40
Using 256 for STEP_SIZE is fine. If the input is 2π, then index_01 is 0 and index_02 is 1, but the interpolation should give the first index a weight of 1 and the second 0. I think x - index_01 * STEP_SIZE is wrong though; what if x is a much larger number than 2π? — G. Sliepen, Jul 14, 2021 at 21:52
I suppose that the arguments will be in range \$\left<0, 2\pi\right>\$. I am sorry that I didn't mention that information in my question. — L3sek, Jul 15, 2021 at 5:01
@L3sek You already transform arguments from [-2π,2π] so that they end up in [0,2π]. G.Sliepen is suggesting a few more transformations so that your values end up in [0,π/2], allowing you to have four times as many values in the table without needing to increase its space. — Teepeemm, Jul 15, 2021 at 22:04

chux - Reinstate Monica · Accepted Answer · 2021-07-15 03:09:33Z

double vs. float

Given the low precision OP's linear interpolation approach, double lut[TABLE_SIZE] could use float lut[TABLE_SIZE]. Or perhaps use a 2x TABLE_SIZE to maintain the same code space, but gain accuracy?

-0.0

To provide sine(-0.0, ...) --> -0.0, use signbit test rather than a x < 0 test.

3.14 vs 3.1415926535897932384626433832795

Use a better machine pi.
Recommend at least the precision of double.

Magic number

Why 255 in #define STEP_SIZE (2*PI/255)? I'd expect #define STEP_SIZE (2*PI/(TABLE_SIZE-1))

Reduce accumulated error

Rather than 255 in (2*PI/255), consider a power-of-two to not introduce additional calculation error.

times vs divide

Consider multiplication is speedier.

// uint8_t index_01 = x/STEP_SIZE;
uint8_t index_01 = x*STEP_SIZE_INVERSE;

L3sek · Accepted Answer · 2021-07-15 06:32:34Z

The STEP_SIZE value is tightly coupled with my doubt. In case the table covers \$\left<0, 2\pi\right>\$ with 256 values it would make more sense to have #define STEP_SIZE (2*PI/256) (because of 256 values divide unit circle into 256 segments of \$\frac{2\pi}{256}\,\mathrm{rad}\$ each) instead of #define STEP_SIZE (2*PI/255).

The problem which I have with #define STEP_SIZE (2*PI/256) is that in case x is 2*PI then uint8_t index_01 overflows to 0 and then uint8_t index_02 contains 1. So the interpolation takes place between wrong look-up table entries. Let's say I would use uint16_t for index_01 and index_02 instead of uint8_t I would have for 2*PI argument index_01 with value 256 (so I would need one additional entry in the look-up table) and then I would need to somehow force wrapping up the index_02 containing 257 to 0 which I could probably do in this manner:

#define TABLE_SIZE 257
#define STEP_SIZE (2*PI/256)

double lut[TABLE_SIZE] = {
    0.0000, 0.0245, 0.0490, 0.0735, 0.0980, 0.1223, 0.1467, 0.1709, 
    0.1950, 0.2190, 0.2429, 0.2666, 0.2901, 0.3135, 0.3367, 0.3597, 
    0.3825, 0.4050, 0.4274, 0.4494, 0.4712, 0.4927, 0.5139, 0.5348, 
    0.5553, 0.5756, 0.5954, 0.6150, 0.6341, 0.6529, 0.6713, 0.6893, 
    0.7068, 0.7240, 0.7407, 0.7569, 0.7727, 0.7881, 0.8029, 0.8173, 
    0.8312, 0.8446, 0.8575, 0.8698, 0.8817, 0.8930, 0.9037, 0.9140, 
    0.9237, 0.9328, 0.9413, 0.9493, 0.9568, 0.9636, 0.9699, 0.9756, 
    0.9806, 0.9852, 0.9891, 0.9924, 0.9951, 0.9972, 0.9988, 0.9997, 
    1.0000, 0.9997, 0.9988, 0.9974, 0.9953, 0.9926, 0.9893, 0.9854, 
    0.9810, 0.9759, 0.9703, 0.9640, 0.9572, 0.9498, 0.9419, 0.9333, 
    0.9243, 0.9146, 0.9044, 0.8937, 0.8824, 0.8706, 0.8583, 0.8454, 
    0.8321, 0.8182, 0.8039, 0.7890, 0.7737, 0.7580, 0.7417, 0.7251, 
    0.7080, 0.6904, 0.6725, 0.6541, 0.6354, 0.6162, 0.5967, 0.5769, 
    0.5566, 0.5361, 0.5152, 0.4941, 0.4726, 0.4508, 0.4288, 0.4065, 
    0.3840, 0.3612, 0.3382, 0.3150, 0.2917, 0.2681, 0.2444, 0.2205, 
    0.1966, 0.1724, 0.1482, 0.1239, 0.0996, 0.0751, 0.0506, 0.0261, 
    0.0016, -0.0229, -0.0475, -0.0719, -0.0964, -0.1208, -0.1451, -0.1693, 
    -0.1934, -0.2174, -0.2413, -0.2650, -0.2886, -0.3120, -0.3352, -0.3582, 
    -0.3810, -0.4036, -0.4259, -0.4480, -0.4698, -0.4913, -0.5125, -0.5334, 
    -0.5540, -0.5743, -0.5942, -0.6137, -0.6329, -0.6517, -0.6701, -0.6881, 
    -0.7057, -0.7229, -0.7396, -0.7559, -0.7717, -0.7871, -0.8020, -0.8164, 
    -0.8303, -0.8437, -0.8566, -0.8690, -0.8809, -0.8923, -0.9031, -0.9133, 
    -0.9230, -0.9322, -0.9408, -0.9488, -0.9563, -0.9632, -0.9695, -0.9752, 
    -0.9803, -0.9849, -0.9888, -0.9922, -0.9950, -0.9971, -0.9987, -0.9996, 
    -1.0000, -0.9998, -0.9989, -0.9975, -0.9954, -0.9928, -0.9895, -0.9857, 
    -0.9813, -0.9762, -0.9706, -0.9644, -0.9577, -0.9503, -0.9424, -0.9339, 
    -0.9249, -0.9153, -0.9051, -0.8944, -0.8832, -0.8714, -0.8591, -0.8463, 
    -0.8330, -0.8191, -0.8048, -0.7900, -0.7747, -0.7590, -0.7428, -0.7262, 
    -0.7091, -0.6916, -0.6736, -0.6553, -0.6366, -0.6175, -0.5980, -0.5782, 
    -0.5580, -0.5374, -0.5166, -0.4954, -0.4740, -0.4522, -0.4302, -0.4080, 
    -0.3854, -0.3627, -0.3397, -0.3166, -0.2932, -0.2696, -0.2459, -0.2221, 
    -0.1981, -0.1740, -0.1498, -0.1255, -0.1011, -0.0767, -0.0522, -0.0277,
     0.0000
};
    
uint16_t index_01 = x/STEP_SIZE;
uint16_t index_02 = index_01 + 1;
if (index_02 == TABLE_SIZE) {
index_02 = 0;
}

The reason for choosing #define STEP_SIZE (2*PI/255) was to avoid the if statement for wrapping-up the index_02. But when I am writting this edit I realized that the second variant is probably better because it doesn't bring the error and it is more comprehensible.

JDługosz · Accepted Answer · 2021-07-15 22:24:40Z

One thing to add to what's been said already:

if (negateTableValue) {
    return -aux;
  } else {
    return aux;
  }

There's a lot of cognitive overhead to read this, when its meaning is very simple. As such, it's much clearer to write as:

return negateTableValue? -aux : aux;

You're not doing two arbitrary different things; you are returning, period. Pushing the condition into the return value itself takes it out of the main-line level of reading the code.

Meanwhile, I note that you wrote this as an exercise... but I wonder, what might you really use this for?

I imagine a situation like a small microcontroller that doesn't have floating-point hardware, and it's being used to decide how to move a robot limb or something so the accuracy doesn't need to be a huge number of significant digits.

As such:

make the domain type a template argument.

You might be using your own fixed-point type, or your own take on "floating" point that only has a few different ranges it works in rather than general exponent, taylored to the needed precision at small angles vs large angles of the physical model.

make the table size a (non-type) template argument.

Likewise, how much precision you actually need can be specified as needed for that particular project. If you are generating the table, using constexpr code, this is especially useful. If you are to explicitly give the table as you do in the original post, then use the actual number of values present to pick up on the table size.

make the domain configurable

If you only have small angle deflections, you don't need to waste space in the table for the entire quarter-circle. Maybe it only needs to handle 10° max, e.g. you're modeling the movement of a clock pendulum. Or, perhaps the project only needs to work with angles between 35° and 55°. So, give it arguments to configure the range of angles it is able to handle.

how to package a configurable library

The simplest thing is to use a class template. The table and helper functions are private members of the class.

using TestSign1 = fast_sign_library<myfixed_t, 256, 0, 90>;  // used degrees so it can take ints easily.  That's independent of what sin(θ) actually takes.

myfixed_t angle= .......;
myfixed_t sn_a = TestSign1::sin(angle);

You can further wrap that in a small inline function to give it a simple name.

fishinear · Accepted Answer · 2021-07-16 15:43:39Z

Additional to what is already said, make explicit that you are rounding downwards with floor, and simplify the interpolation computation as a weighted average:

double index_x = x / STEP_SIZE;
uint8_t index_01 = floor(index_x);
uint8_t index_02 = index_01 + 1;

double aux = lut[index_01] * (index_02 - index_x) + lut[index_02] * (index_x - index_01);

Instead of a boolean negateTableValue you could use a double sign, to make it easier to read (for me, although that may be a personal preference):

double sign = 1;
if (x < 0) {
    sign = -1;
    x = -x;
}
...
return sign * aux;

You should also make sure it works for values outside the -2π..2π range, but I leave that up to you as an exercise.

I too was going to mention something about using fmod() to support the full range of input, but I forgot, so thank you for adding that. — Toby Speight, Jul 16, 2021 at 15:42

Stack Exchange Network

Sine function calculated via look-up table and linear interpolation

5 Answers 5

You only need to store values for \$\left[0, \frac{1}{2}\pi\right)\$

Use a proper value for \$\pi\$

The step size is slightly wrong

Prefer `static constexpr` over `#define`

Avoid hard-coding the table

Interpolation

Using an additional entry

make the domain type a template argument.

make the table size a (non-type) template argument.

make the domain configurable

how to package a configurable library

Your Answer

Not the answer you're looking for? Browse other questions tagged
c++
mathematics
lookup
interpolation
or ask your own question.

Hot Network Questions

Sine function calculated via look-up table and linear interpolation

5 Answers 5

You only need to store values for \$\left[0, \frac{1}{2}\pi\right)\$

Use a proper value for \$\pi\$

The step size is slightly wrong

Prefer static constexpr over #define

Avoid hard-coding the table

Interpolation

Using an additional entry

make the domain type a template argument.

make the table size a (non-type) template argument.

make the domain configurable

how to package a configurable library

Your Answer

Sign up or log in

Post as a guest

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

Related

Hot Network Questions

Prefer `static constexpr` over `#define`

Not the answer you're looking for? Browse other questions tagged
c++
mathematics
lookup
interpolation
or ask your own question.