Compact sin(x) which gets to within 0.1%, x = 0 to 2xPI

[peter-h]

Currently I am using standard sin(x) but math.h takes up about 6k (STM 32F417) and I would like to reduce this. There is a lot of stuff online but most of it uses stuff like pow(x,y) which also uses math.h so is dumb unless you just want more speed. My use of sin(x) is for table precomputation only and it doesn't matter how slow it is. Also a lot of code posted doesn't work

Does anyone have a compact C function which has actually been tested?

It is to generate a table of signed 16 bit ints, scaled so that +/- 30000 (approx) is the max value. So I don't need floats at all really.

One problem is that the table size (64/128/256/etc) is set at compile time, otherwise I could just precompute the table, and say a 128 sized table would be just 256 bytes and that would be it. Yes of course anything is possible; I could get somebody on freelancer.com, most likely in the former USSR, to knock up a prog for £50 which reads the .h file where the table size is defined and generates it, and then it can be #included, but then I have to make sure the chain of dependencies is correctly set up (I am using ST Cube IDE) etc. A float sin(x) would be a lot neater.

I am also aware that one needs to compute only the octant, and the other 7/8 of the circle can be obtained by swapping things around. I did this (the Horn algorithm) on a Z80 in the 1980s, but it was primitive and inaccurate - it was for drawing circles on a 512x512 display.

[Ian.M]

Try googling: CORDIC in C
Wikipedia has a decent page on the CORDIC algorithm itself.

[ataradov]

You also don't need to parse anything if the number of possible sizes is limited (and it would be limited by memory size, of course).  Just generate tables for all sizes and #include the one you need depending  on the define. All the code except for the tables themselves would be hand-written and would not require any special maintenance or external tools in the workflow.

And if you can drop the requirement to < 2%, then there is this absolutely trivial formula https://en.wikipedia.org/wiki/Bhaskara_I%27s_sine_approximation_formula

And as usual, Taylor series is a good approach if you just need something that is simple to implement, but not necessarily the fastest.

[T3sl4co1l]

I already gave the runtime function in the other thread, so calculating the table I guess is the remaining part?
https://www.seventransistorlabs.com/sintables.html
Pop open the console, adjust the obvious parameters and run the update function again for whatever you want.

Tim

[nfmax]

Have a look at how the Yamaha DX7 chip implemented a fast sine function (actually the logarithm of the sine) using the chip resources available in 1983. It may give you some ideas

http://www.righto.com/2021/12/yamaha-dx7-reverse-engineering-part-iii.html

[peter-h]

Thank you all for the suggestions.

I am being lazy and was looking for something that is usable straight off and works

It was CORDIC that I implemented in Z80 assembler. The version published by Mr Horn, c. 1976, was used in the French Thomson-CSF (the forerunner of ST that we know today) 9365 graphics chip.

I have also seen a sine generator which uses a recirculating shift register principle. It generates integer sin(x) to any desired precision, with a direct digital output, but the output is continuous; you can't just feed in any gixen x.

I suspect a Taylor implementation, with say 10 iterations, will do, but I have not found code online which looked like it actually works.

[emece67]

.

[brucehoult]

I'm not sure what the difficulty is here. Sine is rather easy to calculate, if it doesn't have to be fast or accurate to the last bit.

Off the top of my head ...

Code: [Select]

#include <stdio.h>

#define TBLSZ 64

#define PI 3.14159265358979323846
#define SCALE 32767

int my_sin(int i){
  double rad = 2*PI*i/TBLSZ, term=rad, sin=rad, old_sin=0;
  for (int n=3; sin != old_sin; n+=2){
    old_sin = sin;
    term *= (rad*rad) / (-n * (n-1));
    sin += term;
  }
  return (sin * SCALE) + (sin >= 0 ? 0.5 : -0.5);
}

int main(){
  printf("short sintbl[] = {\n");
  for (int i=0; i<TBLSZ; ++i)
    printf("  %6d,\n", my_sin(i));
  printf("};\n");
  return 1;
}

