problem with lookup sine table

[Pallav]

i'm trying to build a lookup sine table to run some code in an
embedded processor without floating point unit. the table is 256
values for 360 degrees.

my code is below. the fixed point base is 22.10.

what i'm seeing is that the code seems to work for small integers say
(-10 - 10), assuming the value is in radians. however as you go up, to
large integers there is significant mismatch in the LUT and what the
sin() function returns.

i'm not sure what the problem is. i thought it could have been
overflow or i'm not performing the right conversion during lookup.
however, i seem to be convinced that i'm not running into either
problem and i'm somewhat stuck.

can anyone point out what is going wrong? the full code is below, the
code to generate the sin table is alwso below. it should compile and
the test program will exit if the difference is greater than 10%.

thanks
pallav

/* fixed.h*/

/*
* Fixed Point Arithmetic for 18.14 format
*/

#ifndef FIXED_H
#define FIXED_H

typedef long fp14;

#define DEGREES 0xFF
/* THIS WILL GIVE ABOUT 3 DECIMAL DIGITS OF PRECISION */
#define FRACTION 10

#define FP_EXP FLOAT2FP(2.71828)
#define FP_ONE INT2FP(1)
#define MY_PI 3.141592654
#define FP_PI FLOAT2FP(MY_PI)
#define FP_MAGIC FLOAT2FP(40.72) // 256/2*pi


#define FPMAX(a,b) (((a) > (b)) ? (a) : (b))
#define FPMIN(a,b) (((a) < (b)) ? (a) : (b))
#define FPABS(a,b) (((a) < 0) ? -(a) : (a))
#define FPSIGN(a) (((a) < 0) ? -1 : (a) > 0 ? 1 : 0)


/* CONVERSTION UTILITIES */
#define INT2FP(a) ((a) << FRACTION)
#define FP2INT(a) ((a) >> FRACTION)
#define FPFRACT(a) ((a) & 0x3FF)
#define FP2FLOAT(a) (((double) (a)) / (1 << FRACTION))
#define FLOAT2FP(a) ((fp14) (a * (1 << FRACTION)))


/* ARITHMETIC FUNCTIONS */
#define FPADD(a1,a2) ((a1) + (a2))
#define FPSUB(a1,a2) ((a1) - (a2))
#define FPMUL(a1,a2) ((fp14) ((((long long int) a1) * ((long long int)
a2)) >> FRACTION))
#define FPDIV(a1,a2) (FLOAT2FP(FP2FLOAT(a1)/FP2FLOAT(a2)))

/* SQUARE ROOT OF X^2 + Y^2 */
#define FPHYPOT(a1,a2) (FPSQRT(FPADD(FPMUL(a1,a1), FPMUL(a2,a2))))

/* SQUARE ROOT OF X^2 - Y^2 */
#define FPHYPOTDIFF(a1,a2) (FPSQRT(FPSUB(FPMUL(a1,a1), FPMUL(a2,a2))))


/* TRIGNOMETRIC FUNCTIONS */
#define FPSIN(a) (sine[(a) & DEGREES])
#define FPCOS(a) (cosine[(a) & DEGREES])
#define FPATAN(a) (invtan[(a) & DEGREES])

extern long sine[];
//extern long cosine[];
//extern long invtan[];

extern fp14 FPSQRT(fp14 a1);
extern fp14 FPEXP(fp14 a1);

#endif /* FIXED_H */

---------------------------------
/* fixed.c */

#include "fixed.h"

