Code used for sine table

[orolo]

Orolo is manually drawing a sine and cosine just like the first basic program post by GK -> https://www.eevblog.com/forum/microcontrollers/code-used-for-sine-table/msg1107610/#msg1107610

The differential equation is a great idea, but I wasn't doing that. To put it briefly, it took the complex number z = cos(1) + I*sin(1), and computed its powers, 1, z, z^2, z^3,... each time multiplying, in the complex plane, by z. Since z^n = cos(n) + I*sin(n), you get the sine (and cosine) table.

In binary arithmetic, sin(1) = 0,0000010001111000... and cos(1)=0,1111111111110110... . I used the fact that 1 = 0,11111111... in binary, also.

The final ' t = (sa >> 2) + ((sa&2)>>1); ' just scales down the results so that the table goes from 0-16384 instead of 0-65535.

Exactly. The ((sa&2)>>1) part adds rounding to the scaling down. Just using t = sa>>2 is more compact, but loses a bit of precision.

Thank you for the interest and the graphical analysis, that was great. 

Edit:

Just to keep  , if a little increase in complexity is allowed, the table can be computed with more precision and half the number of iterations. The idea is using the cosines to compute half the table of sines:

Code: [Select]

#include <stdio.h>
int main()
{
   unsigned short sx, cx, sa, ca, t, t2;
   char c;

   sx = 0b0000010001111000;
   cx = 0b1111111111110110;
   sa = 0;
   ca = 0b1111111111111111;
   
   for (c=0 ; c < 46 ; c++) {
    t = (sa >> 2) + ((sa&2)>>1);
    t2 = (ca >> 2) + ((ca&2)>>1);
    printf("%d %d\n", t, t2);
   
    t  = ((long) (sa*cx + ca*sx))>>16;
    ca = ((long) (ca*cx - sa*sx))>>16;
    sa = t;
   }
}
The idea is to fill the lower half of the table with t, and the upper part with t2, in descending order. With this trick, the precision of the upper part of the table increases remarkably. Both parts of the table meet in the middle with the same value. And this time the table starts at 0! (thanks, Brian).

[BrianHG]

Here you go Orolo,  Green is the original sin(x) command.  The Cyan is your new code with the (long) addition/subtraction between the 2 multiplies.  Red is your new code adding/subtracting only the upper 16 bits from the 2 multiplies as this coding for 16 bit CPU with a hardware multiplier with the upper 16 bit results would be the most compact means generating this sine table.

Code: [Select]

' http://www.freebasic.net/
' for dos, win & linux
' sintable comparison generator by BrianHG.

const pi = 3.141592654  : Rem This is the same thing >>> atn(1) * 4

const xscale  = 10
const xoffset = 50
const yscale  = 30
const yoffset = 800

Declare Sub plot_sin_table()

dim Shared as unsigned integer x,y,n,cr,cg,cb,tmp,sx,cx,sa,ca,t,c,t2
dim Shared sin_table(0 to 90) as integer

ScreenRes 1152,816,32,0,0


' Generate sine_table using the function sin(x).
   for n=0 to 90
    sin_table(n) = sin(pi*(n/90)/2) * 16384
   next n

' Plot sin_table(n) in green
cr=0:cg=255:cb=0:plot_sin_table()


' approximate sine_table 'orolo's NEW code using Long.
' please excuse basic's limitation of shifting bits, I had to replace them with divide...
' It is expected that in the CPU, you would use the bit shifting for speed and code size.

   sx = 1144  ' 0b0000010001111000;
   cx = 65526 ' 0b1111111111110110;
   sa = 0     ' 0;
   ca = 65535 ' 0b1111111111111111;
   
   for c=0 to 46

    t  = int(sa / 4) + ((sa and 2)/2)
    t2 = int(ca / 4) + ((ca and 2)/2)

    sin_table( c    ) = t
    sin_table( 90-c ) = t2

    t  = ((sa*cx) + (ca*sx))/65536
    ca = ((ca*cx) - (sa*sx))/65536
    sa = t

   next c

' Plot sin_table(n) in Cyan
cr=0:cg=255:cb=255:plot_sin_table()


