#include #include /* This is test/prototype code. Compile with gcc. The idea is * to work out the details of these algorithms and have a working * integer-only example to compare against when I try to write it * in assembly language */ /* This algorithm is described in "Computer Arithmetic Algorithms" * Ch 9, by Israel Koren, Prentice Hall, 1993, ISBN 0-13-151952-2 */ /* More info on the cordic algorithm * * http://www.andraka.com/cordic.htm * http://www.andraka.com/files/crdcsrvy.pdf * http://www.voidware.com/cordic.htm */ /* HELP: Info needed.... * * the paper by C. Mazenc, X. Merrheim, J.-M. Muller * IEEE transactions on computers, vol 42 no 1, pg 118-122 * is supposed to have the improvements to the cordic algorithm * for accurate arcsin and full-range arccos. * */ #define VERBOSE //#define PRINT_TABLE union f2l_union { float f; unsigned long l; }; #define NUM_STEPS 25 long table[NUM_STEPS]; double kval=1.0; void mk_table(void) { int i; double n; long sum=0; for (i=0; i> 8) & 255, (table[i] >> 16) & 255, n); #endif } printf("sum = %08lX\n", sum); } long do_trig(long x) { long w, z, w1, z1; unsigned int i, j; w = 0; z = (long)((1.0 / kval) * (double)0x1000000); for (i=0; i= 0) { x -= table[i]; z -= w1; w += z1; } else { x += table[i]; z += w1; w -= z1; } #ifdef VERBOSE printf("trig: x = %8ld, ", x); printf("sin = %08lX (%.7f), cos = %08lX (%.7f)\n", w, (double)w / (double)0x1000000, z, (double)z / (double)0x1000000); #endif } return w; } long arc_tan(long x) { int i, j; long u, v, y, u1, v1; u = x; v = 0x1000000; y = 0; for (i=0; i= 0) { y += table[i]; u -= v1; v += u1; } else { y -= table[i]; u += v1; v -= u1; } #ifdef VERBOSE printf("trig: atan = %8ld (%.8f), ", y, (double)y / (double)0x1000000); printf("u = %08lX, ", u); printf("v = %08lX\n", v); #endif } return u; } // working, but accuracy very poor for input greater than 0.9 // the paper by C. Mazenc, X. Merrheim, J.-M. Muller // IEEE transactions on computers, vol 42 no 1, pg 118-122 // is supposed to have the answer. long arc_sin(long c) { int i, j, d; long x, y, z, x1, y1; x = (long)((1.0 / kval) * (double)0x1000000); y = 0; z = 0; for (i=0; i