BeginPackage["gammaReduce`",{"jTool`","loopTool`","multipleZeta`","laporta2`"}];

greX::usage="";

greL::usage="";

greδ::usage="";

greϵ::usage="";

greSC::usage="";

grePermId::usage="";

greMakePIDFunc::usage="";

greSimplify::usage="";

greRecurrenceId::usage="";

greMakeRIDFunc::usage="";

greSumToL::usage="";

greLToSum::usage="";

greIntToL::usage="";

greLToInt::usage="";

greGPGSumToDL::usage="";

greD::usage="";

greToD::usage="";

greToSC::usage="";

greLToSCInt::usage="";

greSplitL::usage="";

greConvergenceQ::usage="";

greEvaluateSCL::usage="";

grePermId::usage="";

greMakePIDFunc::usage="";

Begin["`private`"]

Unprotect["gammaReduce`*"];

Unprotect["gammaReduce`private`*"];

[約分]

Attributes[greX]={NHoldAll};

greL[a___,{{b___,c_,d___},{e___,c_,f___},x___},g___]:=greL[a,{{b,d},{e,f},x},g];

Derivative[a___,{{b___,0,d___},{e___,0,f___},x___},g___][greL][p__]=.//Quiet;

Derivative[a___,{{b___,0,d___},{e___,0,f___},x___},g___][greL][p__]:=
	Module[{tmp,pos},
		pos=Position[{a,{{b,tmp,d},{e,tmp,f},x},g},tmp];
		Derivative[a,{{b,d},{e,f},x},g][greL]@@ReplacePart[{p},pos->Sequence[]]/;
			Quiet[SameQ@@Extract[{p},pos]]
	];

[PermID等の簡単化]

