/* *******************************************************************
  app_math_math.c

  ====================================================================
  void app_math_init_table()
  --------------------------------------------------------------------
  機能：  	下記の配列に，三角関数のテーブルを作成
				配列名・サイズ：	Gf_math_sin_table[N_SIN_TABLE + 1]
				有効数字は，N_SIN_TABLE が 512 の場合で５桁，2048 の場合で 7 桁．
				その他の数学関数使用開始前に，必ず一度実行すること．
  引数：	無し
  戻り値：	無し
  ====================================================================
  double app_math_sin( double f_rad )
  double app_math_cos( double f_rad )
  --------------------------------------------------------------------
  機能：  	位相（radian）を与えると，三角関数の sin および cos の値を返す．
  引数：	double	f_rad		位相（radian）
  戻り値：	sin および cos の値
  ====================================================================
  double app_math_atan( double f_tangent )
  --------------------------------------------------------------------
  機能：	Tangent 値を与えると，偏角（radian）を返す．
  			偏角は，0 ≦ θ ＜ 2π の範囲で返る．
  引数：	double	f_rad		Tangent 値
  戻り値：	偏角（radian，0 ≦ θ ＜ 2π）
  ====================================================================
  double app_math_log( double f_org )
  --------------------------------------------------------------------
  機能：	値を与えると，その自然対数を返す．
  			値が０以下の場合，可能な限り負に大きい数を返す（本来，-Inf.）．
  引数：	double	f_org		自然対数を求める値（ f_org > 0 ）
  戻り値：	自然対数
  ====================================================================
  double app_math_log10( double f_org )
  --------------------------------------------------------------------
  機能：	値を与えると，その常用対数を返す．
  			値が０以下の場合，可能な限り負に大きい数を返す（本来，-Inf.）．
  引数：	double	f_org		常用対数を求める値（ f_org > 0 ）
  戻り値：	常用対数
  ====================================================================
  double app_math_exp( double f_org )
  --------------------------------------------------------------------
  機能：	指数を与えると，その e を底とするべき乗を返す．
  引数：	double	f_org		e を底とするべき乗を求める指数
  戻り値：	e のべき乗
  ====================================================================
  double app_math_pow( double f_base, double f_exp )
  --------------------------------------------------------------------
  機能：	底と指数を与えると，べき乗を返す．
  引数：	double	f_base		べき乗を求める底
  			double	f_exp		べき乗を求める指数
  戻り値：	べき乗
  ====================================================================
  double app_math_sqrt( double f_org )
  --------------------------------------------------------------------
  機能：	値を与えると，その平方根を返す．
  引数：	double	f_org		平方根を求める値（ f_org > 0 ）
  戻り値：	平方根
  ====================================================================
  double app_math_rad2deg( double f_rad )
  --------------------------------------------------------------------
  機能：	Radian を与えると，degree を返す．
			戻り値は 0 ≦ θ ＜ 360 の範囲で返る．
  引数：	double	f_rad		位相（radian）
  戻り値：	位相（degree, 0 ≦ θ ＜ 360）
  ====================================================================
  double app_math_regulate_rad( double f_rad )
  --------------------------------------------------------------------
  機能：	位相（radian）を 0 ≦ θ ＜ 2π の範囲に正規化する．
  引数：	double	f_rad		位相（radian）
  戻り値：	位相（radian, 0 ≦ θ ＜ 2π）
  ====================================================================
 ******************************************************************* */
#include		"app_common.h"

#define 		MATH_PI        (3.1415926536)
#define 		MATH_LN10      (2.302585093)
#define 		MATH_N_NEWTON  10

/*
double app_math_sin( double f_rad );
double app_math_cos( double f_rad );
double app_math_atan( double f_tangent );
double app_math_log( double f_org );
double app_math_log10( double f_org );
double app_math_exp( double f_org );
double app_math_pow( double f_base, double f_exp );
double app_math_sqrt( double f_org );
double app_math_rad2deg( double f_rad );
double app_math_regulate_rad( double f_rad );
*/


// sin テーブルを生成
void app_math_init_table(){
	double	f_rad, f_val, f_tmp;
	double	f_i;
	int	i;
	
	// sin の値が入った配列を生成．
	// 0 ≦ θ < π/2 の範囲を N_SIN_TABLE で分割．
	for( i=0; i<Gd_N_SIN_TABLE; i++ ){
		f_rad = ( MATH_PI / 2 ) / Gd_N_SIN_TABLE * i;
		
		// Taylor 級数により sin を計算
		f_val = f_rad;
		f_tmp = f_rad;
		for ( f_i = 1; f_i < 9; f_i++ ) {
			f_tmp = - f_tmp * f_rad * f_rad / ( ( 2 * f_i ) * ( 2 * f_i + 1 ) );
			f_val += f_tmp;
		}
		
		// 配列に収納
		Gf_math_sin_table[i] = f_val;
	}
	Gf_math_sin_table[Gd_N_SIN_TABLE] = 1.0;
}