' approximate sine_table 'orolo's NEW code using Integer multiply.
' please excuse basic's limitation of shifting bits, I had to replace them with divide...
' It is expected that in the CPU, you would use the bit shifting for speed and code size.

   sx = 1144  ' 0b0000010001111000;
   cx = 65526 ' 0b1111111111110110;
   sa = 0     ' 0;
   ca = 65535 ' 0b1111111111111111;
   
   for c=0 to 46

    t  = int(sa / 4) + ((sa and 2)/2)
    t2 = int(ca / 4) + ((ca and 2)/2)

    sin_table( c    ) = t
    sin_table( 90-c ) = t2

    t  = int((sa*cx)/65536) + int((ca*sx)/65536)
    ca = int((ca*cx)/65536) - int((sa*sx)/65536)
    sa = t

   next c

' Plot sin_table(n) in Red
cr=255:cg=0:cb=0:plot_sin_table()



' Wait until a key if pressed before ending.
  while (inkey$="") : sleep 1 : wend
END

Sub plot_sin_table()
color RGB(cr,cg,cb)
  for n=0 to 90
   ? sin_table(n),
  next n
   ?:?

  for n=0 to 89
   line (n*xscale+xoffset, yoffset-(sin_table(n)/yscale)) - ((n+1)*xscale+xoffset, yoffset-(sin_table(n+1)/yscale) ), RGB(cr,cg,cb)
  next n
End Sub

As you can see in the resulting image, you will now need to download and zoom into the image, the cheap integer version, now red, has drastically improved over my previous integer version which had green segments showing through the cyan throughout the upper half of the sine wave.

[BrianHG]

If you mean this:

Code: [Select]

n = 0;
while (n < 91) {
sin_table[n] = sin(pi*(n/90)/2) * 16384;
n++;
}


I believe that the C compiler should simplify this one out, but, just in case it doesn't:
--------------------------------------------------------------
sin_table[n] = sin(pi*(n/90)/2) * 16384;
===is equivalent to===
sin_table[n] = sin( n*0.01745329252 ) * 16384;
--------------------------------------------------------------

[jesuscf]

One possible solution is to store the second differences.   After some quick spread sheet calculations, all the second differences fit in only three bits.  If they are stored in 32-bit constants, 10 3-bit differences can be stored per 32-bit constant.  That reduces the size of the original table from 182 bytes to just 36 bytes:

Code: [Select]

unsigned long packed_differences[] ={
   0x12451288, 0x142da492, 0x1551c4db,
   0x2c764724, 0x26965af3, 0x2dd6db6d,
   0x2ebb5bb5, 0x3dbf5db6, 0x05f7df76
}
Here is the program I used to test the whole idea.  It both generates the compressed table of second differences and un-compresses it to produce the original table.

Code: [Select]

#define WORD short int
#define TSIZE 91

WORD sin_table[TSIZE] ={
    0,   285,   571,   857,  1142,  1427,  1712,  1996,  2280,  2563,
    2845,  3126,  3406,  3685,  3963,  4240,  4516,  4790,  5062,  5334,
    5603,  5871,  6137,  6401,  6663,  6924,  7182,  7438,  7691,  7943,
    8191,  8438,  8682,  8923,  9161,  9397,  9630,  9860, 10086, 10310,
    10531, 10748, 10963, 11173, 11381, 11585, 11785, 11982, 12175, 12365,
    12550, 12732, 12910, 13084, 13254, 13420, 13582, 13740, 13894, 14043,
    14188, 14329, 14466, 14598, 14725, 14848, 14967, 15081, 15190, 15295,
    15395, 15491, 15582, 15668, 15749, 15825, 15897, 15964, 16025, 16082,
    16135, 16182, 16224, 16261, 16294, 16321, 16344, 16361, 16374, 16381,
    16384
};

WORD diff1[91];
WORD diff2[91];

long pack_diff2[10];
int y=0, d1=285, d2;