grePermId[n_Integer?NonNegative,d:{_Integer?NonNegative..},e:{_Integer?NonNegative..}]/;Length[d]==Length[e]==(2^n):=
Module[{K,tmp,i,j,k,ident},
tmp=Table[greX[i-1,j+k-1],{i,(2^n)},{j,2},{k,1,2{d,e}[[j,i]],2}];
tmp=K@@MapThread[Append,{tmp,Array[ident[greX[#,0]]&,(2^n),0]}];
tmp=tmp[[#]]&/@Ordering[IntegerDigits[Range[0,(2^n)-1],2,n][[1;;,#]]]&/@Permutations[Range[n]];
greL@@List@@@Partition[#,3]&/@Tuples[Permutations/@Sequence@@@#]&/@tmp/.ident->Identity//Flatten
];

grePermId[n_Integer?NonNegative,d:{_Integer?NonNegative..}]/;Length[d]==(2^n):=grePermId[n,d,d];

greMakePIDFunc[in_,f_]:=
Module[{tmp,var,len,rep,lhs,rhs},
SetDelayed@@tmp/;
Quiet@Check[
tmp=in/.greX[i_,j_]:>ToExpression["Hold[greX"<>ToString[i]<>"s"<>ToString[j]<>"]"];
len=Log2@Length@in[[1]];
rep=Ordering/@List@@@in[[1;;,2^Range@Log2@Length@in[[1]]/2+1,3,1]];
rep=Thread[Array[greX,len]->greX/@#]&/@rep;
rhs=MapThread[Hold[ReplaceAll[#1,#2]]&,{tmp,rep}];
lhs=Replace[Replace[tmp[[1]],Hold[a_]:>a_,{3}],Hold[a_]:>a___,{2}]//Quiet;
tmp=Block[{f},Replace[Hold[Evaluate[f@@lhs,rhs]],Hold[a_]:>a,{2,∞}]];
True,
False
]
];

(*[Derivativeなどの簡単化]*)

(* Derivative[__][greL][__]同士の関係を求める。例えば、
Derivative[{{},{}},{{0},{2}}][greL][{{},{}},{{1},{1}}]-
2*Derivative[{{},{}},{{1},{1}}][greL][{{},{}},{{1},{1}}]-
Derivative[{{},{}},{{2},{0}}][greL][{{},{}},{{1},{1}}]==0
のように、ある微分係数は、同じ引数 (この場合1)を持つ部分の微分の次数の合計 (この場合2)を変えずに、
組み合わせを変えたものの線形結合で表される場合がある。derrel[分子の個数,分母の個数,次数の合計]でこの関係を求める。
上の例に対応する関係式はderrel[1,1,2]で求められる。*)

$derrel//Clear;
$derrel[n_,m_,k_]:=
	Module[{x,y,z,sumrep,cdiff,patt,tmp,cc},
sumrep=greLToSum[greL[{{},{},{}},{Array[x,n],Array[x,m,n+1],{}}]];
cdiff=With[{a=z[sumrep,Sequence@@Array[{x[#],y[#]}&,n+m]]},a/._x->0&]/.z->D/.y->Slot;
(* cdiff --> D[sumrep,{x[1],#1},{x[2],#2},...,{x[n+m],# (n+m)}]& *)
patt=Join@@(Permutations/@IntegerPartitions[k,{n+m},Range[0,k]]);
tmp=Expand[(cdiff@@##cc@@##&/@patt//Total)/cdiff@@ConstantArray[0,m+n]];
tmp=DeleteCases[Flatten[CoefficientList[tmp,Union[Cases[tmp,_PolyGamma,{0,∞}]]]],0]//Union;
processLaportaFirst[tmp,cc,"sortFunc"->Reverse];
tmp=processLaportaThird[getVariable[tmp,cc]];
tmp=Expand[ss@@##cc@@##&/@patt/.Rule@@@tmp//Total];
tmp=tmp/.#->1/._cc->0&/@getVariable[tmp,cc];
processLaportaFirst[tmp,ss,"sortFunc"->({-Min[Count[#[[1;;n]],0],Count[#[[n+1;;]],0]],Reverse[#]}&)];
processLaportaThird[getVariable[tmp,ss]]
];

(* $drrrelで求めた関係を用い、微分係数の数を減らす。ここで、上の例は
Derivative[{{},{}},{{0,0},{2,1}}][greL][{{},{}},{{1,2},{1,3}}]-
2*Derivative[{{},{}},{{1,0},{1,1}}][greL][{{},{}},{{1,2},{1,3}}]-
Derivative[{{},{}},{{2,0},{0,1}}][greL][{{},{}},{{1,2},{1,3}}]==0
のような場合にも適用できることに注意。*)

dersimp[in_]:=
in//.Derivative[a___][greL][x___,{s_,t_,u_},y___]/;Intersection[s,t]=!={}:>
Module[{intt,fres,npos,dpos,derv,dord,rsid,argv,argr,rpos,resr},
intt=Intersection[s,t];
fres={};
Do[
npos=Position[s,f]//Flatten;
dpos=Position[t,f]//Flatten;
derv={a}[[Length[{x}]+1]];
dord=Total[derv[[1,npos]]]+Total[derv[[2,dpos]]];
rsid=$derrel[Length[npos],Length[dpos],dord];
argv=Flatten[{derv[[1,npos]],derv[[2,dpos]]}];
argr=ss@@argv/.Rule@@@rsid;
If[ss@@argv===argr,Continue[]];
rpos=Join[{1,#}&/@npos,{2,#}&/@dpos];
resr=argr/.ss[v__]:>ReplacePart[ss@@derv,Thread[rpos->{v}]];
fres=resr/.ss[v__]:>Derivative[Sequence@@ReplacePart[{a},Length[{x}]+1->{v}]][greL][x,{s,t,u},y];
If[fres=!={},Break[]],
{f,intt}
];
fres/;fres=!={}
]//Expand;

(* permutation ID に相当する簡約化 *)
applyPermId[greL[b__],g_]:=
Module[{tmp,pos,tmp2},
tmp=g[b];
If[tmp[[0]]===g,
greMakePIDFunc[grePermId[Log[2,Length[{b}]],
Length/@{b}[[1;;,1]],Length/@{b}[[1;;,2]]],g];
tmp=g[b];
];
greL[Sequence@@SortBy[tmp,List@@#/.Thread[Variables[List@@#]->0]&][[1]]]
];

applyPermId2[greL[b__],g_]:=
Module[{tmp,pos,tmp2},
tmp=g[b];
If[tmp[[0]]===g,
greMakePIDFunc[grePermId[Log[2,Length[{b}]],
Length/@{b}[[1;;,1]],Length/@{b}[[1;;,2]]],g];
tmp=g[b];
];
greL[Sequence@@#]&/@SortBy[tmp,List@@#/.Thread[Variables[List@@#]->0]&]
];

applyPermId[Derivative[a__][greL][b__],g_]:=
Module[{tmp,pos,tmp2},
tmp=Block[{greL=List},g[b]];
If[tmp[[0]]===g,
greMakePIDFunc[grePermId[Log[2,Length[{b}]],
Length/@{b}[[1;;,1]],Length/@{b}[[1;;,2]]],g];
tmp=Block[{greL=List},g[b]];
];
pos=Position[tmp,SortBy[tmp,List@@#/.Thread[Variables[List@@#]->0]&][[1]]]//Flatten;
tmp2=Block[{greL=List},g@@Join[{a},{b}[[Length[{a}]+1;;]]]];
Derivative[Sequence@@Sort[tmp2][[1]]][greL][Sequence@@tmp[[pos[[1]]]]]
];

applyPermId2[Derivative[a__][greL][b__],g_]:=
Module[{tmp,pos,tmp2},
tmp=Block[{greL=List},g[b]];
If[tmp[[0]]===g,
greMakePIDFunc[grePermId[Log[2,Length[{a}]],
Length/@{a}[[1;;,1]],Length/@{a}[[1;;,2]]],g];
tmp=Block[{greL=List},g[b]];
];
tmp2=Block[{greL=List},g@@Join[{a},{b}[[Length[{a}]+1;;]]]];
tmp=SortBy[Transpose[{tmp,tmp2}],List@@#/.Thread[Variables[List@@#]->0]&];
Derivative[Sequence@@#[[2]]][greL][Sequence@@#[[1]]]&/@tmp
];

applyPermIdAll[in_,g_]:=in//.a:(Derivative[__][greL][__]|_greL):>applyPermId[a,g];

greSimplify[in_,g_]:=FixedPoint[dersimp[applyPermIdAll[#,g]]&,in]//Expand;

[漸化式]

id1[n_Integer?Positive,d:{_Integer?NonNegative..}]/;Length[d]==(2^n):=
Module[{in,tmp,id1,id2,id3,id4,id5,id6,p,m,v,dummy,K},
in=K@@Table[greX[i-1,j+k-1],{i,(2^n)},{j,2},{k,1,2d[[i]],2}];tmp=Function[f,Select[Range[0,2^n-1],Reverse[IntegerDigits[#,2,n]][[f]]==1&]]/@Range[n]+1;
(
id6=Times@@(p/@Flatten[List@@in[[#]]])-1+v@@#&/@tmp;
tmp=Select[tmp[[1;;]],Length[in[[#[[1]],1]]]>0&];
id5=Union@@Function[v,
Times@@m/@DeleteCases[Flatten[List@@in[[v,1]]],#[[1]]](1-(#[[1]]-1)m[#[[1]]])-
Times@@p/@DeleteCases[Flatten[List@@in[[v,2]]],#[[2]]](1-(#[[2]]-1)p[#[[2]]])
&/@Tuples[in[[v[[1]]]]]
]/@tmp;
id1=p[#[[1]]]-m[#[[2]]]-(#[[1]]-#[[2]]+1)&/@Tuples[#]&/@List@@in//Flatten;
id2=p[#[[1]]]-p[#[[2]]]-(#[[1]]-#[[2]])&/@Subsets[#,{2}]&/@List@@in[[1;;,1]]//Flatten;
id3=m[#[[1]]]-m[#[[2]]]-(#[[1]]-#[[2]])&/@Subsets[#,{2}]&/@List@@in[[1;;,2]]//Flatten;
tmp=Union@@Function[f,Prepend[#,f]&/@Select[IntegerPartitions[f,{2,∞}],Sort[#]==Union[#]&]]/@Range[(2^n)-1]+1;
id4=-p[#[[1]]]+Total[p/@#[[2;;]]]+#[[1]]-Total[#[[2;;]]]&/@Tuples[List@@in[[#,1]]]&/@tmp//Flatten;
tmp=Union[id1,id2,id3,id4,id5,id6]in//Expand;
p[j_]K[s__]^:=K@@({s}/.j->j+1);
m[j_]K[s__]^:=K@@({s}/.j->j-1);
v[j__]K[s__]^:=
Module[{arg,pos},
arg={s};
pos={j};
arg[[1]]=Join@@@Transpose[{arg[[1]],arg[[pos[[1]]]]}];
Do[
arg[[{pos[[i]],pos[[i]]-pos[[1]]+1}]]={Join@@@Transpose[arg[[{pos[[i]],pos[[i]]-pos[[1]]+1}]]],dummy},
{i,2,Length[pos]}
];
K@@Delete[arg,pos[[1]]]/.dummy->Sequence[]
];
Collect[tmp/.K->greL,greL[__],Simplify]
)/;And@@(Union[in[[#]]]=!=K[{{},{}}]&/@tmp)
	];

srule={
Hypergeometric2F1[1,b_,c_,z_]:>-(((b-c) Hypergeometric2F1[1,b,1+c,z])/c)+(b Hypergeometric2F1[1,1+b,1+c,z])/c,
Hypergeometric2F1[1,b_,c_,z_]:>(1-z) Hypergeometric2F1[1,1+b,c,z]-((1-c) z Hypergeometric2F1[1,1+b,1+c,z])/c,
Hypergeometric2F1[1,a1_,a2_,z_]:>(1-a2)/(1+a1-a2)+(a1 Hypergeometric2F1[1,1+a1,a2,z])/(1+a1-a2)-(a1 z Hypergeometric2F1[1,1+a1,a2,z])/(1+a1-a2)
};

id2[n_Integer?Positive,d:{_Integer?NonNegative..}]/;Length[d]==(2^n):=
Module[{in,tmp,p1,p2,p3,p4,greY},
in=greL@@Table[greY[i-1,j+k-1],{i,2^n},{j,2},{k,1,2d[[i]],2}];
tmp=Function[f,Select[Range[0,2^n-1],Reverse[IntegerDigits[#,2,n]][[f]]==1&]]/@Range[n]+1;
(
p1=in//greLToInt;
p2=Position[p1,_Hypergeometric2F1];
p3=Tuples[{p2,Range[Length[srule]]}];
p4=Function[v,MapAt[#/.srule[[v[[2]]]]&,p1,v[[1]]]//Expand//greIntToL]/@p3;
Select[Numerator[Together[in - p4]],FreeQ[#,greX]&]/.greY->greX
)/;And@@(Union[in[[#]]]=!=K[{{},{}}]&/@tmp)
  ];

greRecurrenceId[n_Integer?Positive,d:{_Integer?NonNegative..}]/;Length[d]==(2^n):=
Module[{in,tmp,K},
in=K@@Table[greX[i-1,j+k-1],{i,2^n},{j,2},{k,1,2d[[i]],2}];tmp=Function[f,Select[Range[0,(2^n)-1],Reverse[IntegerDigits[#,2,n]][[f]]==1&]]/@Range[n]+1;{in/.K->greL,Collect[Join[id1[n,d],id2[n,d]],greL[__],Simplify]}/;
And@@(Union[in[[#]]]=!=K[{{},{}}]&/@tmp)
];

greMakeRIDFunc[in_,f_]:=
Module[{tmp},
tmp=in/.greX[i_,j_]:>ToExpression["Hold[greX"<>ToString[i]<>"s"<>ToString[j]<>"]"];
tmp=Block[{f},
Join[
Replace[Hold[Evaluate[f@@tmp[[1]]]],Hold[a_]:>a_,{-2}]//Quiet,
Replace[Hold[Evaluate[tmp[[2]]]],Hold[a_]:>a,{-2}]
]
];
SetDelayed@@tmp;
];

[表示の変換]

fff[f_,n_,comb_]:=
Reverse[
Fold[
Function[{u,v},
{Sequence@@u[[1;;-2]],Select[u[[-1]],FreeQ[#[[Sequence@@n]]-v,greX]&],Select[u[[-1]],!FreeQ[#[[Sequence@@n]]-v,greX]&]}],
{f},Reverse[Rest[comb]]
]
];

greSumToL//Clear;
greSumToL[in_]:=
Module[{data,var,nvar,comb,i,tmp,numGamma,denGamma,numPoly,denPoly,numzsum,numExp,nt1,nt2,coeffs,res},
data=in//Expand;
data={Replace[data,Plus->List,{1},Heads->True]}//Flatten;
data=Flatten[{Replace[#,Times->List,{1},Heads->True]}]&/@data;
var=Cases[in,_greX,{0,∞}]//Union;
nvar=Array[greX,Length[var]];
data=data/.Thread[var->nvar];
comb=Total/@SortBy[Subsets[Array[greX,Length[nvar]]],Total[#/.greX[a_]:>(a^2)]&];

res/;Catch[
res=Plus@@Table[
tmp=data[[i]];

numGamma=Cases[tmp,((x_Gamma)^(a_.))/;IntegerQ[a]&&a>0];
denGamma=Cases[tmp,((x_Gamma)^(a_.))/;IntegerQ[a]&&a<0];
numPoly=Cases[tmp,(((x:_Plus|_greX))^(a_.))/;IntegerQ[a]&&a>0];
denPoly=Cases[tmp,(((x:_Plus|_greX))^(a_.))/;IntegerQ[a]&&a<0];
(*numzsum=Cases[tmp,((x_hplZC)^(a_.))/;IntegerQ[a]&&a>0&&MatchQ[x[[4]],{0..}]];*)
numzsum=Cases[tmp,(((x:hplZC[{_},∞,{_},{_}]))^(a_.))/;IntegerQ[a]&&a>0];
numExp=Cases[tmp,(((b_))^(a_))/;!FreeQ[{a},greX]];
coeffs=Times@@Cases[tmp,a_/;FreeQ[{numGamma,denGamma,numPoly,denPoly,numzsum,numExp},a,{2}]];
numExp=Times@@numExp;
nt1=PowerExpand[(numExp/.#->1)/(numExp/.#->0)]&/@nvar;
nt2=numExp/Times@@(nt1^nvar)//PowerExpand;
nt2=nt2/.((x_)^(y_)):>(x^(Expand[y]));

If[!FreeQ[{coeffs,nt2},greX],Throw[False]];
coeffs*=Simplify[nt2];
numExp=nt1;

{numGamma,denGamma,numPoly,denPoly,numzsum}=
Flatten[Cases[#,((a_)^(b_.)):>ConstantArray[a,Abs[b]]]]&/@{numGamma,denGamma,numPoly,denPoly,numzsum};
numGamma=fff[numGamma,{1},comb];
denGamma=fff[denGamma,{1},comb];
numPoly=If[Append[DeleteCases[Coefficient[#,nvar],0],1][[1]]==-1,coeffs*=-1;-#,#]&/@numPoly;
denPoly=If[Append[DeleteCases[Coefficient[#,nvar],0],1][[1]]==-1,coeffs*=-1;-#,#]&/@denPoly;
numPoly=fff[numPoly,{All},comb];
denPoly=fff[denPoly,{All},comb];
(*numzsum=fff[numzsum,{2},comb];*)
numzsum=fff[numzsum,{4,1},comb];

numGamma=MapThread[Function[{u,v},#[[1]]-v&/@u],{numGamma,comb}];
denGamma=MapThread[Function[{u,v},#[[1]]-v&/@u],{denGamma,comb}];
numPoly=MapThread[Function[{u,v},#-v&/@u],{numPoly,comb}];
denPoly=MapThread[Function[{u,v},#-v&/@u],{denPoly,comb}];
(*numzsum=MapThread[Function[{u,v},{#[[2]]-v,#[[1]],#[[3]]}&/@u],{numzsum,comb}];*)
numzsum=MapThread[Function[{u,v},{#[[4,1]]-v,#[[1,1]],#[[3,1]]}&/@u],{numzsum,comb}];
numExp=Function[u,If[MatchQ[u,greX[_]],numExp[[u[[1]]]],1]]/@comb;

coeffs greL@@MapThread[
Function[{p,q,r,s,t,u},{Join[p,r+1,s],Join[q,r,s+1],t,u}],
{numGamma,denGamma,numPoly,denPoly,numzsum,numExp}
]
,{i,Length[data]}
];
True
]
];

greLToSum//Clear;
greLToSum[in_]:=
Module[{lvar,svar,len,xvar,tmp1,tmp2,tmp3},
lvar=Cases[in,_greL,{0,∞}]//Union;
len=Log2[Length/@lvar];
(
len=len[[1]];
xvar=Total/@SortBy[Subsets[Array[greX,len]],Total[#/.greX[a_]:>(a^2)]&];
svar=Map[
Function[v,
tmp1=Times@@(Gamma/@Flatten[v[[1;;,1]]+xvar])/Times@@(Gamma/@Flatten[v[[1;;,2]]+xvar]);
(*tmp2=MapThread[Function[{a,b},hplZC[#[[2]],#[[1]]+b,#[[3]],0#[[2]]]&/@a],{v[[1;;,3]],xvar}];*)tmp2=MapThread[Function[{a,b},hplZC[{#[[2]]},∞,{#[[3]]},{#[[1]]+b}]&/@a],{v[[1;;,3]],xvar}];
tmp3=v[[1;;,4]]^xvar//PowerExpand;
tmp1 Times@@Flatten[tmp2]Times@@tmp3
],
List@@@lvar
];
in/.Thread[lvar->svar]
)/;Length[Union[len]]==1
];

Options[greLToInt]={hplI->False};

greLToInt//Clear;
greLToInt[in_,OptionsPattern[]]:=
Module[{lvar,svar,len,len2,rlen,t1,t2,i,F,j,tmp1,var,tmp2,tmp3,comb,tmpz,xx,xren,ir,p1,p2},
lvar=Cases[in,_greL,{0,∞}]//Union;
lvar=Select[lvar,And@@Equal@@@Map[Length,(List@@#)[[2;;,1;;2]],{2}]&&FreeQ[#,greX]&];
If[OptionValue[hplI]=!=True,lvar=Select[lvar,FreeQ[#,greL[___,{_,_,Except[{}],_},___]]&];];
svar=Map[
Function[v,

len=Length[v];
len2=Ceiling[Log2[len]];
xren=Array[xx,len2];

comb=Total/@SortBy[Subsets[xren],Total[#/.xx[a_]:>(a^2)]&];
tmpz=Table[
(*hplZC@@@Thread[{v[[i,3,1;;,2]],comb[[i]]+v[[i,3,1;;,1]],v[[i,3,1;;,3]],0v[[i,3,1;;,2]]}],*)
hplZC[{#1},∞,{#2},{#3}]&@@@Thread[{v[[i,3,1;;,2]],v[[i,3,1;;,3]],comb[[i]]+v[[i,3,1;;,1]]}],
{i,Length[v]}
]//Flatten;
tmpz={Replace[Times@@tmpz,Times->List,{1},Heads->True]}//Flatten;

tmpz=Times@@Table[hplZ2I[tmpz[[i]],greδ,greX[-i,#]&,"shuffle"->hplZ],{i,Length[tmpz]}];
tmpz=tmpz/.a_ hplI[b_,c___]/;FreeQ[a,hplI]:>hplI[a b,c];

tmpz=Expand[tmpz hplI[1]]//.hplI[a_,b___]hplI[c_,d___]:>hplI[a c,b,d];

Catch[
If[!FreeQ[tmpz,hplZ2I],Throw[v]];
CheckAbort[tmp1=gxSort[Rest[List@@v]],Throw[v]];
tmp1=Times@@Times@@@Table[
(1-greX[i,j])^(-1-Subtract@@tmp1[[i,1;;,j]]+p1)greX[i,j]^(-1+tmp1[[i,1,j]]+p2)/
Gamma[-Subtract@@tmp1[[i,1;;,j]]],{i,Length[tmp1]},{j,Length[tmp1[[i,1]]]}
];
var=Cases[tmp1,_greX,{0,∞}]//Union;
tmp2=Times@@@Table[Select[var,IntegerDigits[#[[1]],2,len2][[i]]==1&],{i,len2}];
rlen=Range[len];
tmp2*=Times@@@Table[v[[Select[rlen,IntegerDigits[#,2,len2][[i]]==1&]+1,4]],{i,len2}];
tmp2*=xren;
tmp1=tmp1/Times@@(1-tmp2) Times@@Gamma[v[[1,1]]]/Times@@Gamma[v[[1,2]]]/.p1|p2->0;
ir=Sequence@@({#,0,1}&/@var);
If[OptionValue[hplI]===True,
Quiet[Check[tmp1=tmpz/.hplI[a_ ,b___]:>hplI[(tmp1/.Thread[xren->xbase[a,xren]])(a/._xx->0),b,ir],Throw[v],Power::infy]];
tmp1=tmp1/.hplI[a_Times,c___]:>Select[a,FreeQ[#,greX]&]hplI[Select[a,!FreeQ[#,greX]&],c],
(*else*)
tmp1/._xx->1
]
]

]
,
lvar
];
in/.Thread[lvar->svar]
];

(* gxSort[ List[{{(Subscript[a, 1,1]),(Subscript[a, 1,2]),...},{(Subscript[α, 1,1]),(Subscript[α, 1,2]),...},(Subscript[v, 1])}, {{(Subscript[a, 2,1]),(Subscript[a, 2,2]),...},{(Subscript[α, 2,1]),(Subscript[α, 2,2]),...},(Subscript[v, 2])}, ...] ] は、
各 {(Subscript[α, i,j])} の並べ替えを行い、
(i) (Subscript[a, i,j])(Subscript[およびα, i,j])-(Subscript[a, i,j])が0または負の整数ではなく、
(ii) (Subscript[α, i,j])-(Subscript[a, i,j])が変数を含む場合、(Subscript[α, i,j])が2以上であり、
(iii) (Subscript[greϵを含んでいないα, i,j]),(Subscript[a, i,j])で、(Subscript[全ての変数を0としたとき整数とならないα, i,j])-(Subscript[a, i,j])ものをできるだけ少なく、
(iv) (Subscript[α, i,j])-(Subscript[a, i,j])が整数となるものが最も多く、
(v) 変数を0としたときに整数となる(Subscript[α, i,j])-(Subscript[a, i,j])が最も多い
組み合わせを返す。
ただし、 (i),(ii)を満たす組み合わせが存在しない場合、失敗する。
頭部は任意であり、また (Subscript[v, i]) は無視される。
また、(Subscript[a, i,j])に0または負の整数がないことを確認し、あった場合は失敗する。*)
gxSort[in_]:=
Module[{F,i,t1,t2,t3,r1,r2,r3},
If[MemberQ[Flatten[in[[1;;,1]]],_?NonPositive],
Print[in];Print["gxSortError"];
Abort[]
];
Map[
Function[ttt,
t3=Reap[
t1=ttt[[1]]/.Thread[Variables[ttt[[1]]]->0];
Do[
t2=i/.Thread[Variables[i]->0];
If[MemberQ[i-ttt[[1]],_Integer?NonPositive],Continue[]];
(*If[MemberQ[t2-t1,_Integer?Negative],Continue[]];*)
If[MemberQ[Transpose[{i,ttt[[1]]}],Variables[Subtract@@#]=!={}&&((#[[1]]>2)/.Alternatives@@Variables[#]->0)&],Continue[]];
r1=Count[i-ttt[[1]],_?(FreeQ[#,greϵ]&&!IntegerQ[#/.Alternatives@@Variables[#]->0]&)];
r2=-Count[i-ttt[[1]],_Integer];
r3=-Count[t2-t1,_Integer];
Sow[F[{r1,r2,r3},ttt[[1]],i]],
{i,Permutations[ttt[[2]]]}
];
][[2]]//Flatten//Sort;
If[t3==={},
Print[in];Print["gxSortError"];
Abort[]
];
t3//First//Rest
],
in
]/.F->List
];

(* xbase[in,list] は、式 in に含まれる (a^(b+...)) の形の因子を探し、a を返す。
ただし b は list の各要素であり、対応する a をリストにして返す *)
Clear[xbase];
xbase[a_,c_List]:=
Module[{pat,fac,base,rem},
pat=Alternatives@@c;
fac=Times@@Select[
If[a[[0]]===Times,List@@a,{a}],
MatchQ[#,((x_)^(y_))/;FreeQ[x,pat]&&!FreeQ[y,pat]]&
];
base=(fac/.#->1)/(fac/.#->0)&/@c//PowerExpand;
rem=fac/Times@@(base^c)//PowerExpand;
rem=rem/.((x_)^(y_)):>(x^(Expand[y]));
base/;FreeQ[rem,pat]
];

yln[s_,n_Integer?NonNegative]:=
Module[{a,δ,t1},
t1=Expand[((n!SeriesCoefficient[Gamma[a+δ],{δ,0,n}])/(Gamma[a]))]/.
((PolyGamma[b_,a])^(c_.)):>F[{1,b+1,c,a}]//.
F[x__]F[y__]:>F[x,y]/.F[x___]:>F@@Sort[{x}];
{t1,F[{0,n,1,a}]}/.a->s
];
yld[s_,n_Integer?NonNegative]:=
Module[{a,δ,t1},
t1=Expand[((n!SeriesCoefficient[1/Gamma[a+δ],{δ,0,n}])/(1/Gamma[a]))]/.
((PolyGamma[b_,a])^(c_.)):>F[{1,b+1,c,a}]//.
F[x__]F[y__]:>F[x,y]/.F[x___]:>F@@Sort[{x}];
{t1,F[{0,n,-1,a}]}/.a->s
];

Clear[ylsort];
ylsort[in_,exps_]:=
Module[{tmp},
If[in[[1,1]]==1,
Return[{1,0,Length[in],Max[in[[1;;,2]]],Plus@@Times@@@in[[1;;,1;;2]],in}]
];
(*{0,Length[in],Count[in[[1;;,3]],-1],Total[in[[1;;,2]]],in},*)
tmp=GatherBy[List@@in,#[[4]]&];
tmp=exps[[#[[1,4]]]]-Total[#[[1;;,3]]]&/@tmp//Abs//Total;
{0,tmp,Length[in],Count[in[[1;;,3]],-1],Total[in[[1;;,2]]],in}
];

greGPGSumToDL//Clear;
Options[greGPGSumToDL]={"output"->greL};
greGPGSumToDL[in_,method:OptionsPattern[]]:=
Module[{gammaPat,tmp,vars,exps,replRule,dexp,texp,tmp2,pp,pm,d,
data,lres,curpart,firstrun,δ,fff,res},
res/;Catch[
If[!MemberQ[{SeriesCoefficient,greL,"output"},"output"/.{method}],Throw[False]];
gammaPat=Cases[in,a:Times[((_Gamma)^(_.)),___]:>Select[a,MatchQ[#,((_Gamma)^(_.))]&],{0,∞}]//Union;
If[gammaPat==={}&&FreeQ[in,_Gamma],gammaPat={1}];
If[Length[gammaPat]≠1,Throw[False]];
gammaPat=gammaPat//First;
data=in/gammaPat//Expand;
If[!FreeQ[data/._PolyGamma->1,greX],Throw[False]];
vars=Union[Cases[data,Gamma[a_]|PolyGamma[___,a_]/;!NumericQ[a]:>a,{0,∞}]];
exps=Exponent[gammaPat,Gamma[#]]&/@vars;
replRule=Thread[vars->Range[Length[vars]]];
data=data/.{
Gamma[a_]:>Gamma[Replace[a,replRule]],
((PolyGamma[a_,b_])^(c_.))/;MemberQ[vars,b]:>F[{1,a+1,c,Replace[b,replRule]}],
PolyGamma[a_,b_?NumericQ]:>FunctionExpand[PolyGamma[a,b]]
}//Expand;
data=data//.F[a__]F[b__]:>F[a,b]/.F[a___]:>F@@Sort[{a}];
dexp=Append[Cases[data,F[a___]:>Plus@@Times@@@{a}[[1;;,2;;3]],{0,∞}],0]//Max;
pp=Table[yln[d,i],{i,dexp}];
pm=Table[yld[d,i],{i,dexp}];
texp=If[#==0,{1,-1},ConstantArray[Sign[#],Abs[#]]]&/@exps;
firstrun=True;
curpart={};
While[
True,
(*do*)
tmp=Join@@Permutations/@Join@@(IntegerPartitions[#,{Length[texp]},Range[0,dexp]]&/@Range[dexp]);
tmp=Union@@Table[
tmp2=tmp[[i]];
tmp2=MapThread[Join@@Permutations/@IntegerPartitions[#1,{Length[#2]},Range[0,#1]]&,{tmp2,texp}]//Tuples;
tmp2=Table[fff@@@Thread[{k,tmp2[[j,k]],texp[[k]]}],{j,Length[tmp2]},{k,Length[texp]}];
Times@@Flatten[#]&/@tmp2/.fff[_,0,_]->1//Union,
{i,Length[tmp]}
];
tmp=Complement[tmp,curpart];
If[tmp==={},Break[]];
curpart=Join[curpart,tmp];
tmp=tmp/.{
fff[a_,b_,1]:>(pp[[b]]/.d->a),
fff[a_,b_,-1]:>(pm[[b]]/.d->a)
}//Expand;
tmp=tmp//.{((F[x__])^(y_)):>F@@Join@@ConstantArray[{x},y],F[x__]F[y__]:>F[x,y]};
tmp=tmp/.F[a___]:>F@@SortBy[{a},#[[{1,2,4,3}]]&]//.
F[a___,{1,b_,c_,d_},{1,b_,e_,d_},f___]:>F[a,{1,b,c+e,d},f]/.F[a___]:>F@@Sort[{a}];
If[firstrun,processLaportaFirst[Equal@@@tmp,F,"sortFunc"->(ylsort[#,exps]&)],processLaportaSecond[Equal@@@tmp]];
lres=processLaportaThird[Cases[{tmp,data},_F,{0,∞}]//Union];
(*Print[lres/.F[a__]:>Times@@(F/@{a})/.{
F[{1,a_,b_,c_}]:>((PolyGamma[a-1,vars[[c]]])^b),
F[{0,a_,b_,c_}]:>a!greSC[((Gamma[vars[[c]]+greδ])^b)((Gamma[vars[[c]]])^(-b)),{greδ,0,a}]
}];*)
data=data/.Dispatch[Rule@@@lres]//Expand;
tmp=Cases[data,a:F[{1,___},___]:>List@@a[[1;;,-1]],{0,∞}]//Flatten//Union;
If[tmp==={},Break[]];
(texp[[#]]=Join[texp[[#]],{1,-1}])&/@tmp;
firstrun=False;
];
tmp=Cases[data,a:F[{1,___},___]:>List@@a[[1;;,-1]],{0,∞}]//Flatten//Union;
If[tmp=!={},Throw[False]];
tmp=data gammaPat/.F[a___]:>Times@@MapIndexed[F[Sequence@@#1,#2[[1]]]&,{a}];
tmp=F[1]tmp/.F[0,a_,b_,c_,d_]:>a!F[((Gamma[vars[[c]]+δ[d]])^b),{δ[d],0,a}]((Gamma[vars[[c]]])^(-b))//Expand;
tmp=tmp//.{F[a_,b___]F[c_,d___]:>F[a c,b,d],F[a_,b___](c:((_Gamma)^(_.))):>F[a c ,b]};
res=If[("output"/.{method})===SeriesCoefficient,
tmp/.F[a_]:>a/.{δ->greδ,F->greSC},
tmp/.F[a_,b___]:>F[greSumToL[a],b]/.F[a_]:>a/.F->SeriesCoefficient
]/.greL[{{},{},{}}]->1;
True
]
];

greToD0[Derivative[a___][h_][b___],r_]:=
Module[{tmp,pos,i=0,j},
pos=Position[{a},Except[0],{-1},Heads->False];
tmp={b};
Do[
tmp[[Sequence@@pos[[j]]]]+=r[++i],
{j,Length[pos]}
];
greD[h@@tmp,Sequence@@Transpose[{Array[r,i],Extract[{a},pos]}]]
];

greToD0[a_,_]:=greD[a];

greToD[a_,r_]:=a/.v:Derivative[__][_][__]:>greToD0[v,r];

greToSC[a_,r_]:=
a/.v:Derivative[__][_][__]:>
Module[{tmp},
tmp=greToD0[v,r];
Times@@(#!&/@tmp[[2;;,2]])greSC[tmp[[1]],Sequence@@({#[[1]],0,#[[2]]}&/@tmp[[2;;]])]
]

greSplitL[in_greL]:=
Module[{tmp},
tmp=greLToSum[in];
tmp=If[tmp[[0]]===Times,List@@tmp,{tmp}];
tmp=tmp//.{a___,b_,c___,d_,e___}/;
Intersection[
Cases[b,greX[x_]:>x,{0,∞}],
Cases[d,greX[x_]:>x,{0,∞}]
]=!={}:>
{a,b d,c,e};
Times@@(greSumToL/@tmp)
];

greSplitL[in:Derivative[__][greL][__]]:=
Module[{tmp,r},
tmp=greToSC[in,r];
tmp=tmp/.a_greL:>greSplitL[a];
SeriesCoefficient@@tmp
];

If[!ValueQ[replTable],replTable=<<"gammaReduce_replTable.m"];

greLToSCInt[in:(_greL|Derivative[__][greL][__]),g_]:=
Module[{tmp,tmp2,i},
tmp=applyPermId2[in,g];
Catch[
Do[
tmp2=tmp[[i]]/.replTable;
If[tmp2=!=tmp[[i]],Throw[tmp2]],
{i,Length[tmp]}
];
in
]
];

mLimit[in_,rr_]:=
Module[{tmp},
If[FreeQ[in,rr[[1]]],Return[in]];
tmp=in/.rr//Quiet;
If[tmp=!=Indeterminate,Return[in]];
tmp=Limit[in,rr];
If[tmp[[0]]===Limit,Return[Indeterminate]];
tmp
];

fck[in_,r_,var_]:=
Module[{data,tmp,rep},
data=If[in[[0]]===Plus,List@@in,{in}];
tmp=data/.r//Quiet;
If[!MemberQ[tmp,0|_DirectedInfinity],Return[True]];
rep=Thread[var->#]&/@Tuples[{1,0},Length[var]];
Catch[
Scan[
Function[v,
tmp=data/.v/.{hpl2F1[_,_,_,0]->1,hpl2F1[a_,b_,c_,1]:>((Gamma[c-a-b]Gamma[c])/(Gamma[c-a]Gamma[c-b]))};
tmp=Fold[mLimit,#,r]&/@tmp;
If[MemberQ[tmp,_DirectedInfinity|_Indeterminate],Throw[False]];
],
rep
]//Quiet;
True
]
];

[収束性、零判定]

$gp[n_]:=$gp[n]=Function[f,Select[Range[0,(2^n)-1],Reverse[IntegerDigits[#,2,n]][[f]]==1&]]/@Range[n]+1;

atInf[in_,z_]:=in/.{
Gamma[a_.+b_.z]/;FreeQ[{a,b},z]:>(Sqrt[2π])(E^(-b z))(z^(a+b z-1/2))(b^(a+b z-1/2)),
PolyGamma[0,a_.+b_.z]/;FreeQ[{a,b},z]:>Re[Log[b z]],
PolyGamma[n_Integer?Positive,a_.+b_.z]/;FreeQ[{a,b},z]:>(n-1)!(z^(-n))(b^(-n))
};

greConvergenceQ[in0:_greL|Derivative[__][greL][__],
assump:((Assumptions->_)|(Assumptions:>_)):(Assumptions:>$Assumptions)]:=
Module[{tmp,ret,var,chk},
ret/;Catch[
tmp=in0//greLToSum;
var=Cases[tmp,_greX,{0,∞}]//Union;
Do[
chk=atInf[tmp,greX[i]];
chk=Limit[chk greX[i],greX[i]->∞,assump];
If[chk[[0]]===Limit||chk===Indeterminate,Throw[False]];
If[chk=!=0,ret=False;Throw[True]],
{i,Length[var]}
];
ret=True;
True
]
];

greZeroQ[in:(greL|Derivative[__][greL])[{_List,_List,{0...}}...]]:=MemberQ[(in/.greϵ->0)[[1,2]],0];

[2重和の評価]

hplIExpand[in_]:=nestedDistribute[in,hplI,Alternatives@@#[[2;;,1]]&,fx]

intExpand[in_]:=nestedDistribute[in,int,Alternatives@@Join[List@@#[[2;;-2,1]],#[[-1]]]&,fx]

fx[in_]:=MatchQ[in,(1+_)^_?(!NumericQ[#]&)]

nestedDistribute[in_,head_,fun_,fun2_]:=
Module[{F},
FixedPoint[
#/.orig:head[a_,b___]:>
With[{ret=Distribute[head[Expand[a/.t_?fun2:>F[t]]]/.F->Identity,Plus,head]/.head[x_]:>head[x,b]/.
orig2:head[a2_Times,b2___]:>Select[a2,FreeQ[#,fun[orig2]]&]*
head[Select[a2,!FreeQ[#,fun[orig2]]&],b2]},
ret/;ret=!=orig
]&,
in
]
];

hplIApplyKummerAll//Clear;
hplIApplyKummerAll[in_hplI,δ_,Δ_]:=
Module[{ts,intVar,intVarTo1,allZsumVar,firstZsumVar,gammaVar,zsumVar,denVar,
gammaVar1,gammaVar2,gammaVar3,zsumVar1,zsumVar2,zsumVar3,sVar,tVar,fff,
zCoeff,allZc,gc1,gc2,gc3,F,vc,vcp,ts2,ts3,pr,qr,dum,i,dRep},
ts=in/.(((-1+a_))^(b_.)):>(((1-a))^b)(((-1))^b);
intVar=(List@@ts)[[2;;,1]];
intVarTo1=Cases[ts,{a_,0,1}:>a];
allZsumVar=Cases[ts,{a_,0,b:Except[1]}:>{a,b}]//Flatten//Union;
firstZsumVar=Intersection[intVarTo1,allZsumVar];
gammaVar=Complement[intVarTo1,firstZsumVar];
zsumVar=Flatten[NestWhileList[Cases[ts,{a_,0,#[[1]]}:>a]&,{#},#=!={}&]]&/@firstZsumVar;
denVar=Cases[ts[[1]],(1/((1-a_Times))):>List@@a,{0,∞}];
denVar=PadRight[denVar,2,{{}}];
gammaVar1=SortBy[Intersection[Complement@@denVar,gammaVar],
Boole[!FreeQ[ts[[1]],((((1-#))^(a_.))/;!IntegerQ[a/.Thread[Variables[a]->0]]||!FreeQ[a,δ])|(((#)^(a_.))/;!FreeQ[a,δ])]]&];
gammaVar2=SortBy[Intersection[Complement@@Reverse[denVar],gammaVar],
Boole[!FreeQ[ts[[1]],((((1-#))^(a_.))/;!IntegerQ[a/.Thread[Variables[a]->0]]||!FreeQ[a,δ])|(((#)^(a_.))/;!FreeQ[a,δ])]]&];
gammaVar3=SortBy[Intersection[Intersection@@denVar,gammaVar],
Boole[!FreeQ[ts[[1]],((((1-#))^(a_.))/;!IntegerQ[a/.Thread[Variables[a]->0]]||!FreeQ[a,δ])|(((#)^(a_.))/;!FreeQ[a,δ])]]&];
zsumVar1=Select[zsumVar,Intersection[#,denVar[[1]]]=!={}&&Intersection[#,denVar[[2]]]==={}&];
zsumVar2=Select[zsumVar,Intersection[#,denVar[[1]]]==={}&&Intersection[#,denVar[[2]]]=!={}&];
zsumVar3=Select[zsumVar,Intersection[#,denVar[[1]]]=!={}&&Intersection[#,denVar[[2]]]=!={}&];
zCoeff=Cases[ts[[1]],(a:Alternatives@@allZsumVar)-p_:>{a,p},{0,∞}];
allZc=If[#[[0]]===Times,1/List@@#,1/#]&/@zCoeff[[1;;,2]]//Flatten//Union;
gc1=Times@@Complement[Complement@@denVar,intVar,allZc];
gc2=Times@@Complement[Complement@@Reverse[denVar],intVar,allZc];
gc3=Times@@Complement[Intersection@@denVar,intVar,allZc];
ts2=Fold[
hplI[(#1[[1]]/.(1/((#2[[1]]-#2[[2]]))):>((-1)/(#2[[2]]))(1/(1-dum))/.#2[[1]]->Times@@#2/.dum->#2[[1]])#2[[2]],
Sequence@@(#1[[2;;]]/.{#2[[1]],0,p_}:>{#2[[1]],0,p/#2[[2]]})]&,
ts,
zCoeff
]//PowerExpand;
ts2[[1]]=ts2[[1]]Times@@(((#)^(Δ[#]))(((1-#))^(Δ[1-#]))&/@intVar);
sVar=Join[gammaVar1,gammaVar2,zsumVar1,zsumVar2,gammaVar3,zsumVar3]//Flatten;
tVar={};
If[#=!={},
Attributes[Δ]={HoldAll};
ts2=Developer`ReplaceAllUnheld[ts2,hplIApplyKummer[#[[-1]]]]//hplIExpand;
Attributes[Δ]={};
ts2=ts2/.Δ[#[[-1]]]->Δ[#[[-1]]]-Δ[1-#[[-1]]];
fff=DeleteCases[Select[denVar,Function[p,!FreeQ[p,#[[-1]]]]]//Flatten,#[[-1]]];
Function[p,ts2=ts2/.Δ[p]->Δ[p]+Δ[#[[-1]]]]/@fff;
If[Length[fff]==1,ts2=ts2/.Δ[1-fff[[1]]]->Δ[1-fff[[1]]]-Δ[1-#[[-1]]]];
#=Most[#];
]&/@{Unevaluated[gammaVar2],Unevaluated[gammaVar1]};
vc=List[
If[gammaVar3=!={},Thread[F[Most[gammaVar3], Rest[gammaVar3]]]//Reverse,{}],
If[zsumVar3=!={}&&gammaVar3=!={},Thread[F[First[zsumVar3],First[gammaVar3]]],{}],
If[zsumVar3=!={},MapThread[Thread[F[#1,#2]]&,{Rest[zsumVar3],Last/@Most[zsumVar3]}],{}],
If[gammaVar2=!={}&&zsumVar3=!={},F[Last[gammaVar2],Last@Last@zsumVar3],{}],
If[gammaVar2=!={}&&zsumVar3==={}&&gammaVar3=!={},F[Last[gammaVar2],First[gammaVar3]],{}],
If[gammaVar2=!={},Thread[F[Most[gammaVar2], Rest[gammaVar2]]]//Reverse,{}],
If[zsumVar2=!={}&&gammaVar2=!={},Thread[F[First[zsumVar2],First[gammaVar2]]],{}],
If[zsumVar2=!={}&&gammaVar2==={}&&zsumVar3=!={},Thread[F[First[zsumVar2],Last@Last@zsumVar3]],{}],
If[zsumVar2=!={}&&gammaVar2==={}&&zsumVar3==={}&&gammaVar3=!={},Thread[F[First[zsumVar2],First[gammaVar3]]],{}],
If[zsumVar2=!={},MapThread[Thread[F[#1,#2]]&,{Rest[zsumVar2],Last/@Most[zsumVar2]}],{}],
If[gammaVar1=!={}&&zsumVar3=!={},F[Last[gammaVar1],Last@Last@zsumVar3],{}],
If[gammaVar1=!={}&&zsumVar3==={}&&gammaVar3=!={},F[Last[gammaVar1],First[gammaVar3]],{}],
If[gammaVar1=!={},Thread[F[Most[gammaVar1], Rest[gammaVar1]]]//Reverse,{}],
If[zsumVar1=!={}&&gammaVar1=!={},Thread[F[First[zsumVar1],First[gammaVar1]]],{}],
If[zsumVar1=!={}&&gammaVar1==={}&&zsumVar3=!={},Thread[F[First[zsumVar1],Last@Last@zsumVar3]],{}],
If[zsumVar1=!={}&&gammaVar1==={}&&zsumVar3==={}&&gammaVar3=!={},Thread[F[First[zsumVar1],First[gammaVar3]]],{}],
If[zsumVar1=!={},MapThread[Thread[F[#1,#2]]&,{Rest[zsumVar1],Last/@Most[zsumVar1]}],{}]
]//Flatten;
vc=vc/.F->List;
Attributes[Δ]={HoldAll};
ts3=ts2/.a_hplI:>Fold[hplI[(Developer`ReplaceAllUnheld[#1[[1]],#2[[1]]->Divide@@#2])/#2[[2]],Sequence@@(#1[[2;;]]/.#2[[1]]->Divide@@#2)]&,a,vc];
Attributes[Δ]={};
ts3=ts3//.{{((a_)/(b_)),0,1}:>{a,0,b},{((a_)/(b_)),0,((c_)/(b_))}:>{a,0,c},{((a_)/(p:gc1|gc2|gc3)),0,b_}:>{a,0,p b}}//PowerExpand;
ts3=ts3/.a_hplI:>hplISortIntRegion[a]/.(((1-((a_)/(b_))))^(c_.)):>(((b-a))^c)(b^(-c));
dRep=Δ[#[[2]]]->Δ[#[[2]]]+Δ[#[[1]]]+Δ[1-#[[1]]]&/@vc;
ts3=Fold[ReplaceAll,ts3,dRep];
fff=Cases[ts3,_Δ,{0,∞}]//Union;
ts3=ts3/.Thread[fff->Array[Δ,Length[fff]]];
{ts3,sVar,tVar}
];

hplISortIntRegion[hplI[a_,b___]]:=
Module[{tmp},
tmp={b};
tmp=tmp//.{x___,{s_,t__},y___,{u_,v__},z___}/;!FreeQ[{t},u]:>{x,{u,v},{s,t},y,z};
hplI[a,Sequence@@tmp]
];

hplIApplyKummer[x_]:={
hplI[(((1-x))^(b_.)) (x^(a_.))(((1-x y_))^(c_.))j_.,t___,{x,0,1},s___]/;FreeQ[{t,s,j},x]&&(!IntegerQ[1+a+b]||1+a+b>0)&&
(!IntegerQ[1+a]||1+a>0)&&(!IntegerQ[1+b]||1+b>0):> ((Gamma[1+a] Gamma[1+b]hplI[j(((1-y))^(1+b+c)) (y^(-1-a-b))(((1-x))^(-1-b)) (x^(1+a+b+c)) (((-x+y))^(-1-c)),t,s,{x,0,y}])/(Gamma[-c] Gamma[2+a+b+c])),
hplI[(((1-x))^(b_.)) (((1-x y_))^(c_.))j_.,t___,{x,0,1},s___]/;FreeQ[{t,s,j},x]&&(!IntegerQ[1+b]||1+b>0):>(( Gamma[1+b]hplI[j(((1-y))^(1+b+c)) (y^(-1-b))(((1-x))^(-1-b)) (x^(1+b+c)) (((-x+y))^(-1-c)),t,s,{x,0,y}])/(Gamma[-c] Gamma[2+b+c])),
hplI[ (x^(a_.))(((1-x y_))^(c_.))j_.,t___,{x,0,1},s___]/;FreeQ[{t,s,j},x]&&(!IntegerQ[1+a]||1+a>0):>((Gamma[1+a] hplI[j(((1-y))^(1+c)) (y^(-1-a)) (((1-x))^(-1)) (x^(1+a+c)) (((-x+y))^(-1-c)),t,s,{x,0,y}])/(Gamma[-c] Gamma[2+a+c]))
};

hplIApplySTransform[in:{_,_List,{}}]:=
Module[{tmp,dummy},
tmp=in[[1]]/.a_hplI:>Fold[hplSTrf[#1,dummy,dummy,#2,True]&,a,Intersection[in[[2]],List@@a[[2;;,1]]]];
tmp=tmp/.a_ hplI[b_,c___]:>hplI[a b,c];
tmp=tmp/.{(((4a_)/(((1+a_))^2))-((4b_)/(((1+b_))^2))):>((4(a-b)(1-a b))/((((1+a))^2)(((1+b))^2))),{((4a_)/(((1+a_))^2)),0,((4b_)/(((1+b_))^2))}:>{a,0,b},{((4a_)/(((1+a_))^2)),((4b_)/(((1+b_))^2)),1}:>{a,b,1},{((4a_)/(((1+a_))^2)),0,1}:>{a,0,1}};
tmp=tmp//PowerExpand;
tmp/.((a_)^(b_)):>(a^(Expand[b]))
];

hplIApplyTTransform[in:{_,_List,{}}]:=
Module[{tmp,dummy},
tmp=in[[1]]/.a_hplI:>Fold[hplTTrf[#1,dummy,dummy,#2,True]&,a,Intersection[in[[2]],List@@a[[2;;,1]]]];
tmp=tmp/.a_ hplI[b_,c___]:>hplI[a b,c];
tmp=tmp/.{(((((1+a_))^2)/(4a_))-((((1+b_))^2)/(4b_))):>-(((a-b)(1-a b))/(4a b)),{((((1+a_))^2)/(4a_)),0,((((1+b_))^2)/(4b_))}:>{a,-1,b},{((((1+a_))^2)/(4a_)),((((1+b_))^2)/(4b_)),1}:>{a,b,1},{((((1+a_))^2)/(4a_)),0,1}:>{a,-1,1}};
tmp=tmp//PowerExpand;
tmp/.((a_)^(b_)):>(a^(Expand[b]))
];

hplISplitIntRegion[a_hplI]:=
Module[{reg,start,p3,tmp,tmp2,i,p12},
reg=(List@@a)[[2;;]];
start={_,0,1};
p3={};
p3=
Reap[
While[True,
tmp=Cases[reg,start];
If[Length[tmp]≠1,Break[]];
Sow[tmp[[1,1]]];
reg=DeleteCases[reg,tmp[[1]]];
start={_,0,tmp[[1,1]]};
];
][[2]]//Flatten;
p12={{},{}};
If[Length[tmp]==2,
p12=Table[
Reap[
Sow[tmp[[i,1]]];
reg=DeleteCases[reg,tmp[[i]]];
start={_,0,tmp[[i,1]]};
While[True,
tmp2=Cases[reg,start];
If[Length[tmp2]≠1,Break[]];
Sow[tmp2[[1,1]]];
reg=DeleteCases[reg,Alternatives@@tmp2];
start={_,0,tmp2[[1,1]]};
];
][[2]]//Flatten,
{i,2}
]
];
If[reg=!={},
Print[a];
Print["hplISplitIntRegionError"];
Abort[]
];
{p3,p12[[1]],p12[[2]]}
];

hplIFindSingularVar[in_hplI,δ_]:=
Module[{facl,ivarQ,tmp},
facl=If[in[[1,0]]===Times,List@@in[[1]],{in[[1]]}];
ivarQ=MatchQ[#,Alternatives@@in[[2;;,1]]]&;
tmp=Cases[in,((((a_?ivarQ)|(1-a_?ivarQ)|(b_?ivarQ-a_?ivarQ)))^(b_))/;!FreeQ[b,δ]:>a,{0,∞}]
];

$intRegularize//Clear;
$intRegularize[int_,i_,k_]:=$intRegularize[int,i,k]=
Module[{s,t,integrand,d,dummy,res,bef,aft},
s=int[[i+1,1]];
t=If[int[[i+1,2]]=!=0,int[[i+1,2]],If[int[[i+1,3]]=!=1,int[[i+1,3]],dummy]];
integrand=int[[1]](s^(d[1]))(((1-s))^(d[2]))(((1+s))^(d[3]))(((t-s))^(d[4]))(((1-t s))^(d[5]))//.x_.(((a_-1))^(p_.)):>(x(((-1))^p))(((1-a))^p)/.
a_.((s_)^(x_.))(((1-s_))^(y_.))(((1+s_))^(v_.))(((t_-s_))^(z_.))(((1-t_ s_))^(w_.))/;FreeQ[a/._Log|_hplHC->1,s]:>{a,{s,t,x,y,v,z,w}};
If[integrand[[0]]===List,
Return[$intRegularize2[Sequence@@(integrand/._d->0),Sequence@@int[[2;;]],k]]
];
Abort[];
];

$intRegularize2//Clear;
$intRegularize2[a_,{s_,t_,x_,y_,v_,z_,w_},region___,k_]:=$intRegularize2[a,{s,t,x,y,v,z,w},region,k]=
Module[{xi,yi,vi,zi,wi,integ,tmp,xp,yp,vp,zp,wp,ord,ord2,dir,
reg1x,reg1y,reg1v,reg1z,reg1w,reg2x,reg2y,reg2v,reg2w,reg2z,d},
{xi,yi,vi,zi,wi}={x,y,v,z,w}/.Alternatives@@k->0;
{xp,yp,vp,zp,wp}={x,y,v,z,w};
integ=a (s^xp)(((1-s))^yp)(((1+s))^vp)(((t-s))^zp)(((1-t s))^wp);
If[!NumericQ[#],Print["intRegularize2Error"];Abort[]]&/@{xi,yi,vi,zi,wi};
dir=Cases[{region},{s,l_,u_}/;l==0||l==-1||u==1:>Union[Cases[{l,u},_Integer]]];
If[Length[dir]≠1,Print["intRegularize2Error"];Abort[]];
dir=dir[[1]];
ord=Exponent[a/.Log[1+_.s]:>Log[1+s]/.hplHC[j_,(s^(p_.)),_]:>Log[1+s]^(p Count[j,Except[0]]),Log[1+s]];
ord2=Exponent[a/.f:hplHC[{0..},_,_]:>hplHPLToLog[a],Log[s]];
Switch[dir,
{0},
{reg1x=xi>-ord-1,reg1y=yi≥-1,reg1v=vi≥-1,reg1z=zi≥0,reg1w=wi≥0};(*積分しても問題にならないか*)
{reg2x=xi>-ord||x===0,reg2y=yi≥0||y===0,reg2v=vi≥0||v===0,reg2z=zi≥1||z===0,reg2w=wi≥1||w===0};(*部分積分しても問題にならないか*),
(*else*)
{1},
{reg1x=True,reg1y=yi≥If[ord2≥1,-1,0],reg1v=vi≥-1,reg1z=zi≥0,reg1w=wi≥0};(*積分しても問題にならないか*)
{reg2x=True,reg2y=yi≥If[ord2≥1,0,1]||y===0,reg2v=vi≥0||v===0,reg2z=zi≥1||z===0,reg2w=wi≥1||w===0};(*部分積分しても問題にならないか*),
(*else*)
{0,1},
If[!((FreeQ[a,Log[1+j_.s]/;Variables[j]=!={}])&&z===w===0),Print["intRegularize2Error"];Abort[]];
{reg1x=xi>-ord-1,reg1y=yi>-ord2-1,reg1v=True,reg1z=True,reg1w=True};(*積分しても問題にならないか*)
{reg2x=xi>-ord||x===0,reg2y=yi>-ord2||y===0,reg2v=True,reg2z=True,reg2w=True};(*部分積分しても問題にならないか*),
(*else*)
True,Print["intRegularize2Error"];Abort[]
];
tmp=Which[
reg1x&&reg1y&&reg1v&&reg1z&&reg1w,hplI[integ,region],
reg2x&&reg2y&&reg2v&&reg2z&&!reg1w,intByParts[hplI[integ ,region],s,1-t s],
reg2x&&reg2y&&reg2v&&!reg1z&&reg2w,intByParts[hplI[integ ,region],s,t-s],
reg2x&&reg2y&&!reg1v&&reg2z&&reg2w,intByParts[hplI[integ ,region],s,1+s],
reg2x&&!reg1y&&reg2v&&reg2z&&reg2w,intByParts[hplI[integ ,region],s,1-s],
!reg1x&&reg2y&&reg2v&&reg2z&&reg2w,intByParts[hplI[integ ,region],s,s],
reg2x&&reg2y&&reg2v&&(!reg1z||!reg1w),-hplI[((integ t(t-s))/((1-t)(1+t))),region]+hplI[((integ(1-t s))/((1-t)(1+t))),region],
reg2x&&reg2y&&reg2z&&(!reg1v||!reg1w),hplI[((integ t(1+s))/(1+t)),region]+hplI[((integ(1-t s))/(1+t)),region],
reg2x&&reg2v&&reg2z&&(!reg1y||!reg1w),-hplI[((integ t(1-s))/(1-t)),region]+hplI[((integ(1-t s))/(1-t)),region],
reg2y&&reg2v&&reg2z&&(!reg1x||!reg1w),hplI[integ  t s,region]+hplI[integ(1-t s),region],
reg2x&&reg2y&&reg2w&&(!reg1v||!reg1z),hplI[((integ (1+s))/(1+t)),region]+hplI[((integ(t-s))/(1+t)),region],
reg2x&&reg2v&&reg2w&&(!reg1y||!reg1z),hplI[((integ (1-s))/(1-t)),region]-hplI[((integ(t-s))/(1-t)),region],
reg2y&&reg2v&&reg2w&&(!reg1x||!reg1z),hplI[((integ  s)/t),region]+hplI[((integ(t-s))/t),region],
reg2x&&reg2z&&reg2w&&(!reg1y||!reg1v),(1/2)hplI[integ (1-s),region]+(1/2)hplI[integ(1+s),region],
reg2y&&reg2z&&reg2w&&(!reg1x||!reg1v),-hplI[integ  s,region]+hplI[integ(1+s),region],
reg2v&&reg2z&&reg2w&&(!reg1x||!reg1y),hplI[integ  s,region]+hplI[integ(1-s),region],
reg2x&&reg2y&&(!reg1v||!reg1z||!reg1w),hplI[((integ (t-s) (1-s t))/(((1+t))^2)),region]+hplI[((integ (1+s) (1-s t))/((1-t) (((1+t))^2))),region]-hplI[((integ (t^2)(1+s)  (t-s))/((1-t) (((1+t))^2))),region],
reg2x&&reg2v&&(!reg1y||!reg1z||!reg1w),-hplI[((integ (t-s) (1-s t))/(((1-t))^2)),region]+hplI[((integ (1-s) (1-s t))/((((1-t))^2) (1+t))),region]+hplI[((integ (t^2) (1-s) (t-s))/((((1-t))^2) (1+t))),region],
reg2y&&reg2v&&(!reg1x||!reg1z||!reg1w),hplI[((integ (t-s) (1-s t))/t),region]+hplI[((integ s (1-s t))/((1-t) t (1+t))),region]-hplI[((integ (t^2)s(t-s))/((1-t) (1+t))),region],
reg2x&&reg2z&&(!reg1y||!reg1v||!reg1w),(1/2)hplI[((integ (1+s) (1-s t))/(1-t)),region]+(1/2)hplI[((integ (1-s) (1-s t))/( 1+t)),region]-hplI[((integ (t^2)(1-s) (1+s))/((1-t) (1+t))),region],
reg2y&&reg2z&&(!reg1x||!reg1v||!reg1w),hplI[integ (1+s) (1-s t),region]-hplI[((integ s (1-s t))/(1+t)),region]+hplI[((integ (t^2)s (1+s))/(1+t)),region],
reg2v&&reg2z&&(!reg1x||!reg1y||!reg1w),hplI[integ (1-s) (1-s t),region]+hplI[((integ s (1-s t))/(1-t)),region]-hplI[((integ (t^2)s(1-s) )/(1-t)),region],
reg2x&&reg2w&&(!reg1y||!reg1v||!reg1z),-(1/2)hplI[((integ (1+s) (t-s))/(1-t)),region]+(1/2)hplI[((integ (1-s) (t-s))/(1+t)),region]+hplI[((integ (1-s) (1+s))/((1-t) (1+t))),region],
reg2y&&reg2w&&(!reg1x||!reg1v||!reg1z),hplI[((integ (1+s) (t-s))/t),region]-hplI[((integ s (t-s))/(1+t)),region]+hplI[((integ s (1+s))/(t (1+t))),region],
reg2v&&reg2w&&(!reg1x||!reg1y||!reg1z),hplI[((integ (1-s) (t-s))/t),region]-hplI[((integ s (t-s))/(1-t)),region]+hplI[((integ s(1-s))/(t(1-t))),region],
reg2z&&reg2w&&(!reg1x||!reg1y||!reg1v),hplI[integ (1-s) (1+s),region]+(1/2)hplI[integ s (1+s),region]-(1/2)hplI[integ s(1-s),region],
reg2x&&(!reg1y||!reg1v||!reg1z||!reg1w),-(1/2)hplI[((integ (1+s) (t-s) (1-s t))/(((1-t))^2)),region]+(1/2)hplI[((integ (1-s) (t-s) (1-s t))/(((1+t))^2)),region]+hplI[((integ (1-s) (1+s) (1-s t))/((((1-t))^2) (((1+t))^2))),region]+hplI[((integ (t^3)(1-s) (1+s)  (t-s))/((((1-t))^2) (((1+t))^2))),region],
reg2y&&(!reg1x||!reg1v||!reg1z||!reg1w),hplI[((integ (1+s) (t-s) (1-s t))/t),region]-hplI[((integ s (t-s) (1-s t))/(((1+t))^2)),region]+hplI[((integ s (1+s) (1-s t))/(t(1-t)(((1+t))^2))),region]-hplI[((integ (t^3)s (1+s)(t-s))/((1-t) (((1+t))^2))),region],
reg2v&&(!reg1x||!reg1y||!reg1z||!reg1w),hplI[((integ (1-s) (t-s) (1-s t))/t),region]-hplI[((integ s (t-s) (1-s t))/(((1-t))^2)),region]+hplI[((integ s(1-s)(1-s t))/(t(((1-t))^2) (1+t))),region]+hplI[((integ (t^3) s(1-s)(t-s))/((((1-t))^2) (1+t))),region],
reg2z&&(!reg1x||!reg1y||!reg1v||!reg1w),hplI[integ (1-s) (1+s) (1-s t),region]+(1/2)hplI[((integ s (1+s) (1-s t))/(1-t)),region]-(1/2)hplI[((integ s(1-s)(1-s t))/(1+t)),region]-hplI[((integ (t^3)s(1-s)(1+s))/((1-t) (1+t))),region],
reg2w&&(!reg1x||!reg1y||!reg1v||!reg1z),hplI[((integ (1-s) (1+s) (t-s))/t),region]-(1/2)hplI[((integ s (1+s) (t-s))/(1-t)),region]-(1/2)hplI[((integ s(1-s)(t-s))/(1+t)),region]+hplI[((integ s(1-s)(1+s))/(t(1-t)(1+t))),region],
True,hplI[((integ (1-s) (1+s) (t-s) (1-s t))/t) ,region]-(1/2)hplI[((integ s (1+s) (t-s) (1-s t))/(((1-t))^2)),region]-(1/2)hplI[((integ s(1-s)(t-s) (1-s t))/(((1+t))^2)),region]+hplI[((integ s(1-s)(1+s) (1-s t))/(t(((1-t))^2)(((1+t))^2))),region]+hplI[((integ (t^4) s(1-s)(1+s)(t-s))/((((1-t))^2) (((1+t))^2))),region]
];
tmp=Distribute@hplI[Distribute[tmp (t^(d[1]))(1-t)^(d[2])(1+t)^(d[3])]//.c_.(d_-1)^(p_.):>(c(-1)^p)(1-d)^p];
tmp=tmp//.f_ hplI[g_,h___]:>hplI[f g,h]/.hplI[hplI[f___]]:>hplI[f];
tmp=tmp/.hplI[n_,m___]:>hplI[n//rationalize,m];
hplIExpand[tmp]/._d->0
]

series2//Clear;
series2[in_,x_]:=
Module[{tmp,tmp2,coeff},
Fold[
Function[{u,v},
tmp=u;
tmp2=tmp/.a_hplZC:>(a/.∞->hplN)/.v->0//Quiet;
If[!FreeQ[tmp2,DirectedInfinity|Indeterminate]||tmp2===0//Quiet,
If[tmp[[0]]===Times,
coeff=Select[tmp,FreeQ[#,v]&];
tmp2=Select[tmp,!FreeQ[#,v]&],
(*else*)
coeff=1;
tmp2=tmp
];
tmp2=coeff SeriesCoefficient[tmp2,{v,0,0}]//Normal//Together,
(*else*)
tmp2=tmp/.v->0
];
tmp2
],
in//Together,Reverse[x]
]
];

rationalize[in_]:=
in//.r_.(((t_.((x_)^(p_))+s_))^(q_Integer?Negative))/;IntegerQ[Log2[p]]:>
With[{j=Expand[(((1-(t/s)(x^p)))^(-q))]},If[j[[0]]===Plus,r (s^q)(((1-((t^2)/(s^2))(x^(2p))))^q)#&/@j,r (s^q)(((1-((t^2)/(s^2))(x^(2p))))^q)j]];

intByParts[in_hplI,var_,fac_]:=
Module[{region,pos,tmp,f1,f2,i,δ},
region=List@@in[[2;;]];
pos=Position[region[[1;;,1]],var];
If[pos==={},Print[{in,var,fac}];Print["intByPartsError"];Abort[]];
pos=pos[[1,1]];
tmp=If[in[[1,0]]===Times,List@@in[[1]],{in[[1]]}];
f1=Select[tmp,MatchQ[#,(fac^(_.))]&][[1]];
f2=Times@@DeleteCases[tmp,f1];
tmp=(hplI[Integrate[f1,var]f2,Sequence@@Delete[region,pos]]/.var->region[[pos,3]]-δ)-
(hplI[Integrate[f1,var]f2,Sequence@@Delete[region,pos]]/.var->region[[pos,2]]+δ)-
(hplI[Integrate[f1,var]D[f2,var],Sequence@@region])-
Sum[
(hplI[Integrate[f1,var]f2 D[region[[i,3]],var],Sequence@@Delete[region,i]]/.region[[i,1]]->region[[i,3]]-δ)-
(hplI[Integrate[f1,var]f2 D[region[[i,2]],var],Sequence@@Delete[region,i]]/.region[[i,1]]->region[[i,2]]+δ),
{i,pos+1,Length[region]}
]//hplIExpand;
tmp=tmp/.((((_.δ))^(_.))((Log[_.δ])^(_.))->0)/.δ->0/.hplI[a_]:>a;
If[!FreeQ[tmp,Indeterminate],Print[{in,var,fac}];Print["intByPartsError"];Abort[]];
tmp/.((0^(a_))->0)/.hplI[0,___]->0/.hplI[___,{_,0,0},___]->0
];

logCoeff[f_,0]=0;
logCoeff[f_,1]:=-(1/f);
logCoeff[f_,n_Integer?NonNegative]:=(((((-1))^n)n!)/f)hplZC[ConstantArray[1,n-1],f-1,ConstantArray[1,n-1],ConstantArray[0,n-1]]

sumSimp[in_]:=
FixedPoint[
#/.a_ sum[b_,c___]:>sum[a b,c]/.
sum[a_,b___]:>Distribute@sum[Expand[a],b]/.
sum[sum[a_,b___],c___]:>sum[a,c,b]&,
in
]

pnlim[in_,x_]:=
Module[{tmp,nn},
tmp=PowerExpand[in]/.((a_)^(b_)):>(a^(Expand[b]));
tmp=Refine[tmp,hplN∈Integers];
tmp=tmp/.(x^(_.+s_.hplN))/;s>0->0/.((p_)^(q_))/;Abs[p]==1&&!FreeQ[q,hplN]:>(p^(q/.hplN->nn))//Expand;
tmp=If[tmp[[0]]===Plus,List@@tmp,{tmp}];
tmp=If[
Limit[#/hplN/._hplZC->1,hplN->∞]==0&&!FreeQ[#,hplZC[_,_?(FreeQ[#,hplN]&),_,_?(!FreeQ[#,hplN]&)]],
0,#
]&/@tmp;
tmp=Normal[Series[#,{hplN,∞,0}]]&/@tmp//Total;
If[!FreeQ[tmp,nn],Print[{in,x}];Print["pnlimError"];Abort[]];
tmp
];

Clear[ts1,$ts1];
ts1[in_hplI]:=$ts1[in]=
Module[{var,tmp,c,d,bef,aft,cond,exp,ks,coeff,lcm,r1,r2,pat},
If[Length[in]>2,Print[in];Print["ts1Error"];Abort[]];
var=DeleteCases[Variables[List@@in[[2,2;;]]],in[[2,1]]];
If[!MatchQ[var,{_Symbol}],Print[in];Print["ts1Error"];Abort[]];
var=var[[1]];
tmp=Cases[in[[1]],((Log[1+x_])^(_.))/;!FreeQ[x,in[[2,1]]],{0,∞}]//Union;
c=d=1;
bef=tmp//.{
{a___,Log[p_],b___,((Log[p_])^(q_)),c___}:>{a,b,((Log[p])^q),c},
{a___,((Log[p_])^(q_)),b___,Log[p_],c___}:>{a,b,((Log[p])^q),c}
};
aft=bef/.((Log[1+p_])^(q_.)):>With[{cc=c++},logCoeff[j[cc],q]PowerExpand[(((-p))^(j[cc]))]];
tmp=tmp/.Thread[bef->aft];
tmp=in[[1]] Times@@aft/Times@@bef;
tmp=hplIntToSum[tmp,VerifyConvergence->False];
tmp=tmp/.(((1+s_))^(p_Integer?Negative)):>With[{cc=c++},(Pochhammer[j[cc],-1-p] (((-1))^(j[cc]-1)) (s^(j[cc]-1)))/(-1-p)!];tmp=tmp/.x_hplZC:>With[{cd=d++},Block[{Sum=sum,newvarList=(Array[k[cd],#2]&)},hplZSum@@x]];
tmp=tmp//.((a_sum)^(b_Integer?Positive)):>Product[With[{cd=d++},a/.k[_]->k[cd]],{b}];
tmp=Distribute@sum[tmp//Expand,Sequence@@Array[{j[#],1,∞}&,c-1]];
tmp=tmp//sumSimp//PowerExpand;
tmp=tmp/.sum[a_,b___]:>Distribute@sum[hplIntegrate[a,Sequence@@in[[2;;]]]//hplLogExtract//Expand,b];
tmp=tmp/.((0^(_))|hplHC[_,0,_]->0)/.sum[0,__]->0/.c:hplHC[{0..},__]:>hplHPLToLog[c];
tmp=hplIntToSum[tmp,VerifyConvergence->False];
tmp=tmp/.x_hplZC:>With[{cd=d++},Block[{Sum=sum,newvarList=(Array[k[cd],#2]&)},hplZSum@@x]];
tmp=tmp//.((a_sum)^(b_Integer?Positive)):>Product[With[{cd=d++},a/.k[_]->k[cd]],{b}]//sumSimp;
tmp=tmp/.a_sum:>
With[{vk=Union[Cases[a,j[_]|k[_][_],{0,∞}]]},
a/.Join[Thread[vk->Array[k,Length[vk]]]]];
tmp=tmp/.sum[a_,b__]:>sum[a,Sequence@@SortBy[{b},First]];
tmp=PowerExpand[tmp]/.((a_)^(b_)):>(a^(Expand[b]));
bef=Cases[tmp,_sum,{0,∞}]//Union;
aft=Table[
exp=Exponent[bef[[i,1]],var];
ks=Union[Cases[exp,k[_],{0,∞}]];
coeff=Coefficient[exp,ks];
If[Or@@Thread[coeff<0],Print[in];Print["ts1Error"];Abort[]];
lcm=LCM[Sequence@@Denominator[coeff],1];
bef[[i]]/.Thread[ks->ks /coeff/lcm],
{i,Length[bef]}
];
tmp=tmp/.Thread[bef->aft];
tmp=tmp//.x:sum[___,{p_k,1,_k n_-1}/;IntegerQ[1/n],___]:>(x/.p->p n);
tmp=tmp//.sum[a_,b___,{p_k n_,1,q_},c___]/;IntegerQ[1/n]:>
sum[a projω[1/n,0,1][p],b,{p,1,Expand[q/n]},c]/.sum[a_]:>a;
bef=Cases[tmp,_sum,{0,∞}]//Union;
aft=Table[
exp=Exponent[bef[[i,1]],var];
ks=Union[Cases[exp,k[_],{0,∞}]];
r1=Cases[bef[[i,2;;]],{j:Alternatives@@ks,_,_}];
pat=sum[bef[[i,1]],Sequence@@Join[r1,DeleteCases[List@@bef[[i,2;;]],Alternatives@@r1]]];
If[r1=={},
pat,
pat/.r1[[-1,1]]->2r1[[-1,1]]-Total[ks]
],
{i,Length[bef]}
];
tmp=tmp/.Thread[bef->aft]/.((a_)^(b_)):>(a^(Expand[b]));
tmp=tmp/.
sum[a_,b:PatternSequence[{_,1,∞}..],{s_Plus,1,∞},c___]:>
Module[{h,r,p},
h=Total[{b}[[1;;,1]]]+s;
r=Thread[{b}[[2;;,1]]->Differences[{b}[[1;;,1]]]];
p=Thread[{{b}[[1;;,1]],Prepend[1+{b}[[1;;-2,1]],1],h-1}];
sum[a/.r,{h,1,∞},Sequence@@p,c]
];
tmp=tmp/.sum[a_,b_,c__]/;!FreeQ[a,(var^(t_))/;!FreeQ[t,b[[1]]]]:>sum[sum2[a,c,b[[1]]],b];
tmp=tmp/.sum2[a_,b___]:>sum2[a,Sequence@@Join[#,DeleteCases[{b},Alternatives@@#]]&[Cases[{b},{_,_,∞}]]];
tmp=tmp/.sum[a__]:>sum[a,hplN]/.sum2->sum;
bef=Cases[tmp,__k,{0,∞}]//Union;
aft=newvarList[tmp,Length[bef]];
tmp=tmp/.Thread[bef->aft]//multipleZeta`private`myExpand;
tmp=tmp/.sum[a_,hplN]:>(a/.sum[b___]:>sum[b,hplN]);
ts2[tmp,var]
];

ts1P[in:hplI[a_,{b_,0,c_}]]:=
Module[{tmp},
tmp=$ts1[in/.{b->$x,c->$y}]/.{$x->b,$y->c};
If[tmp===Null||tmp[[0]]===$ts1,
tmp=ts1[in/.{b->$x,c->$y}]/.{$x->b,$y->c};
wri[{in,tmp}/.{b->$x,c->$y}];
];
tmp
];

Clear[ts2];
ts2[in_,v_]:=
Module[{bef,aft,tmp,F,pat,δ},
tmp=in;
bef=Cases[tmp,a_sum/;!FreeQ[a[[1]],sum],{0,∞}]//Union;
aft=bef/.sum[a_,b__]:>sum[a/.sum[x__]:>hplSum[x,Simplify->False,Expand->False],b];
aft=aft/.sum[a_,b__]:>sum[pnlim[a/.hplδ->δ,v]//multipleZeta`private`myExpand,b];
If[!FreeQ[aft[[1;;,1]],hplZC[_,_,_,_?(!FreeQ[#,hplN]&)]],Print[{in,v}];Print["ts2Error"];Abort[]];
aft=aft/.sum[a_,b_,c_]:>sum[a/.x_hplZC/;!FreeQ[x[[2]],hplN]:>hplMZVSimplify[x/.b[[1]]->0],b,c];
If[!FreeQ[aft,sum[_?(!FreeQ[#,hplN]&),___]],Print[{in,v}];Print["ts2Error"];Abort[]];
pat=MapThread[{#1/.((a_)^(b_.#2))/;Abs[(a^b)]>0->1,#2}&,{cdenom/@aft[[1;;,1]],aft[[1;;,2,1]]}];
pat=ArgMax[{#[[2]],If[#[[1,0]]===Times,Select[#[[1]],Function[p,!FreeQ[p,#[[2]]]]]/._greδ->0,#[[1]]]==0},#[[2]],Integers]&/@pat//Quiet;
If[!FreeQ[pat,ArgMax]||pat===Indeterminate,Print[{in,v}];Print["ts2Error"];Abort[]];
pat=If[IntegerQ[#]&&#>0,#,-∞]&/@pat;
aft=MapThread[
If[#2[[2,2]]≤#3,
ReplacePart[#2,{2,2}->#3+1]+Sum[#1[[1]],Evaluate[{#2[[2,1]],#2[[2,2]],#3}]],
#2
]&,
{bef,aft,pat}
]//Quiet;
aft=aft/.sum[a__,_Integer]:>sum[a,hplN]//multipleZeta`private`myExpand;
aft=aft/.sum[a_,b__]:>sum[pnlim[a,v],b];
tmp=tmp/.Thread[bef->aft];
bef=Cases[tmp,a_sum,{0,∞}]//Union;
aft=bef/.sum[a_,b_,hplN]:>sum[Refine[a,b[[1]]∈Integers]//multipleZeta`private`myExpand,b,hplN];
aft=Distribute/@aft;
aft=aft/.sum[x__]:>hplSum[x,Simplify->False,Expand->False]//multipleZeta`private`myExpand;
tmp=tmp/.Thread[bef->aft]/.x_hplZC/;!FreeQ[x[[2]],hplN]:>hplMZVSimplify[x];
tmp=pnlim[tmp,v];
Check[tmp=tmp/.a_hplZC:>(a/.hplN->∞),Print[{in,v}];Print["ts2Error"];Abort[]];
tmp=hplSumToInt[tmp,VerifyConvergence->False]//PowerExpand;
tmp=tmp/.hplHC[a__,1,b:{0..}]:>hplHC[a v^-Exponent[a[[-1]],v],v^-Exponent[a[[-1]],v],b];
tmp=pnlim[tmp,v];
Coefficient[Coefficient[tmp,δ,0],hplδ,0]//Expand
];

cdenom[in_]:=If[in[[0]]===Plus,Times@@Denominator[List@@in],Denominator[in]];

Clear[ts3,$ts3];
ts3[in:hplI[a_,{b_,0,1}]]:=$ts3[in]=
Module[{δ,tmp,cc,dd,bef,aft},
tmp=in/.hplI[x_,{y_,0,1}]/;!FreeQ[x,(y^(t_))/;!IntegerQ[t]&&IntegerQ[2t]]:>hplI[2y(x/.y->(y^2)),{y,0,1}];
tmp=tmp/.{(Sqrt[((x_)^2)]):>x,(1-((x_)^2)):>(1-x)(1+x),(1+((x_)^2)):>(1-I x)(1+I x),hplHC[x_,((y_)^2),z:{0..}]:>hplArgTransform[x,2,y]};
tmp=tmp//.Log[x_ y_]:>Log[x]+Log[y]//PowerExpand;
(* fix me *)
cc=If[(Denominator[tmp[[1]]]/b/.b->0)=!=0,0,δ]//Quiet;
dd=If[(Denominator[tmp[[1]]]/.b->1)=!=0,1,1-δ]//Quiet;
tmp=Check[hplIntegrate[tmp[[1]],{b,cc,dd}],hplIntegrate[$Failed]]//Expand;
If[FreeQ[tmp,hplIntegrate],
(* then *)
If[cc=!=0||dd=!=1,tmp=hplLimit[tmp,δ->0]],
(* else *)
(* fix me *)
tmp=hplLimit[ts1P[hplI[a,{b,0,1-δ}]],δ->0]
];
If[!FreeQ[tmp,δ],Print[in];Print["ts3Error"];Abort[]];
bef=Cases[tmp,_hplHC,{0,∞}]//Union;
aft=bef//hplLogExtract//hplIntToSum//hplMZVSimplify;
tmp/.Thread[bef->aft]//Expand
];

ts3P[in:hplI[a_,{b_,0,1}]]:=
Module[{tmp},
tmp=$ts3[in/.{b->$x}]/.{$x->b};
If[tmp===Null||tmp[[0]]===$ts3,
tmp=ts3[in/.{b->$x}]/.{$x->b};
wri[{in,tmp}/.{b->$x}];
];
tmp
];

ts4P[in:hplI[a_,b__]]:=
Module[{tmp,var,rule1,rule2},
var=Variables[{b}];
If[Length[var]>4,Print[in];Print["ts4PError"];Abort[]];
rule1=Thread[var->Take[{$x,$y,$z,$w},Length[var]]];
rule2=Thread[Take[{$x,$y,$z,$w},Length[var]]->var];
tmp=$ts4[in/.rule1]/.rule2;
If[tmp===Null||tmp[[0]]===$ts4,
tmp=ts4[in/.rule1]/.rule2;
wri[{in,tmp}/.rule1];
];
tmp
];

Clear[ts4,$ts4];
ts4[in:hplI[a_,b__,c_]]:=$ts4[in]=
Module[{tmp,coeff,bef,aft},
tmp=If[a[[0]]===Times,List@@a,{a}];
coeff=Times@@Select[tmp,FreeQ[#,c[[1]]]&];
tmp=Times@@Select[tmp,!FreeQ[#,c[[1]]]&];
tmp=ts1P[hplI[tmp,c]];
tmp=Expand[coeff tmp//hplShuffleAll];
tmp=Distribute@hplI[tmp,b];
bef=Cases[tmp,_hplI,{0,∞}]//Union;
aft=Switch[#,hplI[_,{_,0,1}],ts3P[#],hplI[_,_],ts1P[#],_,ts4P[#]]&/@bef;
tmp=tmp/.Thread[bef->aft];
Expand[tmp]
];

hf[hplI[a_,x___]]:=
Module[{var,fac,tmp,H},
var=Alternatives@@Variables[{x}];
fac=Cases[a,(((p_)^(q_))/;!FreeQ[p,var]&&!NumericQ[q])|(p:Log[q_/;!FreeQ[q,var]]|p_hplHC|p_G):>p,{0,∞}]//Union;
fac=Join[fac,Variables[{x}]];
tmp=a/.p:(((Alternatives@@fac))^(_.)):>H[p]//.H[p_]H[q_]:>H[p q]//Expand;
hplI[Collect[tmp,_H,F[Expand[PowerExpand[FunctionExpand[#]]]]&]//Expand,x]/.G|H->Identity//Distribute
];

greSCLToInt//Clear;
greSCLToInt[greSC[in:greL[_,{a_,b_,{},1},{c_,d_,{},1},{e_,f_,{},1}],p___],δ_,Δ_]/;
Length[a]==Length[b]≥1&&Length[c]==Length[d]≥1&&Length[e]==Length[f]≥1:=
Module[{integ,integ0,var1,var2,G,order,regparam},
integ=integ0=in/.
x_greL:>greLToInt[x,hplI->True]/.
x_hplI:>G@@hplIApplyKummerAll[x,δ,Δ]/.
x_ G[y_,z__]:>G[x y,z]/.
x_ hplI[y_,z___]:>hplI[x y,z];
order=hplISelectIntOrder[integ[[1]],δ];
integ=integ/.δ->Δ[0];
regparam=Table[
Cases[Cases[integ,((x_)^(y_))/;!FreeQ[x,Alternatives@@order[[i,1;;,1]]]&&!FreeQ[y,_Δ]:>y,{0,∞}],_Δ,{0,∞}]//Union,
{i,Length[order]}
];
If[MemberQ[Tally[Flatten[regparam]][[1;;,2]],_?(#>2&)],Print[{in,δ,Δ}];Print["greSCLToIntError"];Abort[]];
regparam=Table[Complement@@regparam[[i;;]],{i,Length[regparam]}];
If[!FreeQ[regparam[[2;;]],Δ[0]],Print[{in,δ,Δ}];Print["greSCLToIntError"];Abort[]];
integ=integ/.G[hplI[x_,y___],z__]:>hplIApplySTransform[{hplI[Fold[SeriesCoefficient,x,{p}],y],z}];
integ=integ//.{(((x_-y_))^(z_))(((y_-x_))^(w_.))/;IntegerQ[w]:>(((x-y))^(z+w))(((-1))^w),x_.(((y_-1))^(z_.)):>(x(((-1))^z))(((1-y))^z)}//hf;
integ=integ/.hplI[x_,y___]/;{y}===List@@integ0[[1,2;;]]:>hplI[x,order];
integ=integ/.hplI[x_,y:({_List..}..)]:>Fold[int[#1,Sequence@@#2[[1]],#2[[2]]]&,x,Reverse@Transpose[{y,regparam}]];
integ=integ//intExpand;
integ/.Thread[#->newvarList[integ,Length[#]]]&[Cases[integ,_greX,{0,∞}]//Union]
];

hplISelectIntOrder[in_hplI,δ_]:=
Module[{regseq,singular,i,pair,tmp1,tmp2,ret,F},
regseq=hplISplitIntRegion[in];
singular=hplIFindSingularVar[in,δ];
ret=Which[
(*if*)
regseq[[1]]==={}||(regseq[[2]]==={}&&regseq[[3]]==={}),
(*then*)
Join@@Table[
If[MemberQ[singular,regseq[[i]]=!={}&&regseq[[i,-1]]],
{#[[2]],#[[1]],1}&/@Partition[{0,Sequence@@Reverse[regseq[[i]]]},2,1],
{#[[2]],0,#[[1]]}&/@Partition[{1,Sequence@@regseq[[i]]},2,1]
],
{i,3}
],
(*if*)
MemberQ[singular,regseq[[2,-1]]],
(*then*)
Join[
{#[[2]],#[[1]],1}&/@Partition[{0,Sequence@@Reverse[Join[regseq[[1]],regseq[[2]]]]},2,1],
{#[[2]],0,#[[1]]}&/@Partition[{regseq[[1,-1]],Sequence@@regseq[[3]]},2,1]
],
(*if*)
MemberQ[singular,regseq[[3,-1]]],
(*then*)
Join[
{#[[2]],#[[1]],1}&/@Partition[{0,Sequence@@Reverse[Join[regseq[[1]],regseq[[3]]]]},2,1],
{#[[2]],0,#[[1]]}&/@Partition[{regseq[[1,-1]],Sequence@@regseq[[2]]},2,1]
],
(*default*)
True,
(*then*)
Join[
{#[[2]],0,#[[1]]}&/@Partition[{1,Sequence@@Join[regseq[[1]],regseq[[2]]]},2,1],
{#[[2]],0,#[[1]]}&/@Partition[{regseq[[1,-1]],Sequence@@regseq[[3]]},2,1]
]
];
If[Length[regseq]≥3&&!FreeQ[regseq[[3;;]],δ],Print[{in,δ}];Print["hplISelectIntOrderError"];Abort[]];
pair=SortBy[#,Position[ret,#][[1,1]]&]&/@{
If[Length[regseq[[1]]]>1&&Length[regseq[[2]]]==1,{regseq[[1,-1]],regseq[[2,1]]},Unevaluated[Sequence[]]],
If[Length[regseq[[1]]]>1&&Length[regseq[[3]]]==1,{regseq[[1,-1]],regseq[[3,1]]},Unevaluated[Sequence[]]]
};
Scan[(
tmp1=Cases[ret,{#[[2]],_,_}][[1]];
ret=DeleteCases[ret,tmp1];
tmp2=Position[ret,{#[[1]],_,_},{1},1]+{1};
ret=Insert[ret,tmp1,tmp2]
)&,pair
];
If[Variables[Fold[#1/.#2[[1]]->#2[[2;;]]&,Variables[ret],Reverse[ret]]]=!={},
Print[{in,δ}];Print["hplISelectIntOrderError"];Abort[]];
If[pair=!={},pair={Union@@pair}];
If[Length[regseq[[2]]]>1,AppendTo[pair,regseq[[2,-2;;]]]];
If[Length[regseq[[3]]]>1,AppendTo[pair,regseq[[3,-2;;]]]];
ret=F/@ret;
Scan[(
tmp1=Position[ret,F[{Alternatives@@#,_,_}]];
If[tmp1[[1;;,1]]=!=Range[tmp1[[1,1]],tmp1[[-1,1]]],Print[{in,δ}];Print["hplISelectIntOrderError"];Abort[]];
tmp2=Extract[ret,tmp1];
ret=Delete[ret,tmp1];
ret=Insert[ret,Join@@tmp2,tmp1[[1,1]]]
)&,pair
];
ret/.F->List
];

Clear[$zzz];
intToZSum[in_int]:=
Module[{res,tmp,bef,aft,len,i,j,coeff,d,tmp2,var,lg},
ifMatchQ[in[[1]],a_F b_.,
coeff=a;tmp=hplI[b,Sequence@@in[[2;;-2]]],
coeff=1;tmp=hplI@@Most[in]
];
res=$zzz[tmp];
If[res[[0]]===$zzz,
$zzz[tmp]=res=Block[{$rnum=0},
len=Length[in]-1;
Do[
While[
bef=Cases[tmp,a_hplI/;Length[a]==i,{0,∞}]//Union;
aft=$intRegularize[#,j,Last[in]]&/@bef;
bef=!=aft,
(* do *)
tmp=tmp/.Thread[bef->aft];
tmp=tmp/.t:(1+a_._Symbol)/;Abs[a]==1:>G[t];
tmp=Collect[tmp,_hplI,Together];
tmp=tmp/.G->Identity;
tmp=tmp//.{x_.(((a_Symbol-1))^(c_.)):>(x(((-1))^c))(((1-a))^c),(((a_Symbol-b_Symbol))^(c_))(((b_-a_))^(k_.))/;IntegerQ[k]:>(((a-b))^(c+k))(((-1))^k)};
],
{i,len,2,-1},{j,i-1,1,-1}
];
tmp
];
];
res=res/.hplI[a_,b___]:>hplI[a coeff,b]/.F->Identity;
res=res/.t:(1+a_._Symbol)/;Abs[a]==1:>G[t];
res=res/.a_ hplI[b_,c___]:>hplI[a b,c];
res=res/.hplI[a_Times,b___]:>hplI[F[Select[a,!FreeQ[#,Alternatives@@in[[-1]]]&]]*
Select[a,FreeQ[#,Alternatives@@in[[-1]]]&],b];
bef=Cases[res,_F,{0,∞}];
aft=bef/.F[t_]:>multipleZeta`private`seriesx[t,Reverse@in[[-1]]];
aft=aft//FunctionExpand//PowerExpand;
res=res/.Thread[bef->aft];
res=res/.Log[G[a_]]:>Log[a]/.Log[a_-b_]:>Log[a]+Log[1-b/a];
res=res//.{x_.(((a_-1))^(c_.)):>(x(((-1))^c))(((1-a))^c),(((a_Symbol-b_Symbol))^(c_))(((b_-a_))^(k_.))/;IntegerQ[k]:>(((a-b))^(c+k))(((-1))^k)};
res=res/.p:(((Log[1+n_.x_Symbol/;NumericQ[n]])^(_.))):>hplLogToHPL[p,x];
res=res//hplIExpand//hplShuffleAll;
res=res/.hplI[a_,b___]:>hplI[a/.G->Identity,b]/.{hplI[1]->1,hplI[0,___]->0}//hplIExpand;
bef=Cases[res,_hplI,{0,∞}]//Union;
aft=Switch[#,hplI[_,{_,0,1}],ts3P[#],hplI[_,_],ts1P[#],_,ts4P[#]]&/@bef;
res=res/.Thread[bef->aft];
res=res//.{x_.(((a_-1))^(c_.)):>(x(((-1))^c))(((1-a))^c),(((a_Symbol-b_Symbol))^(c_))(((b_-a_))^(k_.))/;IntegerQ[k]:>(((a-b))^(c+k))(((-1))^k)};
res=Expand[res]/.Log[a_]:>Log[a/.G->Identity]/.((G[a_])^(b_)):>G[(a^b)]//.G[a_]G[b_]:>G[a b];
var=Variables[List@@in[[2;;-2]]];
var=Union[var,Cases[res,Log[a_]/;!FreeQ[a,Alternatives@@var],{0,∞}]];
res=Collect[res,{_G,Alternatives@@var,_hplHC},F[Together[#]]&]//Expand;
res

];

greEvaluateSCL[in:greSC[_greL,___],δ_:δ,Δ_:Δ]:=
Module[{ibr,jl1,vt},
vt=If[Length[in]==1,in[[1]],SeriesCoefficient@@in];
vt=vt/.a:(_greL|Derivative[__][greL][__]):>greLToSCInt[a,$intPID];
If[!FreeQ[vt,hplSCInt],(*Print["A mode"];*)Return[vt/.a_hplSCInt:>hplSCIntToMZV[a,δ]/.hplδ->0]];
(*Print["B mode"];*)
ibr=greSCLToInt[in,δ,Δ];
ibr=If[ibr[[0]]===Plus,List@@ibr,{ibr}];
jl1=Sum[
jl1=ibr[[i]];
ClearSystemCache[];Share[];
While[
!FreeQ[jl1,int],
jl1=intExpand[jl1/.a_int/;FreeQ[List@@a,int]:>intToZSum[a/.G->Identity]];
ClearSystemCache[];Share[];
];
jl1=jl1/.F[a_]:>PowerExpand[FunctionExpand[a]]//Expand;
ClearSystemCache[];Share[];
jl1,
{i,Length[ibr]}
]/.Δ[0]->δ;
jl1//hplIntToSum//hplMZVSimplify//FunctionExpand//PowerExpand//Expand
];

greEvaluateSCL[in_greL,δ_:δ,Δ_:Δ]:=greEvaluateSCL[greSC[in],δ,Δ];

greMakePIDFunc[grePermId[2,{0,1,1,1}],$intPID]
greMakePIDFunc[grePermId[2,{0,1,1,2}],$intPID]
greMakePIDFunc[grePermId[2,{0,1,0,1}],$intPID]
greMakePIDFunc[grePermId[2,{0,1,0,2}],$intPID]
greMakePIDFunc[grePermId[2,{0,0,1,1}],$intPID]
greMakePIDFunc[grePermId[2,{0,0,1,2}],$intPID]
greMakePIDFunc[grePermId[2,{0,0,0,1}],$intPID]
greMakePIDFunc[grePermId[2,{0,0,0,2}],$intPID]
greMakePIDFunc[grePermId[1,{0,1}],$intPID]
greMakePIDFunc[grePermId[1,{0,2}],$intPID]
greMakePIDFunc[grePermId[1,{0,3}],$intPID]

Protect["gammaReduce`*"];

Protect["gammaReduce`private`*"];

Unprotect["gammaReduce`$*"];

Unprotect["gammaReduce`private`$*"];

End[];

EndPackage[];