// pval.C  Version 7.42  Copyleft (c) 杉本典夫
// Time-stamp: <2025-11-03 14:58:19 tac>
//
// 各種確率分布のp値変換と逆変換
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W pval pval.C -I/usr/local/include
//               -L/usr/local/lib -ldans -lfun

#include <dans.H>               // dans用ヘッダファイル(fun.H含む)

using namespace dans;           // dans名前空間使用(名前空間matvec含む)

const vecstr USAGE =            // プログラム使用法表示用
{
  "pval  Version 7.42           Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:pval",
  "機能:各種確率分布のp値変換・逆変換を行います。",
  "オプション:なし"
};

const string Title = "=== 各種確率分布のp値変換・逆変換 ==="; // 表題
constexpr double Emax = 100.0; // 項の指数最大値(オーバーフロー防止用)
constexpr double Tmax = 1.0e+100; // 項最大値(オーバーフロー防止用)
constexpr double Pmin = FDmin * 1.0e+100; // 確率最小値
constexpr double Peps = 1.0e-10; // 収束限界
constexpr double Delta  = Peps * 1.0e+5; // きざみ幅δ
constexpr double Rerror = 1.0e-8; // 収束限界
constexpr double Feps = FDeps * 10.0; // F計算誤差
constexpr int DFlarg = 100000;        // 自由度最大値

constexpr double Xvec[16] =     // Legendre-Gauss積分の積分点
{
  0.989400934991649932596154173450,
  0.944575023073232576077988415534,
  0.865631202387831743880467897712,
  0.755404408355003033895101194847,
  0.617876244402643748446671764048,
  0.458016777657227386342419442983,
  0.281603550779258913230460501460,
  0.095012509837637440185319335424,
  -0.095012509837637440185319335424,
  -0.281603550779258913230460501460,
  -0.458016777657227386342419442983,
  -0.617876244402643748446671764048,
  -0.755404408355003033895101194847,
  -0.865631202387831743880467897712,
  -0.944575023073232576077988415534,
  -0.989400934991649932596154173450
};

constexpr double Wvec[16] =     // Legendre-Gauss積分の重み係数
{
  0.027152459411754094851780572456,
  0.062253523938647892862843836994,
  0.095158511682492784809925107602,
  0.124628971255533872052476282192,
  0.149595988816576732081501730547,
  0.169156519395002538189312079030,
  0.182603415044923588866763667969,
  0.189450610455068496285396723208,
  0.189450610455068496285396723208,
  0.182603415044923588866763667969,
  0.169156519395002538189312079030,
  0.149595988816576732081501730547,
  0.124628971255533872052476282192,
  0.095158511682492784809925107602,
  0.062253523938647892862843836994,
  0.027152459411754094851780572456
};

double fgamma(int df1, int df2) // F分布のガンマ関数
{
  if ((df1 % 2) == 0 || df1 == 1 || (df2 % 2) == 0) return 0.0;
  double term = 2.0 / M_PI;
  int dfhalf = (df2 - 1) / 2;
  for (int i = 1; i <= dfhalf; i++) term *= (double)i / ((double)i - 0.5);
  return term;
}

double even1(double f, int df1, int df2) // df1:偶数, df2:偶数又は奇数
{
  double x = (double)df2 / ((double)df1 * f + (double)df2);
  double sum = 1.0, term = 1.0;
  double tc = pow(x, Emax);
  double dfc = 0.0;
  for (int i = 2; i <= df1 - 2; i += 2)
  {
    term *= (1.0 - x) * (double)(df2 + i - 2) / (double)i;
    if (term > Tmax) { term *= tc;  sum *= tc;  dfc += Emax; }
    sum += term;
  }
  if (log(sum) + ((double)df2 / 2.0 - dfc) * log(x) <= log(Feps)) return Feps;
  return sum * pow(x, (double)df2 / 2.0 - dfc);
}

double oddeven(double f, int df1, int df2) // df1:奇数,df2:偶数
{
  double x = (double)df2 / ((double)df1 * f + (double)df2);
  double sum = 1.0, term = 1.0;
  double tc = pow(1.0 - x, Emax);
  double dfc = 0.0;
  for (int i = 2; i <= df2 - 2 && term > Feps; i += 2)
  {
    term *= x * (double)(df1 + i - 2) / (double)i;
    if (term > Tmax) { term *= tc;  sum *= tc;  dfc += Emax; }
    sum += term;
  }
  return 1.0 - sum * pow(1.0 - x, (double)df1 / 2.0 - dfc);
}

double oddodd(double f, int df1, int df2) // df1:奇数,df2:奇数
{
  double theta = atan(sqrt((double)df1 * f / (double)df2));
  double sin_t = sin(theta);
  double cos_t = cos(theta);
  double sum = 0.0, term = 1.0, theta2 = 1.0;
  if (df2 != 1)
  {
    sum = term = cos_t * sin_t;  theta2 = cos_t * cos_t;
    for (int i = 2; i <= df2 - 3 && term > Feps; i += 2)
    {
      term *= theta2 * (double)i / (double)(i + 1);
      sum += term;
    }
  }
  double a = (theta + sum) * 2.0 / M_PI;
  double b = 0.0;
  if (df1 != 1)
  {
    sum = term = fgamma(df1, df2) * sin_t;  theta2 = sin_t * sin_t;
    double tc = pow(cos_t, Emax);
    double dfc = 0;
    for (int i = 2; i <= df1 - 3; i += 2)
    {
      term *= theta2 * (double)(df2 + i - 1) / (double)(i + 1);
      if (term > Tmax) { term *= tc;  sum *= tc;  dfc += Emax; }
      sum += term;
    }
    if (log(sum) + ((double)df2 - dfc) * log(cos_t) <= log(Feps)) b = Feps;
    else b = sum * pow(cos_t, (double)df2 - dfc);
  }
  return 1.0 - a + b;
}

