Lib/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;
}

---