void main (void)
{
int j, k;

// Compute first differences
for(j=0; j<(TSIZE-1); j++)
{
diff1[j]=sin_table[j+1]-sin_table[j];
}

// Compute second differences
for(j=0; j<(TSIZE-2); j++)
{
diff2[j]=diff1[j]-diff1[j+1]+1; // Change the order and add one so differences are all positive
}

// Pack the differences: each second difference uses 3 bits.
printf("unsigned long packed_differences[] ={\n");
for(j=0; j<9; j++)
{
pack_diff2[j]=0;
for(k=0; k<10; k++)
{
pack_diff2[j]<<=3;
pack_diff2[j]|=(diff2[(j*10)+(9-k)]&0x07);
}
printf("   0x%08x,\n", pack_diff2[j]);
}
printf("};\n");

// With the packed differences, rebuilt the original table:
printf("Un-Packed table:\n");
for(j=0; j<9; j++)
{
for(k=0; k<10; k++)
{
printf("%5d ", y);
y+=d1;
d2=pack_diff2[j]&0x7;
pack_diff2[j]>>=3;
d1=d1-d2+1;
}
printf("\n");
}
printf("%d\n", y);
}
This is the rebuilt sine table from the packed table of second differences:

Code: [Select]

    0   285   571   857  1142  1427  1712  1996  2280  2563
 2845  3126  3406  3685  3963  4240  4516  4790  5062  5334
 5603  5871  6137  6401  6663  6924  7182  7438  7691  7943
 8191  8438  8682  8923  9161  9397  9630  9860 10086 10310
10531 10748 10963 11173 11381 11585 11785 11982 12175 12365
12550 12732 12910 13084 13254 13420 13582 13740 13894 14043
14188 14329 14466 14598 14725 14848 14967 15081 15190 15295
15395 15491 15582 15668 15749 15825 15897 15964 16025 16082
16135 16182 16224 16261 16294 16321 16344 16361 16374 16381
16384
This is a complete test program using the compressed table of second differences.  It generates the original table and prints it.  I am sure it can be optimized a bit more...

Code: [Select]

#define WORD short int
#define TSIZE 91

unsigned long pack_diff2[9] = {
   0x12451288, 0x142da492, 0x1551c4db,
   0x2c764724, 0x26965af3, 0x2dd6db6d,
   0x2ebb5bb5, 0x3dbf5db6, 0x05f7df76
};

WORD sin_table[TSIZE];

void main (void)
{
WORD y=0, d1=285;
unsigned long d2;
char i=0, j, k;

// Unpack the table
for(j=0; j<9; j++)
{
d2=pack_diff2[j];
for(k=0; k<10; k++)
{
sin_table[i++]=y;
y+=d1;
d1=d1-(d2&0x7)+1;
d2>>=3;
}
}
sin_table[i]=y;

// Just to make sure, print the table
for(i=0; i<TSIZE; i++)
{
printf("%5d ", sin_table[i]);
if((i%10)==0) printf("\n");
}
printf("\n");
}
Jesus

[jesuscf]

Solving these two coupled differential equations will do it.

dsine/dt = cos
dcos/dt = -sine

I'm not a mathematician either and I'm not necessarily saying that this is an efficient solution, and the code example is in BASIC rather than C. This is how I solved the problem on one lazy day just out of idle curiosity on my Beeb. I can't figure an integration algorithm stripped any barer.

The way the solution is implemented, using Forward Euler integration, produces results that increase in magnitude over time (FE.jpg).  This is because Forward Euler integration introduces "negative" damping.  One can use the Backward Euler integration, but then the opposite problem arises as Backward Euler introduces positive damping and produces results that decrease in magnitude over time (BE.jpg) .  The solution to this problem is to use trapezoidal integration which does not introduce any damping (trap.jpg).  The resulting equations using trapezoidal integration are

f1(x)=2*h*f2(x-h)+f1(x-2h)
f2(x)=-2*h*f1(x-h)+f2(x-2h)

where h is the time step, f1(x) is the sine at time x, f1(x-h) is the sine at the previous time step, and f1(x-2h) is the sine at the time step before that.  Similarly f2(x) is the cosine at time x, f2(x-h) is the cosine at the previous time step, and f2(x-2h) is the cosine before that.  table_trap.jpg shows the numeric results obtained with the equations above and a time step of one degree (pi()/180=0.01745).

[CatalinaWOW]

Lots of good answers here.

For this and all other problems relating to digital implementations of math, one of my first stops is the Handbook of Mathematical Functions by Abramowicz and Stegun.  It has formal definitions, series approximations, rational approximations and tabular data for most functions anyone will encounter.  It is available online from several sources, one of which is here:

http://people.math.sfu.ca/~cbm/aands/

[@rt]

Another thing which may or may not be significant is the nature of the data... your table is constant.  It therefore needs to end up in Flash not get allocated and initialised in RAM.  If you are lucky adding const will allow the compiler to work out it is immutable.

Something I thought about more. Is this a formal Computer Science assertion?
Thinking about it more, if you can’t trust the content of RAM, you can’t trust any software.
You can’t even trust the program counter, so what’s the difference?

[jesuscf]

I have some results.  Using a table of packed second differences to generate the sine table resulted in some code size reduction.  Compiled for an 8051, I wrote the unpacking routine in assembly and changed the table of second differences from long to int.  It takes 108 bytes of code memory to initialize the table in expanded RAM.  Here is the program:

Code: [Select]

#define TSIZE 91

const int pack_diff2[18] = {
0x1288, 0x248a, 0x2492, 0x285b, 0x44db, 0x2aa3,
0x4724, 0x58ec, 0x5af3, 0x4d2c, 0x5b6d, 0x5bad,
0x5bb5, 0x5d76, 0x5db6, 0x7b7e, 0x5f76, 0x0bef
};

xdata unsigned int sin_table[TSIZE];

void expand_table (void) _naked
{
_asm
mov r4, #18 ; loop counter 1
mov dptr, #_pack_diff2
MOV P2, #(_sin_table/0x100) ;  Use _EMI0CN in C8051F38x instead
MOV R0,#(_sin_table%0x100)
; First entry in table is zero thanks to crt0 (skip setting it to zero)
;clr a
;movx @R0, a
;inc R0
;movx @R0, a
;dec R0
; Initialize first D1
mov r2, #(285%0x100)
mov r3, #(285/0x100)
L1:
; Get packed 3-bit D2 x 5 from compressed table.  Store in [r7,r6]
clr a
movc a, @a+dptr
mov r6, a
inc dptr
clr a
movc a, @a+dptr
mov r7, a
inc dptr

mov r5, #5 ;  loop counter 2
L2:
; add d1 to previous value and save in xram table
mov a, r0
add a,#2
mov r1, a

movx a, @R0
add a, r2
movx @r1, a
inc R0
inc r1
movx a, @R0
addc a, r3
movx @r1, a
inc r0

; d1=d1+1-(d2&0x7);
mov a, r6
anl a, #7
mov r1, a

clr c
mov a, r2
subb a, r1
mov r2, a
mov a, r3
subb a, #0
mov r3, a

inc r2
cjne r2, #0, L3
inc r3
L3:
; Shift compressed D2 x 5 3-bits right:
mov r1, #3
L4:
mov a, r7
rrc a
mov r7, a
mov a, r6
rrc a
mov r6, a
djnz r1, L4

djnz r5, L2
djnz r4, L1
ret
_endasm;
}

void main (void)
{
expand_table();
}

The program with the whole table uses 183 bytes of code memory (there is an extra 'ret' instruction that takes one byte added to the main program):

Code: [Select]

#define WORD short int
#define TSIZE 91

const WORD sin_table[TSIZE] ={
    0,   285,   571,   857,  1142,  1427,  1712,  1996,  2280,  2563,
    2845,  3126,  3406,  3685,  3963,  4240,  4516,  4790,  5062,  5334,
    5603,  5871,  6137,  6401,  6663,  6924,  7182,  7438,  7691,  7943,
    8191,  8438,  8682,  8923,  9161,  9397,  9630,  9860, 10086, 10310,
    10531, 10748, 10963, 11173, 11381, 11585, 11785, 11982, 12175, 12365,
    12550, 12732, 12910, 13084, 13254, 13420, 13582, 13740, 13894, 14043,
    14188, 14329, 14466, 14598, 14725, 14848, 14967, 15081, 15190, 15295,
    15395, 15491, 15582, 15668, 15749, 15825, 15897, 15964, 16025, 16082,
    16135, 16182, 16224, 16261, 16294, 16321, 16344, 16361, 16374, 16381,
    16384
};

void main (void)
{
}

Using a table of second differences to generate the sine table saves 75 bytes of code memory in the 8051.

Jesus

[jesuscf]

Attached is the Excel spread sheet I used to test some of the ideas presented in this thread.