double ftop_df(double f, int df1, int df2) // F値→上側p値変換
{
  if (f < Feps || df1 < 1 || df2 < 1) return 1.0;
  double pr = 1.0;
  switch ((df1 % 2) * 10 + (df2 % 2))
  {
  case 0  :
  case 1  : pr = even1(f, df1, df2);
    break;
  case 10 : pr = oddeven(f, df1, df2);
    break;
  case 11 : pr = oddodd(f, df1, df2);
    break;
  }
  return fabs(pr);              // 計算誤差により負の数になる時があるため
}

double qdisinf(double q, int g)
// スチューデント化された範囲のq分布(自由度無限大)
// 理論式に忠実な方法のテスト用(pow(pr, double(g)を利用しない)
{
  constexpr double Qlarg = 20.0;
  constexpr size_t Qdiv = 10;
  constexpr size_t Nleg = 16;
  constexpr double Eps1 = -30.0;
  constexpr double Eps2 = -50.0;
  constexpr double Eps3 = 60.0;
  constexpr double Sqrt2pi = 1.0 / sqrt(2.0 * M_PI);

  double qsqz = q * 0.5;
  double zupl = 0.0;
  double zh = (Qlarg - zupl) / double(Qdiv);
  double zupu = zupl + zh;
  double zlowu = zupl;
  double zlowl = zlowu - zh;
  double prlsum = 0.0, prusum = 0.0;
  double g1 = double(g - 1);
  bool cnt = true;
  for (size_t i = 0; cnt && (i < Qdiv); ++i)
  {
    double prup = 0.0, prlow = 0.0;
    double zupm = 0.5 * (zupu + zupl), zuphh = 0.5 * (zupu - zupl);
    double zlowm = 0.5 * (zlowu + zlowl), zlowhh = 0.5 * (zlowu - zlowl);
    for (size_t j = 0; cnt && (j < Nleg); ++j)
    {
      double zupx = zupm + zuphh * Xvec[j];
      double zupx2 = zupx * zupx;
      if (zupx2 > Eps3)
      { cnt = false; }
      else
      {
        double nedpi = nedp(zupx) - nedp(zupx - q);
        prup += Wvec[j] * exp(-0.5 * zupx2) * pow(nedpi, g1);
        double zlowx = zlowm + zlowhh * Xvec[j];
        double zlowx2 = zlowx * zlowx;
        nedpi = nedp(zlowx) - nedp(zlowx - q);
        prlow += Wvec[j] * exp(-0.5 * zlowx2) * pow(nedpi, g1);
      }
    }
    if (cnt)
    {
      prusum += zuphh * (double)g * Sqrt2pi * prup;
      prlsum += zlowhh * (double)g * Sqrt2pi * prlow;
      zupl = zupu;  zlowu = zlowl;
      zupu += zh;  zlowl -= zh;
    }
  }
  double pr = prlsum + prusum;
  if (pr < 0.0) pr = -pr;
  if (pr > 1.0) return 1.0;
  return pr;
}

double qdis(double q, int g, int df) // スチューデント化された範囲のq分布
{
  constexpr size_t Nmax = 50;
  constexpr size_t Nlegq = 16;
  constexpr double Eps1 = -30.0;
  constexpr double Eps2 = 1.0e-14;
  constexpr int DFhaf = 100;
  constexpr int DFquar = 800;
  constexpr int DFeigh = 5000;
  constexpr double Ulen1 = 1.0;
  constexpr double Ulen2 = 0.5;
  constexpr double Ulen3 = 0.25;
  constexpr double Ulen4 = 0.125;

  double f2 = (double)df * 0.5;
  double f2lf = f2 * log((double)df) - (double)df * log(2.0) - lgamma(f2);
  double f21 = f2 - 1.0;
  double ff4 = (double)df * 0.25;
  double ulen = Ulen4;
  if (df <= DFhaf) ulen = Ulen1;
  else if (df <= DFquar) ulen = Ulen2;
  else if (df <= DFeigh) ulen = Ulen3;
  f2lf += log(ulen);
  double pr = 0.0;
  for (size_t i = 1; i <= Nmax; ++i)
  {
    double otsum = 0.0;
    double twa1 = (double)(2 * i - 1) * ulen;
    for (size_t j = 0; j < Nlegq; ++j)
    {
      double twu = twa1 + Xvec[j] * ulen;
      double t1 = f2lf + f21 * log(twu) - twu * ff4;
      double qsqz = q * sqrt(twu * 0.5);
      otsum += Wvec[j] * exp(t1) * qdisinf(qsqz, g);
    }
    if ((double(i) * ulen >= 1.0) && (otsum < Eps2)) break; // 計算終了
    pr += otsum;
  }
  if (pr > 1.0) return 1.0;
  return pr;
}

double qntop(double q, int g, int df = 0)
// スチューデント化された範囲のq値→上側p値変換
{
  if (g == 2)
  {
    double t = q / sqrt(2.0);
    return (df > 0) ? ttop(t, df) : ztop(t);
  }
  if (df < 1 || df >= DFlarg) return 1.0 - qdisinf(q, g); // df=∞
  return 1.0 - qdis(q, g, df);
}

