(* chiR.m, 2021-01-05 
Example: cdf4[50,100]    r=4, n=50, y=100
Example: Plot[{cdf6[1, y], cdf6[10, y], cdf6[50, y]}, {y, 0, 100}]
   BUG: mathematica の数値積分手法のコントロール方法に不明部分があるため,
      r が一般の場合のコード (cdfEven)と r=4, r=6 専用のコード(cdf4,cdf6)
      で結果が不一致となる.
      xi*[s] の値は同じなのに不思議.
    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 の計算, 一般の実装. BUG: n が小さい時のみ正しく動く. *)
cdfEven[r_, n_, y_] := 
 NIntegrate[((1 - Exp[-y*s])/s)*N[xiEven[r, n, s]], {s, 0, Infinity}]

(* 以下は r=6 の専用のコード. n が100くらいでも動く *)
(* *r6-case2.nb をコピーして答えが違う原因を探る*)
hgvarphi36 = (1/6*Gamma[1/6]*z)*
     HypergeometricPFQ[{1/6}, {1/3, 2/3}, -1/27*wr^3*z^6] + (-1/6*
       Gamma[1/2]*wr*z^3)*
     HypergeometricPFQ[{1/2}, {2/3, 4/3}, -1/27*wr^3*z^6] + (1/12*
       Gamma[5/6]*wr^2*z^5)*
     HypergeometricPFQ[{5/6}, {4/3, 5/3}, -1/27*wr^3*z^6] /. {wr -> 
     Exp[2*Pi*I/6]/2};
hgvarphi3s6[s_] := N[hgvarphi36 /. {z -> s^(-1/6)}]
xi6[n_?NumericQ, s_?NumericQ] := 
 Block[{a, b, theta, f, rr, v0005, v0000, border, upBorder, v8, v9},
  (* s が border より小さいときは線形近似*)
  v0005 = 1.0961216365125503` - 0.6312643298893801` I;
  v0000 = 1.08540188404668` - 0.6266570540070552` I;
  
  border = 0.0005;
  upBorder = 10^8;
  f = Which[s < border, N[v0000*(1 - s/border) + (s/border)*v0005],
    s > upBorder, N[hgvarphi36 /. {z -> N[s^(-1/6)]}],
    True, hgvarphi3s6[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/6 + n*theta]]]]
(* r=6 の時の CDF *)
cdf6[n_,y_]:=NIntegrate[((1 - Exp[-y*s])/s)*N[xi6[n, s]], {s, 0, Infinity}]

(* 以下は r=4 での専用コード *)
hgvarphi34 = (1/4*Gamma[1/4]*z)*
     HypergeometricPFQ[{1/4}, {1/2}, 
      1/4*wr^2*z^4] + (-1/4*Gamma[3/4]*wr*z^3)*
     HypergeometricPFQ[{3/4}, {3/2}, 1/4*wr^2*z^4] /. {wr -> 
     Exp[2*Pi*I/4]/2};
hgvarphi3s4[s_] := N[hgvarphi34 /. {z -> s^(-1/4)}]
xi4[n_?NumericQ, s_?NumericQ] := 
 Block[{a, b, theta, f, rr, v0005, v0000, border, upBorder, v8, v9},
  border = 0.0001;
  f = Which[s < border, N[varphi3s[4, s]],
    True, hgvarphi3s4[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/4 + n*theta]]]]
(* r=4 の時の CDF *)
cdf4[n_,y_]:=NIntegrate[((1 - Exp[-y*s])/s)*N[xi4[n, s]], {s, 0, Infinity}]