As I understand, the 32F417 has an FPU so you don't have to worry about using fixed point to create your table.

[SiliconWizard]

Generating tables at compile time sounds trivial enough. Why would you need to hire someone to do it? It should be a few minutes work to write the utility?

[emece67]

.

[Siwastaja]

What, you want to custom implement sin() only to run at the init, generating a fixed LUT, only because the size of that LUT is a compile-time constant?

I fail to see the problem in just creating the tooling. You'd have done it in the same time it took to write that post. For example, a small C program that #includes the same header (or takes the same command line parameters) where the desired size is. Then run sin() in a loop, write results to a file in C array syntax. #include that generated file where needed, and do a safety check that the size matches (#if DESIRED_SIZE != sizeof sin_table  #error Workflow error). Just add that generator as a step in a makefile before the compilation of the file(s) that need this table. If Cube IDE makes it hard or impossible to modify makefile, then I'm afraid you have no other options than not use it. If some manager demands it, you can "keep using" it without actually using it.

[peter-h]

Sorry, guys, of course you are right. I should just generate the LUTs for 64/128/256/512/1024 values per cycle and #include them according to the SPC (samples per cycle) value which is #defined in the embedded code.

I just don't know how to write a win32 command line program. I should find out how, but right now I don't have a win32 command line C compiler on my PC. Last time I did that sort of thing was in the days of Borland C++ Or I would have used Basic, which came with every "IBM PC". I have found this: https://www.freebasic.net/.

Re "#error Workflow error" that is something else I wondered about: how does one generate warnings, or fatal compilation errors, in GCC. That's very useful.

[langwadt]

Sorry, guys, of course you are right. I should just generate the LUTs for 64/128/256/512/1024 values per cycle and #include them according to the SPC (samples per cycle) value which is #defined in the embedded code.

I just don't know how to write a win32 command line program. I should find out how, but right now I don't have a win32 command line C compiler on my PC. Last time I did that sort of thing was in the days of Borland C++ Or I would have used Basic, which came with every "IBM PC". I have found this: https://www.freebasic.net/.

Re "#error Workflow error" that is something else I wondered about: how does one generate warnings, or fatal compilation errors, in GCC. That's very useful.

if you can live with the dirty feeling, you could just use a spreadsheet ...

[ledtester]

What, you want to custom implement sin() only to run at the init, generating a fixed LUT, only because the size of that LUT is a compile-time constant?

The goal was to create a constexpr function for the sin table, and not all C++ compilers allow you to use the math functions, sin(), cos(), sqrt(), etc., in constexpr expressions/functions.

GCC is an exception is this case, and has implemented the std math functions as constexpr (when it is possible):

https://stackoverflow.com/questions/27744079/is-it-a-conforming-compiler-extension-to-treat-non-constexpr-standard-library-fu

https://stackoverflow.com/questions/17347935/constexpr-math-functions

I just don't know how to write a win32 command line program. I should find out how, but right now I don't have a win32 command line C compiler on my PC. Last time I did that sort of thing was in the days of Borland C++ Or I would have used Basic, which came with every "IBM PC". I have found this: https://www.freebasic.net/.

For one-time compilation there are several online compilers available, e.g.:

https://www.onlinegdb.com/online_c_compiler

https://rextester.com/l/cpp_online_compiler_gcc

[Jan Audio]

I should just generate the LUTs for 64/128/256/512/1024 values per cycle

Ehh ?, dont double things many times, just shift/skip bytes and use only 1 table.

My sine.h file includes all sizes lookuptable with the right define.
I add the file and say how big i want : 2048 maybe 4096 used without interpolation.
Can add to RAM also for better performance.

[Siwastaja]

I just don't know how to write a win32 command line program. I should find out how, but right now I don't have a win32 command line C compiler on my PC. Last time I did that sort of thing was in the days of Borland C++ Or I would have used Basic, which came with every "IBM PC".