double ddisinflow(double d, int g, double r) // 多変量正規分布の下側確率近似分布
{
  constexpr double Dlarg = 20.0;
  constexpr size_t Ddiv = 10;
  constexpr size_t Nleg = 16;
  constexpr double Eps1 = -30.0;
  constexpr double Eps2 = -50.0;
  constexpr double Eps3 = 60.0;
  constexpr double Sqrt2pi = 1.0 / sqrt(2.0 * M_PI);

  double g1 = (double)(g - 1);
  double zupl = 0.0;
  double zh = (Dlarg - zupl) / (double)Ddiv;
  double zupu = zupl + zh;
  double zlowu = zupl;
  double zlowl = zlowu - zh;
  double sqr = sqrt(fabs(r)), sqr1 = sqrt(1.0 - fabs(r));
  double sqrr = sqr / sqr1;
  double rd = d / sqr1;
  double prlsum = 0.0, prusum = 0.0;
  bool cnt = true;
  for (size_t i = 0; cnt && (i < Ddiv); ++i)
  {
    double prup = 0.0, prlow = 0.0;
    double zupm = 0.5 * (zupu + zupl), zuphh = 0.5 * (zupu - zupl);
    double zlowm = 0.5 * (zlowu + zlowl), zlowhh = 0.5 * (zlowu - zlowl);
    for (size_t j = 0; cnt && (j < Nleg); ++j)
    {
      double zupx = zupm + zuphh * Xvec[j];
      double zupx2 = zupx * zupx;
      if (zupx2 > Eps3)
      { cnt = false; }
      else
      {
        double rzupx = sqrr * zupx;
        double nedpi = nedp(rzupx + rd);
        prup += Wvec[j] * exp(-0.5 * zupx2) * pow(nedpi, g1);
        double zlowx = zlowm + zlowhh * Xvec[j];
        double zlowx2 = zlowx * zlowx;
        double rzlowx = sqrr * zlowx;
        nedpi = nedp(rzlowx + rd);
        prlow += Wvec[j] * exp(-0.5 * zlowx2) * pow(nedpi, g1);
      }
    }
    if (cnt)
    {
      prusum += zuphh * Sqrt2pi * prup;
      prlsum += zlowhh * Sqrt2pi * prlow;
      zupl = zupu;  zlowu = zlowl;
      zupu += zh;  zlowl -= zh;
    }
  }
  double pr = prlsum + prusum;
  if (pr > 1.0) return 1.0;
  return pr;
}

double ddislow(double d, int g, int df, double r)
// 多変量t分布の下側確率近似分布
{
  constexpr size_t Nmax = 100;
  constexpr size_t Nlegq = 16;
  constexpr double Eps1 = -30.0;
  constexpr double Eps2 = 1.0e-14;
  constexpr int DFhaf = 100;
  constexpr int DFquar = 800;
  constexpr int DFeigh = 5000;
  constexpr double Hhalf1 = 0.5;
  constexpr double Hhalf2 = 0.25;
  constexpr double Hhalf3 = 0.125;
  constexpr double Hhalf4 = 0.0625;

  double hhalf = Hhalf4;
  if (df <= DFhaf) hhalf = Hhalf1;
  else if (df <= DFquar) hhalf = Hhalf2;
  else if (df <= DFeigh) hhalf = Hhalf3;
  double f2 = (double)df * 0.5;
  double f21 = f2 - 1.0;
  double ff4 = (double)df * 0.25;
  double f2lf = log(hhalf) + f2 * log((double)df) - (double)df * log(2.0) - lgamma(f2);
  double pr = 0.0;
  for (size_t i = 1; i <= Nmax; ++i)
  {
    double pri = 0.0;
    double cmed = (double)(2 * i - 1) * hhalf;
    for (size_t j = 0; j < Nlegq; ++j)
    {
      double cx = cmed + Xvec[j] * hhalf;
      double lfung = f2lf + f21 * log(cx) - cx * ff4;
      double dsx = d * sqrt(cx * 0.5);
      pri += Wvec[j] * exp(lfung) * ddisinflow(dsx, g, r);
    }
    if ((double(i) * hhalf >= 1.0) && (pri < Eps2)) // 計算終了
    {
      break;
    }
    pr += pri;
  }
  if (pr > 1.0) return 1.0;
  return pr;
}

double ddisinf(double d, int g, double r) // 多変量正規分布の両側確率近似分布
{
  constexpr double Dlarg = 20.0;
  constexpr size_t Ddiv = 10;
  constexpr size_t Nleg = 16;
  constexpr double Eps1 = -30.0;
  constexpr double Eps2 = -50.0;
  constexpr double Eps3 = 60.0;
  constexpr double Sqrt2pi = 1.0 / sqrt(2.0 * M_PI);

  double g1 = (double)(g - 1);
  double zlb = 0.0;
  double zinc = (Dlarg - zlb) / (double)Ddiv;
  double zub = zlb + zinc;
  double sqr = sqrt(fabs(r)), sqr1 = sqrt(1.0 - fabs(r));
  double sqrr = sqr / sqr1;
  double rd = d / sqr1;
  double prsum = 0.0;
  bool cnt = true;
  for (size_t i = 0; cnt && (i < Ddiv); ++i)
  {
    double zmed = 0.5 * (zub + zlb), hhalf = 0.5 * (zub - zlb);
    double pri = 0.0;
    for (size_t j = 0; cnt && (j < Nleg); ++j)
    {
      double zx = zmed + hhalf * Xvec[j];
      double zx2 = zx * zx;
      if (zx2 > Eps3)           // 計算終了
      { cnt = false; }
      else
      {
        double rzx = sqrr * zx;
        double nedpi = nedp(rzx + rd) - nedp(rzx - rd);
        pri += Wvec[j] * exp(-0.5 * zx2) * pow(nedpi, g1);
      }
    }
    if (cnt)
    {
      prsum += hhalf * Sqrt2pi * pri;
      zlb = zub;
      zub += zinc;
    }
  }
  double pr = 2.0 * prsum;
  if (pr < 0.0) pr = -pr;
  if (pr > 1.0) return 1.0;
  return pr;
}

double ddis(double d, int g, int df, double r)
// 多変量t分布の両側確率近似分布(2s2乗の代わりにs2乗を用いる方法)
{
  constexpr size_t Nmax = 50;
  constexpr size_t Nlegq = 16;
  constexpr double Eps1 = -30.0;
  constexpr double Eps2 = 1.0e-14;
  constexpr int DFhaf = 100;
  constexpr int DFquar = 800;
  constexpr int DFeigh = 5000;
  constexpr double Hhalf1 = 0.5;
  constexpr double Hhalf2 = 0.25;
  constexpr double Hhalf3 = 0.125;
  constexpr double Hhalf4 = 0.0625;

  double hhalf = Hhalf4;
  if (df <= DFhaf) hhalf = Hhalf1;
  else if (df <= DFquar) hhalf = Hhalf2;
  else if (df <= DFeigh) hhalf = Hhalf3;
  double f2 = (double)df * 0.5;
  double f21 = f2 - 1.0;
  double ff4 = (double)df * 0.25;
  double f2lf = log(hhalf) + f2 * log((double)df) - f2 * log(2.0) - lgamma(f2);
  double pr = 0.0;
  for (size_t i = 1; i <= Nmax; ++i)
  {
    double pri = 0.0;
    double cmed = (double)(2 * i - 1) * hhalf;
    for (size_t j = 0; j < Nlegq; ++j)
    {
      double cx = cmed + Xvec[j] * hhalf;
      double lfung = f2lf + f21 * log(cx) - cx * f2;
      double dsx = d * sqrt(cx);
      pri += Wvec[j] * exp(lfung) * ddisinf(dsx, g, r);
    }
    if (((double)i * hhalf >= 1.0) && (pri < Eps2)) break; // 計算終了
    pr += pri;
  }
  if (pr > 1.0) return 1.0;
  return pr;
}

