cordic2c.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <unistd.h>
#include <libgen.h>


// =================================================================
typedef int Cordic_T;
typedef unsigned int UCordic_T;
typedef double FCordic_T;

// Number of bits in Cordic_T
#define Cordic_T_Bits        (int)(sizeof(Cordic_T) << 3)
#define FCordic_T_Bits       (int)(sizeof(FCordic_T) << 3)
// 3 integer bits 1 sign bin
#define Cordic_T_FractionBits    (int)(Cordic_T_Bits - 3 - 1)

#define Cordic_T_STR "typedef int Cordic_T;"
#define UCordic_T_STR "typedef unsigned int UCordic_T;"
#define FCordic_T_STR "typedef double FCordic_T;"
// =================================================================

static Cordic_T v_atan[Cordic_T_Bits+1];
static Cordic_T X;
static Cordic_T Y;
static Cordic_T Z;

#define Cordic_One      (1UL << Cordic_T_FractionBits)
#define Cordic_K        (Cordic_One * 0.6072529350088812561694)
#define Cordic_INVK     (Cordic_One / 0.6072529350088812561694)
#define Cordic_KP       (Cordic_One * 1.20749706776307212887772)
#define Cordic_INVKP    (Cordic_One * 1/1.20749706776307212887772)
#define Cordic_HalfPI   ((UCordic_T) (Cordic_One * M_PI_2) )
#define Cordic2FP(a)    ( (double) (a) / (double) (Cordic_One)) 
#define FP2Cordic(a)    ((Cordic_T) (Cordic_One * (a)))

char *get_date()
{
    int len;
    char *ptr;
    time_t timev;
        timev = time(0);
        ptr = asctime(localtime(&timev));
    len = strlen(ptr);
    ptr[len-1] = 0;
    return(ptr);
}


void dump_tables(FILE *FO)
{
    int i;
    Cordic_T xx;

// Dump types
    fprintf(FO,"%s /* %u */\n", Cordic_T_STR, Cordic_T_Bits);
    fprintf(FO,"%s /* %u */\n", UCordic_T_STR, Cordic_T_Bits);
    fprintf(FO,"%s /* %u */\n", FCordic_T_STR, FCordic_T_Bits);
    fprintf(FO,"#define Cordic_T_Bits     %d\n", (int) Cordic_T_Bits);
    fprintf(FO,"#define Cordic_T_FractionBits  %d\n", (int) Cordic_T_FractionBits);


// Dump defines
    fprintf(FO,"#define Cordic_One    0x%lx /* %.15le */\n", 
        (long) Cordic_One, Cordic2FP(Cordic_One));
    fprintf(FO,"#define Cordic_K      0x%lx /* %.15le */\n", 
        (long) Cordic_K ,Cordic2FP(Cordic_K));
    fprintf(FO,"#define Cordic_INVK   0x%lx /* %.15le */\n", 
        (long) Cordic_INVK ,Cordic2FP(Cordic_INVK));
    fprintf(FO,"#define Cordic_KP     0x%lx /* %.15le */\n", 
        (long) Cordic_KP ,Cordic2FP(Cordic_KP));
    fprintf(FO,"#define Cordic_INVKP  0x%lx /* %.15le */\n", 
        (long) Cordic_INVKP ,Cordic2FP(Cordic_INVKP));
    fprintf(FO,"#define Cordic_HalfPI 0x%lx /* %.15le */\n", 
        (long) Cordic_HalfPI ,Cordic2FP(Cordic_HalfPI));
    fprintf(FO,"#define Cordic2FP(a)  ( (double) (a) / (double) (Cordic_One)) \n");
    fprintf(FO,"#define FP2Cordic(a)  ((Cordic_T) (Cordic_One * (a)))\n");

    fprintf(FO,"#ifdef CORDIC_TABLE\n");
// Dump normalize atan table
    xx = Cordic_One;
    for(i=0;i<Cordic_T_Bits;++i) {
            v_atan[i] = Cordic_One * atan( Cordic2FP(xx)) * (2.0/M_PI);
            xx >>= 1;
    }
    fprintf(FO,"static const Cordic_T v_atangrad[] = {\n");
    for(i=0;i<=Cordic_T_FractionBits;++i) {
        fprintf(FO,"\t0x%lx, /* %.8le */\n", 
        (long) v_atan[i], Cordic2FP(v_atan[i]));
    }
    fprintf(FO,"\t0\n};\n");
}