// sin テーブルを読み，直線補完する．
double app_math_read_table( double f_rad ){
	int		i_num1, i_num2;
	double	f_val, f_val1, f_val2;
	double	f_delta, f_diff;
	
	// テーブルの配列番号を計算
	i_num1 = (int)( Gd_N_SIN_TABLE * f_rad / ( MATH_PI / 2 ) );
	i_num2 = i_num1 + 1;
	
	// テーブルの読み出し
	f_val1 = Gf_math_sin_table[i_num1];
	f_val2 = Gf_math_sin_table[i_num2];
	
	// 直線補完
	f_delta = ( MATH_PI / 2 ) / Gd_N_SIN_TABLE;
	f_diff = ( f_rad - ( f_delta * i_num1 ) ) / f_delta;
	f_val = f_val1 + ( f_val2 - f_val1 ) * f_diff;
	
	return f_val;
}


double 	app_math_sin( double f_rad )
{
	double 	f_val, f_sign;
// 	double 	f_i;							// **** REVISED:110618 by YK ****
	
	// 位相を 0 <= θ < 2π の範囲にブランチカット
	f_rad = app_math_regulate_rad( f_rad );
	
	// π <= θ < 2π :	θ - π で計算し，-1 倍
	if( f_rad < MATH_PI ){
		f_sign = 1;
	}
	else{
		f_sign = -1;
		f_rad -= MATH_PI;
	}
	
	// 0 <= θ < π/2 :	そのまま
	if( f_rad < MATH_PI / 2 ){
		f_val = app_math_read_table( f_rad ) * f_sign;
	}
	// π/2 <= θ < π :	π - θ で計算
	else{
		f_val = app_math_read_table( MATH_PI - f_rad ) * f_sign;
	}
	
  	return f_val;
}


double 	app_math_cos( double f_rad )
{
  	return app_math_sin( f_rad + MATH_PI / 2 );
}


double 	app_math_atan_base( double f_tangent )
{
  	double	f_tmp1, f_tmp2, f_tmp3, f_val;
  	int	i;

  	f_val = f_tangent;
  	f_tmp1 = f_tangent;
  	f_tmp3 = 3;
  	for ( i=0; i<10; i++ ) {
    	f_tmp1 = - f_tmp1 * f_tangent * f_tangent;
    	f_tmp2 = f_tmp1 / f_tmp3;
    	f_val += f_tmp2;
    	f_tmp3 += 2;
  	}

  return f_val;
}


double 	app_math_atan( double f_tangent ) {
  	double	f_param, f_val;

  	// 負の場合は，いったん符号を反転．
  	if ( f_tangent < 0 ){
    	f_param = - f_tangent;
  	}
  	else {
    	f_param = f_tangent;
  	}

  	// 大きさが 0.5 未満の場合，そのまま Taylor 級数で計算．
  	if ( f_param < 0.5 ) {
    	f_val = app_math_atan_base( f_param );
  	}
  	// 大きさが 1 未満の場合，加法定理に基づき，結果が 1/2 になるように
  	// 計算して２倍する．
  	else if ( f_param < 1 ) {
    	f_param = ( -1 + app_math_sqrt( 1 + f_param * f_param ) ) / f_param;
    	f_val = app_math_atan_base( f_param );
    	f_val *= 2;
  	}
  
  	// 大きさが 2 未満の場合，逆数について，加法定理に基づき，結果が 1/2 に
  	// なるように計算して２倍し，π/2 から差し引く．
  	else if ( f_param < 2 ) {
    	f_param = 1 / f_param;
    	f_param = ( -1 + app_math_sqrt( 1 + f_param * f_param ) ) / f_param;
    	f_val = app_math_atan_base( f_param );
    	f_val *= 2;
    	f_val = MATH_PI / 2 - f_val;
  	}
  	// 大きさが 2 以上の場合，逆数について計算し，π/2 から差し引く．
  	else {
    	f_param = 1 / f_param;
    	f_val = app_math_atan_base( f_param );
    	f_val = MATH_PI / 2 - f_val;
  	}
  	
  	// 負の場合は，結果の符号を再反転．
  	if( f_tangent < 0 ){
    	f_val *= -1;
  	}

  	// 結果を 0 < theta < 2π に正規化する．
  	f_val = app_math_regulate_rad( f_val );
  
  	return f_val;
}


