(* chiR.m, 2021-01-05 
Example: cdfEven[4,50,100]    r=4, n=50, y=100
Example: Plot[{cdfEven[6,1, y], cdfEven[6,10, y], cdfEven[6,50, y]}, {y, 0, 100}]
   開発中のnow-how: cdfEven で N[xiEven[r,n,s]] としてはいけない. N があると値が変.
    Ref: chiR.m.orig1
    Ref: 2021-01-04-chiR.nb
*)

hgEvenRz[r_, z_, wr_] := 
  Block[{r1, i, j}, 
   r1 = r/2; (z/r)*
    Sum[(Gamma[(2*i + 1)/r]/i!)*z^(2*i)*(-wr)^i*
      HypergeometricPFQ[{(2*i + 1)/r}, 
       Complement[Table[(i + j)/r1, {j, 1, r1}], {1}], 
       z^(2*r1)*(-wr)^r1*r1^(-r1)], {i, 0, r1 - 1}]];
(* HG による varphi_3(s) *)
hgEvenR[r_, s_] := Block[{z, r1, wr, f},
  f = hgEvenRz[r, z, wr] /. {wr -> Exp[2*Pi*I/r]/2};
  Return[ f /. {z -> s^(-1/r)}]
  ]

hgOddRz[r_, z_, wr_] := 
  Block[{i, j}, (z/r)*
    Sum[(Gamma[(2*i + 1)/r]/i!)*z^(2*i)*(-wr)^i*
      HypergeometricPFQ[{(2*i + 1)/(2*r), (2*i + 1)/(2*r) + 1/2}, 
       Complement[Table[(i + j)/r, {j, 1, r}], {1}], 
       z^(2*r)*(-wr)^r*r^(-r)*4], {i, 0, r - 1}]];
(* HG による varphi_3(s) *)
hgOddR[r_, s_] := Block[{z, wr, f},
  f = hgOddRz[r, z, wr] /. {wr -> Exp[2*Pi*I/r]/2};
  Return[ f /. {z -> s^(-1/r)}]
  ]

(* 数値積分による varphi_3(s) *)
varphi3s[r_, s_] := 
 NIntegrate[Exp[-s*t^r]*Exp[-Exp[2*Pi*I/r]*t^2/2], {t, 0, Infinity}]

asympPhi3z[r_, zeta_, upto_] := 
 Evaluate[(zeta/r)*(r/2)*
   Sum[(-1)^m*Gamma[(m*r + 1)/2]*zeta^(-(m*r + 1)/2)/m!, {m, 0, upto}]]

(* s=0 での漸近展開 *)
asympPhi3[r_, s_, upto_] := 
  asympPhi3z[r, (s^(-1/r))^2*Exp[2*Pi*I/r]/2, 
    upto]/(s^(-1/r)*Exp[2*Pi*I/r]/2);

myArcTan[a_?NumericQ, b_?NumericQ] := Which[a > 0, N[ArcTan[b/a]],
  (a < 0) && (b >= 0), N[ArcTan[b/a] + Pi],
  (a < 0) && (b < 0), N[ArcTan[b/a] - Pi],
  (a == 0) && (b >= 0), N[Pi/2],
  (a == 0) && (b < 0), N[-Pi/2]]

(* 線形近似のために積分をあらかじめ計算して表にしておく *)
ChiRmax = 8; 
(* ChiRBorder=Table[0.0005,{i,1,ChiRmax}]; *)
ChiRBorder = Table[0.0005, {i, 1, ChiRmax}];
ChiRphi3At0 = Join[{0, 0, 0}, Table[varphi3s[r, 0], {r, 4, ChiRmax}]]
(* Todo, r=3 の数値積分は対策必要 *)
ChiRphi3AtB = 
 Join[{0, 0, 0}, Table[varphi3s[r, ChiRBorder[[r]]], {r, 4, ChiRmax}]]

(* s = 0 の nbd で varphi_3(s) を線形近似 *)
phi3SmallS[r_, s_] := 
 If [r > ChiRmax, N[asympPhi3[r, s, 1]], 
  N[ChiRphi3At0[[r]]*(1 - s/ChiRBorder[[r]]) + (s/ChiRBorder[[r]])*
     ChiRphi3AtB[[r]]]]

(* r even の場合の xi(s) *)
xiEven[r_, n_?NumericQ, s_?NumericQ] := Block[{a, b, theta, f, rr},
  If[Mod[r, 2] == 1, ErrorBox["r odd case is not implemented"]];
  f = Which[s < ChiRBorder[[r]], phi3SmallS[r, s], (*数値積分から補完*)
    s > 10^8, N[hgEvenRz[r, z] /. {z -> N[s^(-1/r)]}],
    True, hgEvenR[r, s]];
  a = Re[f]; b = Im[f]; theta = myArcTan[a, b]; 
  rr = N[Sqrt[a^2 + b^2]];
  Return[N[(1/Pi)*(2/Pi)^(n/2)]*rr^n*N[Sin[n*Pi/r + n*theta]]]]

(* r even の場合の CDF の計算, 一般の実装.  *)
cdfEven[r_, n_, y_] := 
 NIntegrate[((1 - Exp[-y*s])/s)*xiEven[r, n, s], {s, 0, Infinity}]