MINGW is a quick install. Cross-compling between Windows/linux is a bit more tedious to set up, but luckily you don't need to do that.

I really prefer my tooling in C, not just because I like C, but because you can include shared headers with datatypes, constants and macros, even share functions or complete compilation units between PC/MCU side, etc. Of course need to be careful with data type widths (always use stdint), implementation defined areas, endianness etc., but with great responsibility comes great power.

If you have just those few different options and you are quite sure they don't change often, just create the files once by whatever means. Matlab, Octave, Python, spreadsheet, heck, Qbasic or how about 8088 assembly or maybe COBOL? The point is, create all the files once and include the right one automagically (#if LEN == 128 #include lut128.inc #elif LEN == 256... ... ... #elif #error Illegal LEN #endif)

[T3sl4co1l]

As you may've noticed, I have a tendency to write my tooling in JS... I mean, a browser is basically anywhere you are, in fact you're even using it right now!

Tim

[magic]

Do it in Excel and save as CSV
Then it's a simple substitution in Notepad to get C code.

[Ian.M]

Its easy enough to use worksheet formulae to get Excel to emit valid C code for an array as a column, then copy/paste that column into a source file.

[SiliconWizard]

I just don't know how to write a win32 command line program. I should find out how, but right now I don't have a win32 command line C compiler on my PC.[/url]

Any C or C++ compiler for Windows can generate command line programs. If you ever have Visual Studio, for instance, it absolutely can (and btw, its compilers are command line too!)
But you can use many other compilers (GCC, Clang, whatever...) available on Windows. Or you could use an interpreted language if you prefer/find it easier in your case. Some will suggest Python, I personally mostly use Lua instead, and can do this kind of things with a few lines of code.

[Jeroen3]

One problem is that the table size (64/128/256/etc) is set at compile time, otherwise I could just precompute the table, and say a 128 sized table would be just 256 bytes and that would be it. Yes of course anything is possible; I could get somebody on freelancer.com, most likely in the former USSR, to knock up a prog for £50 which reads the .h file where the table size is defined and generates it, and then it can be #included, but then I have to make sure the chain of dependencies is correctly set up (I am using ST Cube IDE) etc. A float sin(x) would be a lot neater.

If you generate all options as const tables, with a little bit if python code for example, and then add all of them to your project, but you only reference one of then, the rest will be removed by the linker.
You could even, if you wanted to, have the python generate the list based on what is in the #define.... if you add the python as a build-step.

Code: (python) [Select]

import math
def calculateTable(length):
    table = []
    for i in range(length):
        value = math.cos( 2 * math.pi * i / length );
        table.append(value)
    return table

n = 256
table = calculateTable(n)

f = open("table"+str(n)+".c", "w", encoding="ansi")
f.write("float table"+str(n)+"["+str(n)+"] = {\n")
for x in table:
    f.write("    "+str(x)+"f,\n")
f.write("};\n")
f.close()
I do that all the time. This is why python is amazing.

[peter-h]

I did it in Freebasic, in about 10 mins

Code: [Select]

dim twopi as double
dim x as double

twopi = 2*3.1415926
for x = 0 to twopi step twopi/1024
print x, sin(x)
next x
Now I just need to generate the C const format.

Took me a while to realise that real men don't use line numbers anymore

[Nominal Animal]

Here is what I'd use, assuming single-precision floating-point hardware is available, noting that
  sinscale(p, s) ≃ (int)round(s * sin(2 * PI * p))

Code: [Select]