double 	log_base( double f_org )
{
  	double	f_x;
  	double	f_tmp1, f_tmp2, f_tmp3, f_val;
  	int		i;

  	f_x = f_org - 1;

  	f_val = f_x;
  	f_tmp1 = f_x;
  	f_tmp3 = 2;
  	for ( i=0; i<10; i++ ) {
    	f_tmp1 = - f_tmp1 * f_x;
    	f_tmp2 = f_tmp1 / f_tmp3;
    	f_val += f_tmp2;
    	f_tmp3++;
  	}
  
  	return f_val;
}


double app_math_log( double f_org )
{
  	double	f_param, f_val;
  	int	i, j;
	
	// 1e30 以上の場合，Inf と見なし可能な限り大きい数を返す（本来 Inf）
	if( f_org >= 1e30 ){
		return 1e30;
	}
	
  	// 0 以下の場合，可能な限り負に大きい数を返す（本来 -Inf）．
  	if( f_org <= 0 ){
    	return -1e30;
  	}

  	// Taylor 級数で正しく求まる 1 付近になるまで平方根を取る．
  	for ( i=0, f_param = f_org; f_param > 1.3; i++ ){
    	f_param = app_math_sqrt( f_param );
  	}

  	f_val = log_base( f_param );

  	// 平方根を取った回数だけ２倍する．
  	for ( j=0; j<i; j++ ) {
    	f_val *= 2;
  	}

  	return f_val;
}


double app_math_log10( double f_org )
{
  	double	f_val;

  	f_val = app_math_log( f_org );
  	f_val /= MATH_LN10;

  	return f_val;
}


double 	app_math_exp( double f_org )
{
	int		i_org_sign;
	double	f_org_abs;
	int		i_times;
	double	f_base;
  	double	f_val, f_tmp1, f_tmp2;
  	double	f_i;
  	int		i;

	// 指数の絶対値，正負
	if( f_org < 0 ){
		i_org_sign = -1;
		f_org_abs = f_org * -1;
	}
	else{
		i_org_sign = 1;
		f_org_abs = f_org;
	}
	
	// 値が１以下となるまで，指数を２で除算
	for( f_base = f_org_abs, i_times = 0; f_base > 1; i_times++ ){
		f_base /= 2;
	}

	// Taylor 級数で exp を計算
  	f_val = 1;
  	f_tmp1 = 1;
  	f_tmp2 = 1;
  	for ( f_i = 1; f_i < 12; f_i++ ){
		f_tmp1 *= f_base;	// Taylor 級数の分子
		f_tmp2 *= f_i;		// Taylor 級数の分母
    	f_val += f_tmp1 / f_tmp2;
  	}
  
  	// ２で除算した回数だけ自乗
  	for( i=0; i<i_times; i++ ){
		f_val *= f_val;
  	}
  
  	// 指数が負の場合，逆数にする
	if( i_org_sign < 0 ){
		f_val = 1 / f_val;
	}
  
  	return f_val;
}


double app_math_pow( double f_base, double f_exp )
{
	// x^y の場合，exp( y * log( x ) ) を返す．
	return app_math_exp( f_exp * app_math_log( f_base ) );
}


double app_math_sqrt( double f_org )
{
  	double	f_x, f_s, f_last;
  	int	i;

	f_x = f_org;

	if( f_x <= 0.0 ) return 0.0;
	if( f_x > 1 )	f_s = f_x;
	else			f_s = 1;
	
	i = 0;
	do {
		i++;
		f_last = f_s;
		f_s = ( f_x / f_s + f_s ) * 0.5;
	} while( f_s < f_last );

  
  	return f_last;
}


/*
double app_math_rad2deg( double f_rad )
{
  	double f_deg;

  	f_deg = f_rad * 180 / MATH_PI;
  	f_deg -= (int)( f_deg / 360 ) * 360;
  	if ( f_deg < 0 ) f_deg += 360;
  
  	return f_deg;
}
*/


double app_math_regulate_rad( double f_rad )
{
	f_rad -= (double)( (int)( f_rad / ( MATH_PI * 2  ) ) ) * MATH_PI * 2;
	
	if ( f_rad < 0 ){
		f_rad += MATH_PI * 2;
	}

  	return f_rad;
}