Ruby 1.9.2p136(2010-12-25revision30365)

- Main Page
- Modules
- Data Structures
- Files

- File List
- Globals

# missing/erf.c

Go to the documentation of this file.

    00001 /* erf.c  - public domain implementation of error function erf(3m)
    00002
    00003 reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
    00004             (New Algorithm handbook in C language) (Gijyutsu hyouron
    00005             sha, Tokyo, 1991) p.227 [in Japanese]                 */
    00006 #include <stdio.h>
    00007 #include <math.h>
    00008
    00009 #ifdef _WIN32
    00010 # include <float.h>
    00011 # if !defined __MINGW32__ || defined __NO_ISOCEXT
    00012 #  ifndef isnan
    00013 #   define isnan(x) _isnan(x)
    00014 #  endif
    00015 #  ifndef isinf
    00016 #   define isinf(x) (!_finite(x) && !_isnan(x))
    00017 #  endif
    00018 #  ifndef finite
    00019 #   define finite(x) _finite(x)
    00020 #  endif
    00021 # endif
    00022 #endif
    00023
    00024 static double q_gamma(double, double, double);
    00025
    00026 /* Incomplete gamma function
    00027    1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt  */
    00028 static double p_gamma(double a, double x, double loggamma_a)
    00029 {
    00030     int k;
    00031     double result, term, previous;
    00032
    00033     if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
    00034     if (x == 0)     return 0;
    00035     result = term = exp(a * log(x) - x - loggamma_a) / a;
    00036     for (k = 1; k < 1000; k++) {
    00037         term *= x / (a + k);
    00038         previous = result;  result += term;
    00039         if (result == previous) return result;
    00040     }
    00041     fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__);
    00042     return result;
    00043 }
    00044
    00045 /* Incomplete gamma function
    00046    1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt  */
    00047 static double q_gamma(double a, double x, double loggamma_a)
    00048 {
    00049     int k;
    00050     double result, w, temp, previous;
    00051     double la = 1, lb = 1 + x - a;  /* Laguerre polynomial */
    00052
    00053     if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
    00054     w = exp(a * log(x) - x - loggamma_a);
    00055     result = w / lb;
    00056     for (k = 2; k < 1000; k++) {
    00057         temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
    00058         la = lb;  lb = temp;
    00059         w *= (k - 1 - a) / k;
    00060         temp = w / (la * lb);
    00061         previous = result;  result += temp;
    00062         if (result == previous) return result;
    00063     }
    00064     fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__);
    00065     return result;
    00066 }
    00067
    00068 #define LOG_PI_OVER_2 0.572364942924700087071713675675 /* log_e(PI)/2 */
    00069
    00070 double erf(double x)
    00071 {
    00072     if (!finite(x)) {
    00073         if (isnan(x)) return x;      /* erf(NaN)   = NaN   */
    00074         return (x>0 ? 1.0 : -1.0);   /* erf(+-inf) = +-1.0 */
    00075     }
    00076     if (x >= 0) return   p_gamma(0.5, x * x, LOG_PI_OVER_2);
    00077     else        return - p_gamma(0.5, x * x, LOG_PI_OVER_2);
    00078 }
    00079
    00080 double erfc(double x)
    00081 {
    00082     if (!finite(x)) {
    00083         if (isnan(x)) return x;      /* erfc(NaN)   = NaN      */
    00084         return (x>0 ? 0.0 : 2.0);    /* erfc(+-inf) = 0.0, 2.0 */
    00085     }
    00086     if (x >= 0) return  q_gamma(0.5, x * x, LOG_PI_OVER_2);
    00087     else        return  1 + p_gamma(0.5, x * x, LOG_PI_OVER_2);
    00088 }
    00089