long sine[] =
{
0, 25, 50, 75, 100, 125, 150, 175,
200, 225, 249, 274, 298, 322, 346, 369,
393, 416, 439, 462, 484, 506, 528, 549,
570, 591, 612, 632, 651, 671, 689, 708,
726, 743, 760, 777, 793, 809, 824, 839,
853, 867, 880, 893, 905, 916, 927, 938,
947, 957, 965, 973, 981, 988, 994, 1000,
1005, 1009, 1013, 1016, 1019, 1021, 1023, 1023,
1023, 1023, 1022, 1020, 1018, 1015, 1011, 1007,
1002, 997, 991, 984, 977, 969, 961, 952,
943, 932, 922, 910, 899, 886, 873, 860,
846, 832, 817, 801, 785, 769, 752, 735,
717, 699, 680, 661, 641, 622, 601, 581,
560, 539, 517, 495, 473, 450, 427, 404,
381, 358, 334, 310, 286, 261, 237, 212,
188, 163, 138, 113, 88, 63, 37, 12,
-12, -37, -63, -88, -113, -138, -163, -188,
-212, -237, -261, -286, -310, -334, -358, -381,
-404, -427, -450, -473, -495, -517, -539, -560,
-581, -601, -622, -641, -661, -680, -699, -717,
-735, -752, -769, -785, -801, -817, -832, -846,
-860, -873, -886, -899, -910, -922, -932, -943,
-952, -961, -969, -977, -984, -991, -997, -1002,
-1007, -1011, -1015, -1018, -1020, -1022, -1023, -1023,
-1023, -1023, -1021, -1019, -1016, -1013, -1009, -1005,
-1000, -994, -988, -981, -973, -965, -957, -947,
-938, -927, -916, -905, -893, -880, -867, -853,
-839, -824, -809, -793, -777, -760, -743, -726,
-708, -689, -671, -651, -632, -612, -591, -570,
-549, -528, -506, -484, -462, -439, -416, -393,
-369, -346, -322, -298, -274, -249, -225, -200,
-175, -150, -125, -100, -75, -50, -25, };


/* SQUARE ROOT OPERATION: NOTE: a1 must be <= 0; no checks done */
fp14 FPSQRT(fp14 a1)
{
long root = 0, val1, val2, rem;
int i;

rem = a1;
for (i = 0; i <= 30; i += 2) {
val1 = (0x40000000L >> i);
val2 = val1 + root;
root >>= 1;
if (val2 <= rem) {
rem -= val2;
root |= (val1);
}
}

/* cut the error in half by rounding */
if (root < rem)
root++;

/* need to shift left 5 bits because FRACTION = 10 */
root <<= 5;

return ((fp14) root);
}




/* E^(x) - GIVES GOOD RESULTS ONLY FOR X in -3 to 3. SHOULD BE GOOD
ENOUGH
FOR FINGERPRINT COMPUTATIONS */
fp14 FPEXP(fp14 a)
{
int i;
fp14 prod, fact, approx, count;

count = FP_ONE;
fact = FP_ONE;
approx = FPADD(a, fact);
prod = a;

/* Use Taylor Series. Note only 11 terms are used for the
approximation
because 10! exceeds the amount of space we have to store the
integer */
for (i = 2; i <= 9; i++) {
count = FPADD(count, FP_ONE);
fact = FPMUL(fact, count);
prod = FPMUL(prod, a);
approx = FPADD(approx, FPDIV(prod,fact));
}

return approx;
}
---------------------------

/* test program */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "fixed.h"

int main(int argc, char **argv)
{
int u;
fp14 fpnx;
double lut, orig;
printf("%.6f\n", FP2FLOAT(FP_MAGIC));
for (u = 0; ; u += 10) {
fpnx = INT2FP(u);
fpnx = FPMUL(fpnx, FP_MAGIC);
lut = FP2FLOAT(FPSIN(FP2INT(fpnx)));
orig = sin(u);

if (fabs(lut - orig) > 0.10) {
printf("at iteration u= %d values differ more than 10%. not
good\n", u);
exit (0);
}

printf("%d %f %f %f\n", u, FP2FLOAT(fpnx), lut, orig);
}

return 0;
}
----------------
/* code to generate sin table */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "fixed.h"

int main(void)
{
int i, j;
double angle;
double sine;
long isine;
double pi = 3.141592654;
double delta = (double) DEGREES / (2 * pi);

printf("long sine[] =\n");
printf("{\n");
for (i = 0, j = 1; i < DEGREES; i++, j++)
{
angle = (double) (i/delta);
sine = sin(angle);
printf("%6d, ", FLOAT2FP(sine));
if (j > 7)
{
j = 0;
printf("\n");
}
}
printf("};\n");

return 0;
}