double dntop(double d, int g, int df = 0, double r = 0.5, bool onet = false)
// ダネットのd値→片側/両側p値変換
{
  double rabs = fabs(r);
  if (g == 2 || rabs >= 1.0 - Feps) return onet ? ttop(d, df) * 0.5 : ttop(d, df);
  double pr = 0.0;
  if (df < 1 || df >= DFlarg)   // df=∞
  {
    pr = onet ? ddisinflow(d, g, r) : ddisinf(d, g, r);
  }
  else
  {
    pr = onet ? ddislow(d, g, df, r) : ddis(d, g, df, r);
  }
  return 1.0 - pr;
}

double guesst(double pr, int df) // t値の近似推定
{
  if (pr <= 0.0 || pr >= 1.0) return 0.0;
  if (df < 1) return guessz(pr);
  if (df == 1) return tan(M_PI * (1.0 - pr) * 0.5);
  if (df == 2) return sqrt(2.0) * tan(asin(1.0 - pr));
  double z = guessz(pr);
  return sqrt((double)df / (double)(df - 2)) * z * (1.0 + (z * z - 3.0) / 4.0 / (double)df);
}

double guessd(double pr, int g, int df, double r, bool onet) // d分布%点近似式
{
  if (pr <= 0.0 || pr >= 1.0 || g < 2) return 0.0;
  if (onet) pr *= 2.0;
  if (g == 2) return guesst(pr, df);
  double pri = 1.0 - pow(1.0 - pr, 1.0 / (double)(g - 1));
  double z0 = guesst(pri, df), z1 = guesst(pr, df);
  double w = pow(r, (double)(g - 1));
  return w * z1 + (1.0 - w) * z0;
}

double ptodn(double pr, int g, int df = 0, double r = 0.5, bool onet = false)
// 片側/両側p値→ダネットのd値逆変換
{
  if (pr <= 0.0 || pr >= 1.0 || g < 1) return 0.0;
  if (g == 2 || r <= -1 || 1.0 <= r) return onet ? ptot(pr * 2.0, df) : ptot(pr, df);
  double d = guessd(pr, g, df, r, onet);
  double ppre = dntop(d, g, df, r, onet);
  if (ppre < Pmin)              // pr初期値が非常に小さい時
  {
    d = (double)g;
    ppre = dntop(d, g, df, r, onet);
  }
  while (sgndiff(ppre, pr) != 0 && fabs(ppre - pr) >= Peps)
  {
    double dd = Delta * d;
    if (sgn(dd) == 0) dd = Delta;
    double slope = log(dntop(d + dd, g, df, r, onet) / ppre) / dd;
    if (sgn(slope) == 0) break; // 収束
    d += log(pr / ppre) / slope;
    if (sgn(d) <= 0) d = dd;
    ppre = dntop(d, g, df, r, onet);
    if (sgn(ppre) <= 0)         // 収束失敗
    {
      d = HUGE_VAL;
      break;
    }
  }
  return d;
}

double guessq(double pr, int g, int df)
// スチューデント化された範囲のq分布%点の近似推定
{
  if (pr <= 0.0 || pr >= 1.0 || g < 2) return 0.0;
  if (g == 2) return guesst(pr, df) * sqrt(2.0);
  double pri = 1.0 - pow(1.0 - pr, 1.0 / (double)(g * (g - 1) / 2));
  return guesst(pri, df) * sqrt(2.0) * 0.97;
}

double ptoqn(double pr, int g, int df = 0)
// 上側p値→スチューデント化された範囲のq値逆変換
{
  if (pr <= 0.0 || pr >= 1.0 || g < 1) return 0.0;
  if (g == 2) return ptot(pr, df) * sqrt(2.0);
  double q = guessq(pr, g, df);
  double ppre = qntop(q, g, df);
  if (ppre < Pmin)              // pr初期値が非常に小さい時
  {
    q = (double)g;
    ppre = qntop(q, g, df);
  }
  while (sgndiff(ppre, pr) != 0 && fabs(ppre - pr) >= Peps)
  {
    double dq = Delta * q;
    if (sgn(dq) == 0) dq = Delta;
    double slope = log(qntop(q + dq, g, df) / ppre) / dq;
    if (sgn(slope) == 0) break; // 収束
    q += log(pr / ppre) / slope;
    if (sgn(q) <= 0) q = dq;
    ppre = qntop(q, g, df);
    if (sgn(ppre) <= 0)         // 収束失敗
    {
      q = HUGE_VAL;
      break;
    }
  }
  return q;
}

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
};

option::option(const char* c_pname) : v_path(), v_pname() // コンストラクタ
{
  if (c_pname == nullptr) return;     // 初期設定せず
  v_pname = getname(c_pname);         // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
}

