/* 
 * sincos_table.c: 三角関数をテーブル引きにして高速化
 * 
 *     [2002/08/16] OSHIRO Naoki. <oshiro@mibai.tec.u-ryukyu.ac.jp>
 *   
 */

/* compile: gcc -O3 -c sincos_table.c */
/* 最適化オプション（-O3 など）を必ず付けること！    */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static int SinCos_num=-1;
static double *SinCos_table=NULL;

#include "sincos_table.h"

/* SinCos_init: ９０度分の cos() テーブルを作成する */
void SinCos_init(int num)
{
    int i;

    if (num<=0) {
	fprintf(stderr, "SinCos: specified init number error (%d).\n", num);
	exit(1);
    }
    SinCos_table=malloc(sizeof(double)*num);
    if (SinCos_table==NULL) {
	fprintf(stderr, "SinCos: not enough memory.\n");
	exit(1);
    }
    SinCos_num=num;
    for (i=0; i<num; i++) {
	SinCos_table[i]=cos(i*M_PI/2/num);
    }
}

/* SinCos テーブル範囲を毎回チェック．定義すると安全だが２０倍は遅くなる */
#undef SINCOS_CHECK

/* Cos() 関数 */
double Cos(double x)
{
    double sign=1;
    int i;

#ifdef SINCOS_CHECK
    if (SinCos_num<=0) {
	fprintf(stderr, "Not initialize SinCos_table.\n");
	exit(1);
    }
#endif
    if (x<0) x=-x;              /*  cos(x)=cos(-x)     */
    x/=M_PI;
    x=(x-(int)(x/2)*2);         /*  cos(x)=cos(x+2*PI) */
    if (x>1) x=2-x;             /*  cos(x)=cos(x-PI)   */
    if (x>.5) {x=1-x; sign=-1;} /* -cos(x)=cos(PI-x)   */
    i=(int)(x*2*SinCos_num);
#ifdef SINCOS_CHECK
    if (i<0 || i>=SinCos_num) {
	fprintf(stderr, "SinCos: specified number out of range (%d).\n", i);
	exit(1);
    }
#endif
    return (sign*SinCos_table[i]);
}

/************************************************************/
#ifdef SINCOS_TABLE_TEST
#include <time.h>

int main(void)
{
    int i, j, num;
    time_t t;
    double x;

    /* SinCosテーブル初期化 */
    /* 分割数と誤差の目安: 分割数 1e+X, 誤差 1e-(X+1) */
    num=1e+5;
    SinCos_init(num);
    printf("SinCos_Table:%d\n", num);
#if 1
    /* 0.01 [rad] 刻みで sin() と Sin(), cos() と Cos() の精度比較をする */
    printf("x,sin(),Sin(),sin()-Sin(x)\n");
    for (x=-10; x<10; x+=0.01) {
	printf("%f %f %f %+e\n", x, sin(x), Sin(x), sin(x)-Sin(x));
    }
    printf("\n");
    printf("x,cos(),Cos(),cos()-Cos(x)\n");
    for (x=-10; x<10; x+=0.01) {
	printf("%f %f %f %+e\n", x, cos(x), Cos(x), cos(x)-Cos(x));
    }
    printf("\n");
#endif

    /* cos() と Cos() の速度比較 */
    printf("Time benchmark: cos(), Cos()\n");
    /* cos() */
    t=time(NULL);
    for (i=0; i<50000; i++)
	for (j=0; j<10000; j++)
	    cos(.1);
    printf("cos:%dsec\n", (time(NULL)-t));

    /* Cos() */
    t=time(NULL);
    for (i=0; i<50000; i++)
	for (j=0; j<10000; j++)
	    Cos(.1);
    printf("Cos:%dsec\n", (time(NULL)-t));
}
#endif /* SINCOS_TABLE_TEST */
/************************************************************/