/* phase = x / (2*PI) */
int sinscale(float phase, float scale)
{
    /* Optional: Handle negative phases */
    if (phase < 0.0f) {
        phase = -phase;
        scale = -scale;
    }

    /* Optional: Handle phases outside a full period */
    if (phase >= 1.0f) {
        phase -= (float)(int)phase;
    }

    /* Second half of each period is the first half negated */
    if (phase >= 0.5f) {
        phase -= 0.5f;
        scale = -scale;
    }

    /* Second quarter is first quarter mirrored: */
    if (phase > 0.25f) {
        phase = 0.50f - phase;
    }

    /* Convert phase to radians. */
    float  x = phase * 6.28318530717958647692528676655900576839433879875021164f;
    const float  xx = x*x;
    x *= scale;

    /* Power series expansion.  First term is x. */
    const float  term1 = x;

    /* term -x^3/3! */
    x *= xx;
    x /= 2*3;
    const float  term3 = x;

    /* +x^5/5! */
    x *= xx;
    x /= 4*5; /* 4*5 */
    const float  term5 = x;

    /* -x^7/7! */
    x *= xx;
    x /= 6*7;
    const float  term7 = x;

    /* +x^9/9! */
    x *= xx;
    x /= 8*9;

    /* Sum terms in increasing order of magnitudes, to minimize rounding errors. */
    x -= term7;
    x += term5;
    x -= term3;
    x += term1;

    /* Rounding. */
    return (x < 0.0f) ? (x - 0.5f) : (x + 0.5f);
}
If you are interested in why one would implement the power series expansion of sin() that way: we want to do the addition of the terms in order of increasing magnitude, to minimize domain cancellation errors.

For N=65536 and S=32000, i.e. sinscale((float)i/(float)N, S) ≃ (int)round(S*sin(2*PI*i/N)), 99% of the values are exact, and the rest only differ by ±1. Same for N=65536 and S=32767.

For N=1048576 and S=262144,  91% of the values are exact, and the rest only differ by ±1.

In general, for S ≤ 262144, I do believe the maximum difference is ±1.

This works well for non-power of two N, including primes.


If you are using C++, then consider expanding what emece67 suggested, into something like the following:

Code: [Select]

#define  HALFPI  1.57079632679489661923132169163975144209858469968755291
#define  PI      3.14159265358979323846264338327950288419716939937510582
#define  TWOPI   6.28318530717958647692528676655900576839433879875021164

template <class T, int N, T M>
class sine_table {
    public:
        T value[N];

        constexpr sine_table(): value() {
            for (int i = 0; i < N; i++) {
                double  x = TWOPI * i / N;
                double  s = (double)M;

                if (x >= PI) {
                    x -= PI;
                    s  = -s;
                }

                if (x >= HALFPI) {
                    x = PI - x;
                }

                const double  xx = x*x;

                /* Power series begins with the x term */
                x *= s;
                const double  t1 = x;

                /* -x**3/3! */
                x *= xx;
                x /= 2*3;
                const double  t3 = x;

                /* +x**5/5! */
                x *= xx;
                x /= 4*5;
                const double  t5 = x;

                /* -x**7/7! */
                x *= xx;
                x /= 6*7;
                const double  t7 = x;

                /* +x**9/9! */
                x *= xx;
                x /= 8*9;
                const double  t9 = x;

                /* -x**11/11! */
                x *= xx;
                x /= 10*11;
                const double  t11 = x;

                /* +x**13/13! */
                x *= xx;
                x /= 12*13;
                const double  t13 = x;

                /* -x**15/15! */
                x *= xx;
                x /= 14*15;
                const double  t15 = x;

                /* +x**17/17! */
                x *= xx;
                x /= 16*17;
                const double  t17 = x;

                /* -x**19/19! */
                x *= xx;
                x /= 18*19;
                const double  t19 = x;

                /* +x**21/21! */
                x *= xx;
                x /= 20*21;

                x -= t19;
                x += t17;
                x -= t15;
                x += t13;
                x -= t11;
                x += t9;
                x -= t7;
                x += t5;
                x -= t3;
                x += t1;

                value[i] = (x < 0.0) ? (x - 0.5) : (x + 0.5);
            }
        }
};

#undef  HALFPI
#undef  PI
#undef  TWOPI
gives about 50-bit approximation at compile time (±8 ULP in double precision).  For example,
    constexpr sine_table<int16_t, 1048573, 32767>  table1;
    constexpr sine_table<int16_t, 65536, 32000>  table2;