void option::set(const string& optstr, char devc)
// オプション指定によってオプションを設定する
{
  string opt[2];                // オプションのopt[0]:項目,opt[1]:内容
  if (getopt(opt, optstr, devc) < 1) return; // オプション指定なし
  int err = 0;
  if (opt[0] == "f")            // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

class pvalclass                 // pvalメインクラス
{
private:
  double get_fdata(const string& message, double low, double up, double def) const;
  // 実数を入力する
public:
  pvalclass()                   // コンストラクタ
  {}
  ~pvalclass()                  // デストラクタ
  {}
  int sel(int kind) const;      // メニュー表示と選択
  void zp() const;              // 正規分布z値→両側p値変換
  void pz() const;              // 正規分布z値←両側p値変換
  void cp() const;              // χ2乗値→上側p値変換
  void pc() const;              // χ2乗値←上側p値変換
  void tp() const;              // t値→両側p値変換
  void pt() const;              // t値←両側p値変換
  void fp() const;              // F値→上側p値変換
  void pf() const;              // F値←上側p値変換
  void dp() const;              // d値→両側p値変換
  void pd() const;              // d値←両側p値変換
  void qp() const;              // q値→上側p値変換
  void pq() const;              // q値←上側p値変換
  void np() const;              // 正規偏位→累積確率変換
  void pn() const;              // 正規偏位←累積確率変換
  void ep() const;              // 誤差関数x値→erf(x)、erfc(x)変換
  void ncclp() const;           // 非心χ2乗値→下側p値変換
  void lpncc() const;           // 非心χ2乗値←下側p値変換
  void nctlp() const;           // 非心t値→下側p値変換
  void lpnct() const;           // 非心t値変換←下側p値
  void ncflp() const;           // 非心F値→下側p値変換
  void lpncf() const;           // 非心F値←下側p値変換
  void test() const;            // テストルーチン実行
};

int pvalclass::sel(int kind) const // メニュー表示と選択
{
  puttitle(cout, Title);
  cout << "\n";
  cout << " 1:正規分布z値→両側p値\n";
  cout << " 2:正規分布z値←両側p値\n";
  cout << " 3:χ2乗値→上側p値\n";
  cout << " 4:χ2乗値←上側p値\n";
  cout << " 5:t値→両側p値\n";
  cout << " 6:t値←両側p値\n";
  cout << " 7:F値→上側p値\n";
  cout << " 8:F値←上側p値\n";
  cout << " 9:d値→両側p値\n";
  cout << "10:d値←両側p値\n";
  cout << "11:q値→上側p値\n";
  cout << "12:q値←上側p値\n";
  cout << "13:正規偏位z値→累積確率\n";
  cout << "14:正規偏位z値←累積確率\n";
  cout << "15:誤差関数x値→erf(x)、erfc(x)\n";
  cout << "16:非心χ2乗値→下側p値\n";
  cout << "17:非心χ2乗値←下側p値\n";
  cout << "18:非心t値→下側p値\n";
  cout << "19:非心t値←下側p値\n";
  cout << "20:非心F値→下側p値\n";
  cout << "21:非心F値←下側p値\n\n";
  cout << "98:テストルーチン実行\n";
  cout << "99:オンラインマニュアル表示\n";
  cout << "0:コマンド終了\n\n";
  return int(get_fdata("番号を入力してください", 0.0, 99.0, (double)kind));
}

double pvalclass::get_fdata(const string& message, double low, double up, double def) const
// 実数を入力する
{
  string sline;
  bool ok = false;
  double fdata = 0.0;
  istream* finp = nullptr;
  finp = &cin;
  do
  {
    cout << message << "(" << low << "-" << up << ",[Return]:" << def << "):";
    getsline(sline, *finp);
    strchomp(sline);
    fdata = sline.empty() ? def : strtodf(sline);
    if (fdata < low || up < fdata) beep();
    else ok = true;
  } while (!ok);
  return fdata;
}

void pvalclass::zp() const      // 正規分布z値→両側p値変換
{
  double z = 0.0;
  while ((z = get_fdata("正規分布z値(0:終了)", -FDmax, FDmax, z)) != 0.0)
  { cout << "正規分布z値=" << z << "  両側p値=" << ztop((ldbl)z) << "\n" << endl; }
}

void pvalclass::pz() const      // 正規分布z値変換←両側p値
{
  double pr = 0.0;
  while ((pr = get_fdata("両側p値(0:終了)", 0.0, 1.0, pr)) > 0.0)
  { cout << "両側p値=" << pr << "  正規分布z値=" << ptoz((ldbl)pr) << "\n" << endl; }
}

void pvalclass::cp() const      // χ2乗値→上側p値変換
{
  double csqr = 0.0;
  int df = 1;
  while ((csqr = get_fdata("χ2乗値(0:終了)", 0.0, FDmax, csqr)) > 0.0)
  {
    df = (int)(get_fdata("自由度", 1.0, (double)INT_MAX, (double)df));
    cout << "χ2乗値=" << csqr << "  自由度=" << df << "  上側p値=" << ctop((ldbl)csqr, df) << "\n" << endl;
  }
}

void pvalclass::pc() const      // χ2乗値←上側p値変換
{
  double pr = 0.0;
  int df = 1;
  while ((pr = get_fdata("上側p値(0:終了)", 0.0, 1.0, pr)) > 0.0)
  {
    df = (int)(get_fdata("自由度", 1.0, (double)INT_MAX, (double)df));
    cout << "上側p値=" << pr << "  自由度=" << df << "  χ2乗値=" << ptoc((ldbl)pr, df) << "\n" << endl;
  }
}

void pvalclass::tp() const      // t値→両側p値変換
{
  double t = 0.0;
  int df = 0;
  while ((t = get_fdata("t値(0:終了)", -FDmax, FDmax, t)) != 0.0)
  {
    df = (int)(get_fdata("自由度(0:∞)", 0.0, (double)INT_MAX, (double)df));
    if (df > 0) cout << "t値=" << t << "  自由度=" << df << "  両側p値=" << ttop((ldbl)t, df) << "\n" << endl;
    else cout << "t値=" << t << "  自由度=∞  両側p値=" << ttop((ldbl)t, df) << "\n" << endl;
  }
}

void pvalclass::pt() const      // t値←両側p値変換
{
  double pr = 0.0;
  int df = 0;
  while ((pr = get_fdata("両側p値(0:終了)", 0.0, 1.0, pr)) > 0.0)
  {
    df = (int)(get_fdata("自由度(0:∞)", 0.0, (double)INT_MAX, (double)df));
    if (df > 0) cout << "両側p値=" << pr << "  自由度=" << df << "  t値=" << ptot((ldbl)pr, df) << "\n" << endl;
    else cout << "両側p値=" << pr << "  自由度=∞  t値=" << ptot((ldbl)pr, df) << "\n" << endl;
  }
}

void pvalclass::fp() const      // F値→上側p値変換
{
  double f = 0.0;
  int df1 = 1, df2 = 1;
  while ((f = get_fdata("F値(0:終了)", 0.0, FDmax, f)) > 0.0)
  {
    df1 = (int)(get_fdata("第1自由度", 1.0, (double)INT_MAX, (double)df1));
    df2 = (int)(get_fdata("第2自由度", 1.0, (double)INT_MAX, (double)df2));
    cout << "F値=" << f << "  第1自由度=" << df1 << "  第2自由度=" << df2 << "  上側p値=" << ftop((ldbl)f, df1, df2) << "\n";
    cout << "ftop_df=" << f << "  第1自由度=" << df1 << "  第2自由度=" << df2 << "  上側p値=" << ftop_df(f, df1, df2) << "\n";
  }
}

void pvalclass::pf() const      // F値←上側p値変換
{
  double pr = 0.0;
  int df1 = 1, df2 = 1;
  while ((pr = get_fdata("上側p値(0:終了)", 0.0, 1.0, pr)) > 0.0)
  {
    df1 = (int)(get_fdata("第1自由度", 1.0, (double)INT_MAX, (double)df1));
    df2 = (int)(get_fdata("第2自由度", 1.0, (double)INT_MAX, (double)df2));
    cout << "上側p値=" << pr << "  第1自由度=" << df1 << "  第2自由度=" << df2 << "  F値=" << ptof((ldbl)pr, df1, df2) << "\n" << endl;
  }
}

void pvalclass::dp() const      // d値→p値変換
{
  double d = 0.0;
  int g = 2;
  int df = 0;
  double r = 0.5;
  while ((d = get_fdata("d値(0:終了)", -FDmax, FDmax, d)) != 0.0)
  {
    g = (int)(get_fdata("群数", 2.0, (double)INT_MAX, (double)g));
    df = (int)(get_fdata("自由度(0:∞)", 0.0, (double)INT_MAX, (double)df));
    r = get_fdata("相関係数", 0.0, 1.0, r);
    double pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
    if (df > 0) cout << "d値=" << d << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  両側p値=" << pr << "\n";
    else cout << "d値=" << d << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  両側p値=" << pr << "\n";
    pr = dntop(d, g, df, r);
    if (df > 0) cout << "dn=" << d << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  両側p値=" << pr << "\n";
    else cout << "dn=" << d << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  両側p値=" << pr << "\n";
    pr = dntop(d, g, df, r, true);
    if (df > 0) cout << "dn=" << d << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  片側p値=" << pr << "\n";
    else cout << "dn=" << d << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  片側p値=" << pr << "\n";
    cout << endl;
  }
}

void pvalclass::pd() const      // d値←両側p値変換
{
  double pr = 0.0;
  int g = 2;
  int df = 0;
  double r = 0.5;
  while ((pr = get_fdata("両側p値(0:終了)", 0.0, 1.0, pr)) > 0.0)
  {
    g = (int)(get_fdata("群数", 2.0, (double)INT_MAX, (double)g));
    df = (int)(get_fdata("自由度(0:∞)", 0.0, (double)INT_MAX, (double)df));
    r = get_fdata("相関係数", 0.0, 1.0, r);
    double d = (double)(ptod((ldbl)pr, g, df, (ldbl)r));
    if (df > 0) cout << "両側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  d値=" << d << "\n";
    else cout << "両側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  d値=" << d << "\n";
    d = ptodn(pr, g, df, r);
    if (df > 0) cout << "両側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  dn=" << d << "\n";
    else cout << "両側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  dn=" << d << "\n";
    d = ptodn(pr, g, df, r, true);
    if (df > 0) cout << "片側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  dn=" << d << "\n";
    else cout << "片側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  相関係数=" << r << "  dn=" << d << "\n";
    cout << endl;
  }
}

void pvalclass::qp() const      // q値→上側p値変換
{
  double q = 0.0;
  int g = 2;
  int df = 0;
  while ((q = get_fdata("q値(0:終了)", -FDmax, FDmax, q)) != 0.0)
  {
    g = (int)(get_fdata("群数", 2.0, (double)INT_MAX, (double)g));
    df = (int)(get_fdata("自由度(0:∞)", 0.0, (double)INT_MAX, (double)df));
    double pr = (double)(qtop(q, g, df));
    if (df > 0) cout << "q値=" << q << "  群数=" << g << "  自由度=" << df << "  上側p値=" << pr << "\n";
    else cout << "q値=" << q << "  群数=" << g << "  自由度=" << df << "  上側p値=" << pr << "\n";
    pr = qntop(q, g, df);
    if (df > 0) cout << "qn=" << q << "  群数=" << g << "  自由度=" << df << "  上側p値=" << pr << "\n";
    else cout << "qn=" << q << "  群数=" << g << "  自由度=" << df << "  上側p値=" << pr << "\n";
    cout << endl;
  }
}

void pvalclass::pq() const      // q値←上側p値変換
{
  double pr = 0.0;
  int g = 2;
  int df = 0;
  while ((pr = get_fdata("上側p値(0:終了)", 0.0, 1.0, pr)) > 0.0)
  {
    g = (int)(get_fdata("群数", 2.0, (double)INT_MAX, (double)g));
    df = (int)(get_fdata("自由度(0:∞)", 0.0, (double)INT_MAX, (double)df));
    double q = (double)(ptoq(pr, g, df));
    if (df > 0) cout << "上側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  q値=" << q << "\n";
    else cout << "上側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  q値=" << q << "\n";
    q = ptoqn(pr, g, df);
    if (df > 0) cout << "上側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  qn=" << q << "\n";
    else cout << "上側p値=" << pr << "  群数=" << g << "  自由度=" << df << "  qn=" << q << "\n";
    cout << endl;
  }
}

void pvalclass::np() const      // 正規偏位z値→累積確率変換
{
  double z = 0.0;
  while ((z = get_fdata("正規偏位z値(0:終了)", -FDmax, FDmax, z)) != 0.0)
  { cout << "正規偏位z値=" << z << "  累積確率=" << nedp((ldbl)z) << "\n" << endl; }
}

void pvalclass::pn() const      // 正規偏位z値←累積確率変換
{
  double pr = 0.0;
  while ((pr = get_fdata("累積確率(0:終了)", 0.0, 1.0, pr)) > 0.0)
  { cout << "累積確率=" << pr << "  正規偏位z値=" << ptoned((ldbl)pr) << "\n" << endl; }
}

void pvalclass::ep() const      // 誤差関数x値→erf(x)変換
{
  double x = 0.0;
  while ((x = get_fdata("誤差関数x値(0:終了)", -FDmax, FDmax, x)) != 0.0)
  {
    cout << "erf(" << x << ")=" << erf(x) << "  erfc(" << x << ")=" << erfc(x) << "\n" << endl;
  }
}

void pvalclass::ncclp() const   // 非心χ2乗値→下側p値変換
{
  double csqr = 0.0, lamda = 0.0;
  int df = 1;
  while ((csqr = get_fdata("χ2乗値(0:終了)", 0.0, FDmax, csqr)) > 0.0)
  {
    df = (int)(get_fdata("自由度", 1.0, (double)INT_MAX, (double)df));
    lamda = get_fdata("非心度", 0.0, FDmax, lamda);
    cout << "χ2乗値=" << csqr << "  自由度=" << df << "  非心度=" << lamda << "  下側p値=" << ncctolp((ldbl)csqr, df, (ldbl)lamda) << "\n" << endl;
  }
}

void pvalclass::lpncc() const   // 非心χ2乗値←下側p値変換
{
  double pr = 0.0, lamda = 0.0;
  int df = 1;
  while ((pr = get_fdata("下側p値(0:終了)", 0.0, 1.0, pr)) > 0.0)
  {
    df = (int)(get_fdata("自由度", 1.0, (double)INT_MAX, (double)df));
    lamda = get_fdata("非心度", 0.0, FDmax, lamda);
    cout << "下側p値=" << pr << "  自由度=" << df << "  非心度=" << lamda << "  χ2乗値=" << lptoncc((ldbl)pr, df, (ldbl)lamda) << "\n" << endl;
  }
}

void pvalclass::nctlp() const   // 非心t値→下側p値変換
{
  double t = 0.0, lamda = 0.0;
  int df = 1;
  while ((t = get_fdata("t値(0:終了)", -FDmax, FDmax, t)) != 0.0)
  {
    df = (int)(get_fdata("自由度", 1.0, (double)INT_MAX, (double)df));
    lamda = get_fdata("非心度", -FDmax, FDmax, lamda);
    cout << "t値=" << t << "  自由度=" << df << "  非心度=" << lamda << "  下側p値=" << ncttolp((ldbl)t, df, (ldbl)lamda) << "\n" << endl;
  }
}

void pvalclass::lpnct() const   // 非心t値←下側p値変換
{
  double pr = 0.0, lamda = 0.0;
  int df = 1;
  while ((pr = get_fdata("下側p値(0:終了)", 0.0, 1.0, pr)) > 0.0)
  {
    df = (int)(get_fdata("自由度", 1.0, (double)INT_MAX, (double)df));
    lamda = get_fdata("非心度", -FDmax, FDmax, lamda);
    cout << "下側p値=" << pr << "  自由度=" << df << "  非心度="<< lamda << "  t値=" << lptonct((ldbl)pr, df, (ldbl)lamda) << "\n" << endl;
  }
}

void pvalclass::ncflp() const   // 非心F値→下側p値変換
{
  double f = 0.0, lamda = 0.0;
  int df1 = 1, df2 = 1;
  while ((f = get_fdata("F値(0:終了)", 0.0, FDmax, f)) > 0.0)
  {
    df1 = (int)(get_fdata("第1自由度", 1.0, (double)INT_MAX, (double)df1));
    df2 = (int)(get_fdata("第2自由度", 1.0, (double)INT_MAX, (double)df2));
    lamda = get_fdata("非心度", 0.0, FDmax, lamda);
    cout << "F値=" << f << "  第1自由度=" << df1 << "  第2自由度=" << df2 << "  非心度=" << lamda << "  下側p値=" << ncftolp((ldbl)f, df1, df2, (ldbl)lamda) << "\n" << endl;
  }
}

void pvalclass::lpncf() const   // 非心F値←下側p値変換
{
  double pr = 0.0, lamda = 0.0;
  int df1 = 1, df2 = 1;
  while ((pr = get_fdata("下側p値(0:終了)", 0.0, 1.0, pr)) > 0.0)
  {
    df1 = (int)(get_fdata("第1自由度", 1.0, (double)INT_MAX, (double)df1));
    df2 = (int)(get_fdata("第2自由度", 1.0, (double)INT_MAX, (double)df2));
    lamda = get_fdata("非心度", 0.0, FDmax, lamda);
    cout << "下側p値=" << pr << "  第1自由度=" << df1 << "  第2自由度=" << df2 << "  非心度=" << lamda << "  F値=" << lptoncf((ldbl)pr, df1, df2, (ldbl)lamda) << "\n" << endl;
  }
}

void pvalclass::test() const    // テストルーチン実行
{
  double d = 0.0, r = 0.5;
  int g = 2, df = 0;
  double pr = dntop(d, g, df, r);
  cout << "\n";
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  pr = dntop(d, g, df, r, true);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  plwr=" << r << "\n";
  g = 3;  df = 20;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  pr = dntop(d, g, df, r, true);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  plow=" << pr<< "\n";
  d = 1.645;  g = 2;  df = 0;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 1.96;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  g = 3;  r = 1.0;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, 0, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 4.303;  df = 2;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 1.916;  g = 3;  df = 0;  r = 0.5;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.212;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.79;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.349;  g = 4;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.442;  g = 5;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.511;  g = 6;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 1.96;  g = 2;  df = 100000;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.212;  g = 3;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.349;  g = 4;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.442;  g = 5;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 4.303;  g = 2;  df = 2;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 5.418;  g = 3;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.379;  g = 3;  df = 20;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.370;  g = 3;  df = 21;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.540;  g = 4;  df = 20;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.470;  g = 4;  df = 31;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.453;  g = 5;  df = 360;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.442;  g = 5;  df = 10000;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.511;  g = 6;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 1.96;  g = 2;  df = 0;  r = 0.7;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 5.238;  g = 3;  df = 2;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.331;  g = 3;  df = 21;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.254;  g = 3;  df = 20;  r = 0.9;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.236;  g = 3;  df = 0;  r = 0.0;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  d = 2.454;  g = 4;  df = 60;  r = 0.0;
  pr = dntop(d, g, df, r);
  cout << "dn=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptodn=" << ptodn(pr, g, df, r) << "\n";
  pr = (double)(dtop((ldbl)d, g, df, (ldbl)r));
  cout << "d=" << d << "  g=" << g << "  df=" << df << "  r=" << r << "  pr=" << pr << " ptod=" << ptod((ldbl)pr, g, df, (ldbl)r) << "\n";
  double q = 0.0;
  g = 2;  df = 0;
  pr = qntop(q, g, df);
  cout << "qn=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoqn=" << ptoqn(pr, g, df) << "\n";
  pr = (double)(qtop((ldbl)q, g, df));
  cout << "q=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoq=" << ptoq((ldbl)pr, g, df) << "\n";
  g = 3;  df = 21;
  pr = qntop(q, g, df);
  cout << "qn=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoqn=" << ptoqn(pr, g, df) << "\n";
  pr = (double)(qtop((ldbl)q, g, df));
  cout << "q=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoq=" << ptoq((ldbl)pr, g, df) << "\n";
  q = 2.772;  g = 2;  df = 0;
  pr = qntop(q, g, df);
  cout << "qn=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoqn=" << ptoqn(pr, g, df) << "\n";
  pr = (double)(qtop((ldbl)q, g, df));
  cout << "q=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoq=" << ptoq((ldbl)pr, g, df) << "\n";
  q = 3.314;  g = 3;
  pr = qntop(q, g, df);
  cout << "qn=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoqn=" << ptoqn(pr, g, df) << "\n";
  pr = (double)(qtop((ldbl)q, g, df));
  cout << "q=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoq=" << ptoq((ldbl)pr, g, df) << "\n";
  q = 3.633;  g = 4;
  pr = qntop(q, g, df);
  cout << "qn=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoqn=" << ptoqn(pr, g, df) << "\n";
  pr = (double)(qtop((ldbl)q, g, df));
  cout << "q=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoq=" << ptoq((ldbl)pr, g, df) << "\n";
  q = 4.387;  g = 9;
  pr = qntop(q, g, df);
  cout << "qn=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoqn=" << ptoqn(pr, g, df) << "\n";
  pr = (double)(qtop((ldbl)q, g, df));
  cout << "q=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoq=" << ptoq((ldbl)pr, g, df) << "\n";
  q = 6.085;  g = 2;  df = 2;
  pr = qntop(q, g, df);
  cout << "qn=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoqn=" << ptoqn(pr, g, df) << "\n";
  pr = (double)(qtop((ldbl)q, g, df));
  cout << "q=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoq=" << ptoq((ldbl)pr, g, df) << "\n";
  q = 3.578;  g = 3;  df = 20;
  pr = qntop(q, g, df);
  cout << "qn=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoqn=" << ptoqn(pr, g, df) << "\n";
  pr = (double)(qtop((ldbl)q, g, df));
  cout << "q=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoq=" << ptoq((ldbl)pr, g, df) << "\n";
  q = 4.588;  g = 9;  df = 49;
  pr = qntop(q, g, df);
  cout << "qn=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoqn=" << ptoqn(pr, g, df) << "\n";
  pr = (double)(qtop((ldbl)q, g, df));
  cout << "q=" << q << "  g=" << g << "  df=" << df << "  pr=" << pr << " ptoq=" << ptoq((ldbl)pr, g, df) << "\n";
}

int main(int argc, char* argv[]) // メインルーチン
{
  pvalclass pval;               // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    int kind = 0;
    while ((kind = pval.sel(kind)) != 0)
    {
      switch (kind)
      {
      case 1 : pval.zp();
        break;
      case 2 : pval.pz();
        break;
      case 3 : pval.cp();
        break;
      case 4 : pval.pc();
        break;
      case 5 : pval.tp();
        break;
      case 6 : pval.pt();
        break;
      case 7 : pval.fp();
        break;
      case 8 : pval.pf();
        break;
      case 9 : pval.dp();
        break;
      case 10 : pval.pd();
        break;
      case 11 : pval.qp();
        break;
      case 12 : pval.pq();
        break;
      case 13 : pval.np();
        break;
      case 14 : pval.pn();
        break;
      case 15 : pval.ep();
        break;
      case 16 : pval.ncclp();
        break;
      case 17 : pval.lpncc();
        break;
      case 18 : pval.nctlp();
        break;
      case 19 : pval.lpnct();
        break;
      case 20 : pval.ncflp();
        break;
      case 21 : pval.lpncf();
        break;
      case 98 : pval.test();
        break;
      case 99 : help(opt.pname(), USAGE);
        break;
      default : disperr(opt.pname(), "指定番号不適切！");
      }
    }
    disperr(opt.pname(), "終了！");
  }
  catch (help_except& helpex)   // ヘルプファイル表示
  { return 0; }
  catch (option_error& opterr)  // オプション指定エラー
  { return 1; }
  catch (runtime_error& rterr)  // 実行時エラー
  {
    disperr(opt.pname(), rterr.what());
    return 2;
  }
  catch (...)                   // 予期せぬ例外
  {
    disperr(opt.pname(), "内容不明なエラー発生！");
    return 255;
  }
  return 0;
}
