psychlops/core/math/okumura/igamma.cpp

   1 /***********************************************************\r
   2         igamma.c -- \95s\8a®\91S\83K\83\93\83}\8aÖ\90\94\r
   3 ***********************************************************/
   4 #include <math.h>\r
   5 #include "igamma.h"
   6 #include "../../ApplicationInterfaces/psychlops_code_exception.h"\r
   7 \r
   8 /************ loggamma(x) -- gamma.c \82æ\82è\8dÄ\8cf *************/\r
   9 \r
  10 #define PI      3.14159265358979324  /* $\pi$ */\r
  11 #define LOG_2PI 1.83787706640934548  /* $\log 2\pi$ */\r
  12 #define N       8\r
  13 \r
  14 #define B0  1                 /* \88È\89º\82ÍBernoulli\90\94 */\r
  15 #define B1  (-1.0 / 2.0)\r
  16 #define B2  ( 1.0 / 6.0)\r
  17 #define B4  (-1.0 / 30.0)\r
  18 #define B6  ( 1.0 / 42.0)\r
  19 #define B8  (-1.0 / 30.0)\r
  20 #define B10 ( 5.0 / 66.0)\r
  21 #define B12 (-691.0 / 2730.0)\r
  22 #define B14 ( 7.0 / 6.0)\r
  23 #define B16 (-3617.0 / 510.0)\r
  24 \r
  25 double loggamma(double x)  /* \83K\83\93\83}\8aÖ\90\94\82Ì\91Î\90\94 */\r
  26 {\r
  27         double v, w;\r
  28 \r
  29         v = 1;\r
  30         while (x < N) {  v *= x;  x++;  }\r
  31         w = 1 / (x * x);\r
  32         return ((((((((B16 / (16 * 15))  * w + (B14 / (14 * 13))) * w\r
  33                     + (B12 / (12 * 11))) * w + (B10 / (10 *  9))) * w\r
  34                     + (B8  / ( 8 *  7))) * w + (B6  / ( 6 *  5))) * w\r
  35                     + (B4  / ( 4 *  3))) * w + (B2  / ( 2 *  1))) / x\r
  36                     + 0.5 * LOG_2PI - log(v) - x + (x - 0.5) * log(x);\r
  37 }\r
  38 \r
  39 /**********************************************************/\r
  40 \r
  41 double q_gamma(double a, double x, double loggamma_a);\r
  42         /* \90é\8c¾\82¾\82¯. \8eÀ\8dÛ\82Ì\92è\8b`\82Í\8cã. */\r
  43 \r
  44 double p_gamma(double a, double x, double loggamma_a)  /* \96{\95¶\8eQ\8fÆ */\r
  45 {\r
  46         int k;\r
  47         double result, term, previous;\r
  48 \r
  49         if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);\r
  50         if (x == 0)     return 0;\r
  51         result = term = exp(a * log(x) - x - loggamma_a) / a;\r
  52         for (k = 1; k < 1000; k++) {\r
  53                 term *= x / (a + k);\r
  54                 previous = result;  result += term;\r
  55                 if (result == previous) return result;\r
  56         }\r
  57         throw new Psychlops::Exception("p_gamma(): the sequence is not convergent.");\r
  58         return result;\r
  59 }\r
  60 \r
  61 double q_gamma(double a, double x, double loggamma_a)  /* \96{\95¶\8eQ\8fÆ */\r
  62 {\r
  63         int k;\r
  64         double result, w, temp, previous;\r
  65         double la = 1, lb = 1 + x - a;  /* Laguerre\82Ì\91½\8d\80\8e® */\r
  66 \r
  67         if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);\r
  68         w = exp(a * log(x) - x - loggamma_a);\r
  69         result = w / lb;\r
  70         for (k = 2; k < 1000; k++) {\r
  71                 temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;\r
  72                 la = lb;  lb = temp;\r
  73                 w *= (k - 1 - a) / k;\r
  74                 temp = w / (la * lb);\r
  75                 previous = result;  result += temp;\r
  76                 if (result == previous) return result;\r
  77         }\r
  78         throw new Psychlops::Exception("q_gamma(): the sequence is not convergent.");\r
  79         return result;\r
  80 }\r
  81 \r
  82 double p_chisq(double chisq, int df)  /* \83J\83C2\8fæ\95ª\95z\82Ì\89º\91¤\8am\97¦ */\r
  83 {\r
  84         return p_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));\r
  85 }\r
  86 \r
  87 double q_chisq(double chisq, int df)  /* \83J\83C2\8fæ\95ª\95z\82Ì\8fã\91¤\8am\97¦ */\r
  88 {\r
  89         return q_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));\r
  90 }\r
  91 \r
  92 #define LOG_PI 1.14472988584940017  /* $\log_e \pi$ */\r
  93 \r
  94 double erf(double x)  /* Gauss\82Ì\8cë\8d·\8aÖ\90\94 ${\rm erf}(x)$ */\r
  95 {\r
  96         if (x >= 0) return   p_gamma(0.5, x * x, LOG_PI / 2);\r
  97         else        return - p_gamma(0.5, x * x, LOG_PI / 2);\r
  98 }\r
  99 \r
 100 double erfc(double x)  /* $1 - {\rm erf}(x)$ */\r
 101 {\r
 102         if (x >= 0) return  q_gamma(0.5, x * x, LOG_PI / 2);\r
 103         else        return  1 + p_gamma(0.5, x * x, LOG_PI / 2);\r
 104 }\r
 105 \r
 106 double p_normal(double x)  /* \95W\8f\80\90³\8bK\95ª\95z\82Ì\89º\91¤\8am\97¦ */\r
 107 {\r
 108         if (x >= 0) return\r
 109                 0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));\r
 110         else return\r
 111                 0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);\r
 112 }\r
 113 \r
 114 double q_normal(double x)  /* \95W\8f\80\90³\8bK\95ª\95z\82Ì\8fã\91¤\8am\97¦ */\r
 115 {\r
 116         if (x >= 0) return\r
 117                 0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);\r
 118         else return\r
 119                 0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));\r
 120 }\r