


#ifndef __dkutil_Math__ // include once
#define __dkutil_Math__

#include <dkutil/dKingyoMacro.h>

#define MINMAX_SWAP(type,min,max) if(min>max){type temp=min;min=max;max=temp;/*min^=max^=min^=max;*/}
#define MINMAX_ABS(min,max) if(min<0){min=abs(min);}if(max<0){max=abs(max);}
#define MINMAX_SAFETY(min,max) if(min == max){max++;}

#define MINMAX_SWAP_FLOAT(type,min,max) if(min>max){type t=min;min=max;max=t;}
#define MINMAX_ABS_FLOAT(type,min,max) if(min<0){min=(type)fabs(min);}if(max<0){max=(type)fabs(max);}

#ifndef PI
#define PI 3.14159265f
#define DKUTIL_PI_SETED
#endif

namespace dkutil{//begin dkutil namespace
	namespace oldmath{//begin math namespace
/*!
等差数列<br>
A Arithimetic Series<br>
@param array [out] = 配列のポインタを代入
@param num [in] = 配列の数
@param a [in] = 初項
@param d [in] = 公差
@return 最後の配列の中身
*/
extern int CreateArthimeticProgression(int *array,DWORD num,int a,int d);

/*!
等比数列<br>
A Geometric Series<br>
実装してから思った…テンプレートにすれば良かった(爆)<br>
調和数列も似たようなものだから削除。どうせ使わない。数列関連ってお遊び関数だし…<br>
@param array [out] = 配列のポインタを代入
@param num [in] = 配列の数
@param a [in] = 初項
@param r [in] = 公比
@return 最後の配列の中身
*/
extern int CreateGeometricProgression(int *array,DWORD num,int a,int r);
	
/*!
乱数列
A Random Number Series<br>
どうせ誰も使わないと思うが乱数列を作成補助するやつ<br>
引数の左はmin右がmax…ってかmin<=x<=maxまでを返す関数のポインタをぶち込む。<br>
@param  array [out] = 配列のポインタを代入
@param num [in] = 配列の数
@param min [in] = 乱数の最小値
@param max [in] = 乱数の最大値
@param rand_func_callback[in] = ランダム変域を返す関数。
@return 最後の配列の中身
*/
extern int CreateRandomNumberProgression(
		int *array,
		DWORD num,
		int min,int max,
		int (*rand_func_callback)(int min_,int max_));

///二つの２次元座標から角度を求める。
///ちなみに256度とかUSHORT_MAX度とかは算出できないと思います。
extern short GetAngle(int mx,int my,int sx,int sy);

#ifdef _MSC_VER
extern ULONGLONG __fastcall GetRDTSC();
#else
extern ULONGLONG GetRDTSC();
#endif

//#ifndef USE_DKINGYO_FOX_1_1
#if _MSC_VER >= 1000
#	pragma warning(disable:4244)
#	pragma warning(disable:4305)
#endif
/*!
@note
	出典：二反田 雅彦 氏のKittyエンジンより

  sin - テイラー展開した近似式によるサイン値を返す
  
  引数：
    float f 角度（ラジアン）
  
  戻り値：
    float サイン値


@param angle of radian
@return sine number
@note
      テイラー展開した、sin級数
      sin(x) ≒ x - (1 /    6) * x^3 + (1 /    120) * x^5
                  - (1 / 5040) * x^7 + (1 / 362880) * x^9
      
      この近似式の有効な範囲は、-π <= r <= πであるため、
      値の調整を行う
*/
DINLINE float sin(float f)
{

    
    // まず、fを-π〜πの範囲に修正する
    int i;
    if (f > PI)
		{
        i  = f / (2.0f * PI);
        f -= i * 2.0f * PI;
        if (f > PI) f -= 2.0f * PI;

    }else if (f < -PI)
		{
        i  = f / (2.0f * PI);
        f -= i * 2.0f * PI;
        if (f < -PI) f += 2.0f * PI;
    }
    
    float ftemp, f2 = f * f;
    float ret;
    
    ret    = f;
    ftemp  = f2 * f; // f^3
    ret   -= (1.0f /      6.0f) * ftemp;
    ftemp *= f2;     // f^5
    ret   += (1.0f /    120.0f) * ftemp;
    ftemp *= f2;     // f^7
    ret   -= (1.0f /   5040.0f) * ftemp;
    ftemp *= f2;     // f^9
    ret   += (1.0f / 362880.0f) * ftemp;
    
    return ret;
}

/*!
@note
	出典：二反田 雅彦 氏のKittyエンジンより

  cos - テイラー展開した近似式によるコサイン値を返す
  
  引数：
    float f 角度（ラジアン）
  
  戻り値：
    float コサイン値

@param f angle of radian
@return cosine number
*/
DINLINE float cos(float f)
{
    /*
      テイラー展開した、cos級数
      cos(x) ≒ 1 - (1 /   2) * x^2 + (1 /    24) * x^4
                  - (1 / 720) * x^6 + (1 / 40320) * x^8
      
      この近似式の有効な範囲は、-π <= r <= πであるため、
      値の調整を行う
    */
    
    // まず、fを-π〜πの範囲に修正する
    if (f < 0.0f) f *= -1.0f;
    int i = f / (2.0f * PI);
    f -= i * 2.0f * PI;
    if (f > PI) f -= 2.0f * PI;
    
    float ftemp, f2 = f * f;
    float ret;
    
    ret    = 1.0f;
    ftemp  = f2; // f^2
    ret   -= (1.0f /     2.0f) * ftemp;
    ftemp *= f2; // f^4
    ret   += (1.0f /    24.0f) * ftemp;
    ftemp *= f2; // f^6
    ret   -= (1.0f /   720.0f) * ftemp;
    ftemp *= f2; // f^8
    ret   += (1.0f / 40320.0f) * ftemp;

    return ret;
}

/*!
@note
	出典：二反田 雅彦 氏のKittyエンジンより

  asin - テイラー展開した近似式によるアークサイン値を返す
         ※誤差が大きいので、使用をお勧めできない
  
  引数：
    float f サイン値
  
  戻り値：
    float 角度（ラジアン）
@param f sine number
@return angle of radian
@note
テイラー展開した、sin-1級数<br>
      sin-1(x) ≒ x + (1 /   6) * x^3 + ( 3 /   40) * x^5<br>
                    + (5 / 112) * x^7 + (35 / 1152) * x^9<br>
*/
DINLINE float asin(float f)
{
    /*
     
    */
    //assert(!"called KM_Basic::asin!!");
    
    float ftemp, f2 = f * f;
    float ret;
    
    ret    = f;
    ftemp  = f2 * f; // f^3
    ret   += ( 1.0f /    6.0f) * ftemp;
    ftemp *= f2;     // f^5
    ret   -= ( 3.0f /   40.0f) * ftemp;
    ftemp *= f2;     // f^7
    ret   += ( 5.0f /  112.0f) * ftemp;
    ftemp *= f2;     // f^9
    ret   += (35.0f / 1152.0f) * ftemp;

    return ret;
}

/*!
@note
	出典：二反田 雅彦 氏のKittyエンジンより

  acos - テイラー展開した近似式によるアークコサイン値を返す
         ※誤差が大きいので、使用をお勧めできない
  
  引数：
    float f コサイン値
  
  戻り値：
    float 角度（ラジアン）
@param f cosine number
@return angle of radian
@note
テイラー展開した、cos-1級数<br>
      cos-1(x) ≒ (1 / 2)π -               x   - (1 /   6) * x^3<br>
                            - ( 3 /   40) * x^5 - (5 / 112) * x^7<br>
                            - (35 / 1152) * x^9<br>
*/
DINLINE float acos(float f)
{
    /*
      
    */
    //assert(!"called KM_Basic::acos!!");
    
    float ftemp, f2 = f * f;
    float ret;
    
    ret    = ( 1.0f /    2.0f) * PI - f;
    ftemp  = f2 * f; // f^3
    ret   -= ( 1.0f /    6.0f) * ftemp;
    ftemp *= f2;     // f^5
    ret   -= ( 3.0f /   40.0f) * ftemp;
    ftemp *= f2;     // f^7
    ret   -= ( 5.0f /  112.0f) * ftemp;
    ftemp *= f2;     // f^9
    ret   -= (35.0f / 1152.0f) * ftemp;
    
    return ret;
}

/*!
@note
	出典：二反田 雅彦 氏のKittyエンジンより

  sqrtf - float型の平方根を求める
          ※sqrtのfloatバージョン
            ボーランドのｺﾝﾊﾟｲﾗには、ないようなので自分で作成

　引数：
　　float f 平方根を求める値、負の値はダメよん

　戻り値：
　　float 平方根、fが負だった場合は-1を返す
@param f 平方根を求める値、負の値はＮＧ
@return 平方根 fが負だった場合は-1を返す
@note
<PRE>
実数の平方根s = √xは、f(s) = s^2 - x = 0をNewton法で解い
て求める。

f(s)  = s^2 - x
f'(s) = 2s
を

xn+1  = xn - f(xn) / f'(xn)

に代入して、
xn+1  = (x / s + s) / 2

Newton法について、詳しくは数学の本でも読んでね♪
</PRE>
*/

//#ifndef _DKUTIL_MSVC6LINKER_POLICY_H_
DINLINE float sqrtf(float f){

    float s, last;
    
    if (f > 1e-10){
        if (f > 1.0f) s = f;
        else s = 1.0f;
        
        do{
            last = s;
            s    = (f / s + s) * 0.5;
        }while (s < last);
        
        return last;
    }
    
    if (f < 0.0f) return -1.0f;
    
    return 0.0f;
}
//#endif
#if _MSC_VER >= 1000
#	pragma warning(default:4244)
#	pragma warning(default:4305)
#endif







/*!
参考文章
<ul>
<li>http://www.bio.nagoya-u.ac.jp:8001/~shionyu/CompJHP/ex/ex4-4prog.html
<li>http://www.bio.nagoya-u.ac.jp:8001/~shionyu/CompJHP/func.html
<li>http://www.d4.dion.ne.jp/~sekiya_z/na/na16mon1.html
<li>http://www.geocities.co.jp/HeartLand-Asagao/6686/statistic.html
<li>http://www.bb.wakwak.com/~bird/private/reseach/analysis.html
<li>http://www.wombat.zaq.ne.jp/nao/c/
<li>http://chiharu.cside4.jp/elearn/stat2002/INDEX.HTM
</ul>
	@Note
カイ二乗分布のパーセント点を求める。
χ2　=　 (Ｏy−Ｅy）２／Ｅy＋(Ｏg−Ｅg）２／Ｅg.
χ2　=　 Σ(Ｏi−Ｅi）２／Ｅi
危険率を１％に下げるにはカイ二乗の値が６.６３以上かを調べればよい。

d金魚曰、
まだ、実験段階中なので使わないでください〜〜〜
*/
///カイ二乗検定らしい
extern double kai2(double o1,double o2,double e1,double e2);
extern float kai2(float o1,float o2,float e1,float e2);


}//end of math namespace
}//end of dkutil namespace


#ifdef DKUTIL_PI_SETED
#undef PI
#endif

#endif//end of include once