/* vi: set tabstop=4 nocindent noautoindent: */

/**
数学用ライブラリ    mt.c

  ヘッダ
     #include "mt.h"

  set tabstop = 4

 */




#include "mt.h"


/**
double  power(float x, float y)

    機能: xの y乗を計算する．
 */
double  power(float x, float y)
{
    if (y==0.0)      return 1.0;
    else if (x==0.0) return 0.0;
    else return (exp((y)*log(x)));
}




/**
float  fact(int n)

    機能: nの階乗を計算する．n>0 でない場合は 1.0 を返す． 
*/
float  fact(int n)
{
	int   i;
	float ret = 1.0;

    if (n>0) for (i=1; i<=n; i++)  ret = ret*(float)i;
	return ret;
}
    



/**
float  perm(int n, int m)

    機能: パームテーション nPm を計算する．

*/
float  perm(int n, int m)
{
	int   i;
    float p;

    if (m>n || n<=0 || m<=0) return 0.;

	p = 1.0;
	for (i=n-m+1; i<=n; i++) p = p*(float)i;

    return p;
}




/**
float  comb(int n, int m)

    機能: コンビネーション nCm を計算する． 

*/
float  comb(int n, int m)
{
	int i;
    float c, p;

    if (m>n || n<=0 || m<=0) return 0.;

	p = c = 1.0;
	for (i=n-m+1; i<=n; i++) p = p*(float)i;
    for (i=1;     i<=m; i++) c = c*(float)i;

    return p/c;
}





#define  MAX_ITRTN_NEWTON_METHOD  30


int    newton_method(PTR_DFFUNC func, PTR_DFFUNC dfunc, double* t0, double eps)
{
	double tt  = *t0;
	double ff  = (*func)(tt);

	int n = 1;
	while (n<MAX_ITRTN_NEWTON_METHOD && Xabs(ff)>eps) {
		//
		double df = (*dfunc)(tt);
		if (Xabs(df)<eps) return 0;

		tt = tt - ff/df;
		ff = (*func)(tt);
		n++;
	}

	if (n==MAX_ITRTN_NEWTON_METHOD) return 0;
	*t0 = tt;
	return n;
}