mt.c

Go to the documentation of this file.

 
#include "mt.h"
 
double   EPS = 1.0e-6;
 
 
double  power(double x, double y)
{
    if (y==0.0)      return 1.0;
    else if (x==0.0) return 0.0;
    else return exp((y)*log(x));
}
 
 
double  fact(int n)
{
    int    i;
    double ret = 1.0;
 
    if (n>0) for (i=1; i<=n; i++)  ret = ret*(double)i;
    return ret;
}
    
 
double  perm(int n, int m)
{
    int    i;
    double p;
 
    if (m>n || n<=0 || m<=0) return 0.;
 
    p = 1.0;
    for (i=n-m+1; i<=n; i++) p = p*(double)i;
 
    return p;
}
 
 
double  comb(int n, int m)
{
    int    i;
    double 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*(double)i;
    for (i=1;     i<=m; i++) c = c*(double)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;
}