BEGIN { # 使用例
    print gxp(1.96)           # 正規分布の上側確率        引数：Ｚ値
    print xxp(1.96*1.96, 1)   # カイ二乗分布の上側確率　　引数：カイ二乗値，自由度
    print txp(1.96, 300)      # ｔ分布の両側確率　　　　　引数：ｔ値，自由度
    print fxp(1, 1, 1)        # Ｆ分布の上側確率　　　　　引数：Ｆ値，自由度，自由度
}

# x の絶対値を求める
# 使用例： x = abs(-10.95)

function abs(x)
{
   return (x < 0) ? -x : x
}

# カイ二乗分布の上側確率

function xxp(x, is,     i, t, w, pi2)
{
   pi2 = 0.398942280401432677940
   if (x == 0.0) {
       return 1
   }
   else if (is == 1) {
       return gxp(sqrt(x))*2.0
   }
   else if (is == 2) {
       return exp(-x/2.0)
   }
   else if ((is%2) == 0) {
       w = t = 1.0
       for (i = 2; i <= is-2; i += 2) {
           t *= x/i
           w += t
       }
       return exp(log(w) - x*0.5)
   }
   else {
       t = w = sqrt(x)
       for (i = 3; i <= is-2; i += 2) {
           t *= x/i
           w += t
       }
       return 2.0*(gxp(sqrt(x))+pi2*exp(log(w)-x*0.5))
   }
}

# 正規分布の上側確率

function gxp(x,    is, i, c, p, y, z, pi2)
{
   pi2 = 0.398942280401432677940
   
   is = -1
   y = abs(x)
   c = y*y
   p = 0.0
   z = exp(-c*0.5)*pi2
   if (y < 2.5) {
       for (i = 20; i > 0; i--) {
           p = i*c/(i*2+1+is*p)
           is = -is
       }
       p = 0.5-z*y/(1.0-p)
   }
   else {
       for (i = 20; i > 0; i--) {
           p = i/(y+p)
       }
       p = z/(y+p)
   }
   return (x < 0.0) ? 1.0-p : p
}

# ｔ分布の両側確率　<a name="txp">txp.awk</a>

function txp(t, df,    i, t1, t2, p, w, m2pi)
{
   m2pi = 0.636619772367581343076
   t1 = abs(t)/sqrt(df)
   t2 = 1.0/(1.0+t1*t1)
   if ((df%2) == 0) {
       w = t1*sqrt(t2)
       p = 1.0-w
       for (i = 2; i <= df-2; i += 2) {
           p -= w *= t2*(i-1)/i
       }
   }
   else {
       p = 1.0-m2pi*atan2(t1, 1)
       if (df >= 3) {
           w = m2pi*t1*t2
           p -= w
           for (i = 3; i <= df-2; i += 2) {
               p -= w *= t2*(i-1)/i
           }
       }
   }
   return (p < 0.0 && abs(p) < 1e-10) ? 0.0 : p
}

# Ｆ分布の上側確率　<a name="fxp">fxp.awk</a>

function fxp(f, df1, df2,     i1, i2, c, p, s, w, y, m2pi)
{
   m2pi = 0.636619772367581343076
   m_pi = 3.14159265358979323846
   i1 = 2-(df1%2)
   i2 = 2-(df2%2)
   w = df2/(df2+df1*f)
   y = 1.0-w
   p = sqrt(w)
   s = sqrt(y)
   if (2*i1-i2 == 0) {
       p = 1.0-s
       c = s*w/2.0
   }
   else if (2*i1-i2 == 1) {
       c = s*p/m_pi
       p = 1.0-atan2(s, p)*m2pi
   }
   else if (2*i1-i2 == 2) {
       p = w
       c = w*y
   }
   else {
       c = y*p/2.0
   }
   for (; df2 > i2; i2 += 2) {
       p -= 2.0/i2*c
       c *= w*(i1+i2)/i2
   }
   for (; df1 > i1; i1 += 2) {
       p += 2.0/i1*c
       c *= y*(i1+i2)/i1
   }
   return (p < 0.0 && abs(p) < 1e-10) ? 0.0 : p
}