void PrintXYZ (char *str)
{
    printf("/* %s\n",str);
    printf("X: %.8lf\n", Cordic2FP(X));
    printf("Y: %.8lf\n", Cordic2FP(Y));
    printf("Z: %.8lf\n", Cordic2FP(Z));
    printf("*/\n");
}


void Circular (Cordic_T x, Cordic_T y, Cordic_T z)
{
    int i;

    X = x;
    Y = y;
    Z = z;

    for (i = 0; i <= Cordic_T_FractionBits; ++i)
    {
        x = X >> i;
        y = Y >> i;

        if(i < 14)
            z = v_atan[i];
        else
            z >>= 1;

        if(Z >= 0) {
            X -= y;
            Y += x;
            Z -= z;
        }
        else {
            X += y;
            Y -= x;
            Z += z;
        }
    }
}

void cordit1(Cordic_T x, Cordic_T y, Cordic_T z, Cordic_T vecmode)
{

    int i;

    X = x;
    Y = y;
    Z = z;

    for (i = 0; i < Cordic_T_FractionBits; ++i) {
        x = X >> i;
        y = Y >> i;
        if(i < 14)
            z = v_atan[i];
        else
            z >>= 1;
        if (vecmode >= 0 && Y < vecmode || vecmode < 0  && Z >= 0) {
            X -= y;
            Y += x;
            Z -= z;
        }
        else {
            X += y;
            Y -= x;
            Z += z;
        }
    }
}

Cordic_T asinCordic(Cordic_T a)
{
    Cordic_T x, y, z;
    

    x = Cordic_K;
    y = 0;
    z = 0;

    int neg = 1;
    if (a < 0) {
        a = -a;
        neg = 0;
    }
        

    cordit1(x, y, z, a);

    if (neg) Z = -Z;
    return Z;
}


int main (int argc, char *argv[])
{
    Cordic_T x1;
    char str[256];
    char *oname;
    char *p;
    double d;
    int a;
    int i;
    FILE *FO;

    for(i=1; i<argc;++i)
    {
        p = argv[i];
        if(*p != '-')
                continue;
        ++p;
        if(*p == 'o')
        {
                oname = argv[++i];
        }
    }
    if(!oname)
    {
        fprintf(stderr,"Usage: %s -o filename [-p]\n",argv[0]);
        fprintf(stderr,"-o filename is CORDIC C table output file\n");
        exit(1);
    }

    FO = fopen(oname,"w");
    if(FO == NULL)
    {
            fprintf(stderr,"Can not open: [%s]\n", oname);
            exit (1);
    }


    fprintf(FO,"#ifndef _CORDIC_INC_H\n");
    fprintf(FO,"#define _CORDIC_INC_H\n");

    fprintf(FO,"/**\n");
    fprintf(FO," @file %s\n", basename(oname));
    fprintf(FO," Generated by:[%s]\n", basename(argv[0]));
    fprintf(FO," On: %s\n", get_date());
    fprintf(FO," By Mike Gore 2015, Cordic C Table\n");
    fprintf(FO,"*/\n");


    dump_tables( FO);

    fprintf(FO,"#else // CORDIC_TABLE\n");
    fprintf(FO,"extern const Cordic_T v_atangrad[];\n");
    fprintf(FO,"#endif // CORDIC_TABLE\n");
    fprintf(FO,"#endif // _CORDIC_INC_H\n");

    printf("// Verify CORDIC table \n");
    for(d=0;d<=1;d+=.1)
    {
        x1 = FP2Cordic(d);
        Circular (Cordic_K, 0L, x1);
        a = (int) 100 * d + 0.000005;   // rounding - 0.1 is not exact
        sprintf(str,"%d Gradians", a);
        PrintXYZ (str);
    }
    printf("// End of CORDIC verify\n");

    return (0);
}