take annoyingly long to compile, but the generated tables have
    table1.value[i] == round(32767*sin(2*TWOPI*i/1048573));
    table2.value[j] == round(32000*sin(2*TWOPI*j/65536));
with exact values including correct rounding (towards nearest integer, halfway away from zero).


If you use precalculated tables, I recommend you define preprocessor macros for the size (N) and scale (S), and use e.g.

Code: [Select]

#if N-0 == 0
#error Sine table size (N) is zero
#elif S-0 == 0
#error Sine table scale (S) is zero

#elif N == 4096 && S == 32000
#include "sin-4096-32000.h"

#elif N == 16384 && S == 32000
#include "sin-16384-32000.h"

#elif N == 65536 && S == 32000
#include "sin-65536-32000.h"

#else
#error Unknown sine table size (N) and/or scale (S).
#endif
with whatever macro names (N and S) you prefer.  In the preprocessor tests, N-0 yields zero for both N=0 and when N is undefined, and generates an error if N is not an integral constant.

[brucehoult]

If you are interested in why one would implement the power series expansion of sin() that way: we want to do the addition of the terms in order of increasing magnitude, to minimize domain cancellation errors.

For N=65536 and S=32000, i.e. sinscale((float)i/(float)N, S) ≃ (int)round(S*sin(2*PI*i/N)), 99% of the values are exact, and the rest only differ by ±1. Same for N=65536 and S=32767.

My much simpler and smaller code in ...

https://www.eevblog.com/forum/programming/compact-sin(x)-which-gets-to-within-0-1-x-0-to-2xpi/msg3880130/#msg3880130

... gets the exact same 16 bit integer value as using the C library sin() function in 100% of cases (N=65536, S=32000 or 32767), despite not bothering with quadrants etc (it just uses the raw radian value up to 2*Pi in the power series) and not bothering to add up smallest terms first.

It seems using a sufficient number of terms is more important than the other factors.

My code is slower (calculates more terms), but the OP says that's not important.

[SiliconWizard]

Now I still think a small tool to generate this at compile time makes more sense if the N is fixed at compile time. But, this can still be a fun exercise to do this efficiently at run time. Or could be mandatory if N can  be dynamic.

Apart from the approaches already shown above, my "favorite" one for incremental generation of sine/cosine is to use a rotation matrix. Pretty efficient, *and* accurate, useful in particular for real-time generation!

For initializing the generation, you still need to compute the two constants: sin(pi/N) and sin(2pi/N). But if the number of different values of N is fixed and reasonably small, those can be computed in advance and put in a small constant array. Benefit of this algorithm is that you get both sin() and cos() at each iteration. You can discard the cos value, or use it.

Regarding the sin(pi/N) and sin(2pi/N) constants, if N is not too small, you can even approximate both by just pi/N and 2pi/N, not requiring any constant array. I computed tables with both approaches and for N >= 64, the error over the whole table is pretty acceptable!

Code: [Select]

#define kPI 3.14159265358979323846

static void MakeSinCosTable(double *pSinTable, double *pCosTable, unsigned int n)
{
if (n == 0)
return;

double x = kPI / (double) n, y = sin(x), Alpha = 2.0 * y * y, Beta = sin(2.0 * x),
Si = 0.0, Ci = 1.0;

for (unsigned int i = 0; i < n; i++)
{
if (pSinTable != NULL)
pSinTable[i] = Si;

if (pCosTable != NULL)
pCosTable[i] = Ci;

double Si1 = Si - (Alpha*Si - Beta*Ci);
double Ci1 = Ci - (Alpha*Ci + Beta*Si);

Si = Si1;
Ci = Ci1;
}
}

As said above, you can replace the two sin() values required for computing coefficients either by the approximation x and 2*x (try it and see what you get!), or using a table for various values of n, as long as n has a finite (and resonable) number of possible values.