(* :Name: StartUp`HypergeometricPFQ` *)

(* :Title: Generalized Hypergeometric Functions *)

(* :Authors: Jerry B. Keiper  and Victor S. Adamchik *)

(* :Summary:
	Provides rules for evaluating generalized hypergeometric and
	regularized hypergeometric functions.
*)

(* :Context: System` *)

(* :Package Version: 2.0 *)

(* :Copyright: Copyright 1991  Wolfram Research, Inc.

	Permission is hereby granted to modify and/or make copies of
	this file for any purpose other than direct profit, or as part
	of a commercial product, provided this copyright notice is left
	intact.  Sale, other than for the cost of media, is prohibited.

	Permission is hereby granted to reproduce part or all of
	this file, provided that the source is acknowledged.
*)

(* :History:
	Version 2.0 by Jerry B. Keiper, November 1990.

        Extensively modified and extanded by 
                       Victor S. Adamchik, February 1991.
*)

(* :Keywords: hypergeometric, regularized *)

(* :Source:
	A.P.Prudnikov, Yu.A.Brychkov, O.I.Marichev, Integrals and Series, 
        Vol. 3: More Special Functions, Gordon and Breach, New York, 
        London, 1989.
*)

(* :Mathematica Version: 2.0 *)

(* :Limitation:
	This package provides only the obvious simplifications.  There
	are many specialized identities that could be added.  For example,
	the generalized Euler transformation could be used to numerically
	evaluate pFq where p-q == 1 and Abs[z-1] > 1.
*)

Needs["System`", "StartUp`Attributes.m"]

Begin["System`"]

Unprotect[ HypergeometricPFQ, HypergeometricPFQRegularized, 
           Hypergeometric1F1, Hypergeometric2F1]

Map[ Clear, {HypergeometricPFQ, HypergeometricPFQRegularized,
             Hypergeometric1F1, Hypergeometric2F1}]

HypergeometricPFQ::hdiv = 
" Warning: Divergent generalized hypergeometric series `1`."

HypergeometricPFQ::usage =
"HypergeometricPFQ[numlist, denlist, z] gives the generalized hypergeometric
function pFq where numlist is a list of the p parameters in the numerator 
and denlist is a list of the q parameters in the denominator."

HypergeometricPFQRegularized::usage =
"HypergeometricPFQRegularized[numlist, denlist, z] gives the regularized
generalized hypergeometric function pFq where numlist is a list of the p
parameters in the numerator and q is a list of the q parameters in the
denominator.  HypergeometricPFQRegularized[numlist, denlist, z] ==
HypergeometricPFQ[numlist, denlist, z]/(Times @@ (Gamma /@ denlist))
except when denlist contains a non-positive integer when it is defined by
analytic continuation to remove the indeterminacy."

Begin["HypergeometricPFQ`Private`"]

(***************************************************************************
*                           General Cases
*
****************************************************************************)
 
 Unprotect[ Hypergeometric2F1, Hypergeometric1F1 ]

 HypergeometricPFQ[ uppar_,lowpar_,0 ] :=  1            

 HypergeometricPFQ[ uppar_,lowpar_,z_Real/;z==0. ] := 1 + z           

 HypergeometricPFQ[ uppar_, {v1___,u_Integer/;u<=0,v2___},arg_ ] :=
   Module[ { negnum, negden },
      negnum = Select[uppar,IntegerQ[#] && #<=0&];
      If[ Length[negnum]==0, 
	  If[ And@@(NumberQ[#]&/@Re[N[uppar]]), 
	      ComplexInfinity, 
	      Indeterminate 
	  ],
          negden = Select[{v1,u,v2},IntegerQ[#] && #<=0&];
	  If[ And@@(NumberQ[#]&/@Re[N[Join[uppar,{v1,v2}]]]),
              If[ Max[negden] > Max[negnum], 
		  ComplexInfinity,
                  Indeterminate
              ],
	     Indeterminate
	  ]
     ]
  ]         

 HypergeometricPFQ[ {v1___,a_,v2___},{v3___,b_,v4___},arg_] :=
   HypergeometricPFQ[ {v1,v2},{v3,v4}, arg ] /; N[Expand[a-b]] == 0.

 HypergeometricPFQ[ {___,0,___},lowpar_,arg_ ] :=  1  
          
 HypergeometricPFQ[ {___,z_Real/;z==0.,___},lowpar_,arg_ ] := 1 + z            


 HypergeometricPFQ[ {v1___,u_Integer?Negative,v2___},lowpar_,arg_ ] :=
   Sum[ arg^i MultPochham[{u,v1,v2},i]/( i! MultPochham[lowpar,i]),
       {i,0,-Max[Select[ {u,v1,v2},IntegerQ[#] && Negative[#]& ]]} ]
   
 HypergeometricPFQ[ uppar_List,lowpar_List,arg_ ] := 
   (
   Message[HypergeometricPFQ::hdiv,
           HoldForm[HypergeometricPFQ[ uppar,lowpar,arg]]];
   HypergeometricPFQ[ uppar,lowpar,arg] /; Fail
   ) /;
 Length[uppar]>Length[lowpar]+1
   
 HypergeometricPFQ[ uppar_List,lowpar_List,arg_ ] := 
    Block[ {inter, HypergeometricCond},
        Off[HypergeometricPFQ::hdiv];
        inter = Expand[
        GeneralHypergeometricPFQ[ Sort[Expand[uppar]],
                                  Sort[Expand[lowpar]],arg ]//.
        PolyGammaRule];
        On[HypergeometricPFQ::hdiv];
        inter = 
        If[ !FreeQ[inter,Gamma] , 
               SimpGammaH[inter],
               inter ];
        inter = 
        If[ !FreeQ[inter,PolyGamma] , 
            SimpPolyGammaH[inter],
            inter ];
        If[ Or@@(Not[FreeQ[inter,#]]&/@{Tan,Cot,Tanh,Coth}),
            inter//.SimpTrigSum ,
            inter]
      ] /;  Accuracy[{uppar,lowpar,arg}] === Infinity && HypergeometricCond 

 HypergeometricCond = True
      
(* Approximate Numerical Evaluation *)
HypergeometricPFQ[a_List, b_List, z_] :=
    Module[{an = a, bn = Append[b, 1], oldsum = 0, sum = 1, term = 1},
        While[oldsum != sum,
            term *= Apply[Times, an++] z/Apply[Times, bn++];
            oldsum = sum;
            sum += term];
        sum
    ] /; (Apply[And, Map[NumberQ, a]] && Apply[And, Map[NumberQ, b]] &&
		NumberQ[z] && Precision[{a, b, z}] < Infinity &&
		(Length[a] - Length[b] < 1 ||
			(Length[a] - Length[b] == 1 && Abs[z] < 1)))
 
(* Derivatives wrt z *)
Derivative[0, 0, n_Integer][HypergeometricPFQ] ^:=
    (Apply[Times, Map[Pochhammer[#, n]&, #1]] *
        HypergeometricPFQ[#1 + n, #2 + n, #3] /
        Apply[Times, Map[Pochhammer[#, n]&, #2]])& /; n > 0
 
 GeneralHypergeometricPFQ[ {___,0,___},lowpar_,arg_ ] :=  1            


 GeneralHypergeometricPFQ[ {v1___,u_Integer?Negative,v2___},
                           lowpar_,arg_ ] :=
 Sum[ arg^i MultPochham[{u,v1,v2},i]/( i! MultPochham[lowpar,i]),
      {i,0,-Max[Select[ {u,v1,v2},IntegerQ[#] && Negative[#]& ]]} ]

 GeneralHypergeometricPFQ[ {v1___,a_,v2___},{v3___,b_,v4___},arg_] :=
   GeneralHypergeometricPFQ[ {v1,v2},{v3,v4}, arg ] /; Expand[a-b] === 0
   
 GeneralHypergeometricPFQ[ {},{a_/;!NumberQ[a]||Denominator[a]!=2},
                         arg_/;Znak[arg] ] :=
   Gamma[a] (-arg)^(-a/2+1/2) BesselJ[Expand[a-1],2 Sqrt[-arg]]

 GeneralHypergeometricPFQ[ {},{a_/;!NumberQ[a]||Denominator[a]!=2},
                         arg_ ] :=
   Gamma[a] arg^(-a/2+1/2) BesselI[Expand[a-1],2 Sqrt[arg]]
   
 GeneralHypergeometricPFQ[ {},{a_},arg_ ] :=
   Hypergeometric0F1[ a,arg ] 

 GeneralHypergeometricPFQ[ {},{},arg_ ] := Exp[arg] 

 GeneralHypergeometricPFQ[ {par_},{},arg_/;arg=!=1 ] := 
   Together[1 - arg]^(-par) 

 GeneralHypergeometricPFQ[ {par_},{},1 ] := 0 /;
 NumberQ[par] && Im[par]==0 && par<0

 GeneralHypergeometricPFQ[ {par_},{},1 ] := 
   DirectedInfinity[] /; NumberQ[par] && Im[par]==0 && par>0

 GeneralHypergeometricPFQ[ {c1___,a_,c2___}, {d1___,a1_,d2___},arg_ ] :=
   GeneralHypergeometricPFQ[ {a-1,c1,c2},{a1,d1,d2},arg] +                             
   arg Apply[Times,{c1,c2}]/( a1 Apply[Times,{d1,d2}])*
   GeneralHypergeometricPFQ[ 1+{a-1,c1,c2},1+{a1,d1,d2},arg ] /;                             
 IntegerQ[Expand[a-a1]] && Expand[a-a1]>0
 
 GeneralHypergeometricPFQ[ {c1___,a_Integer,c2___},
                           {d1___,a1_Integer,d2___},arg_ ] :=
   (a1-1) Apply[Times,{d1,d2}-1] / (arg Apply[Times,{c1,c2}-1] ) *
   (GeneralHypergeometricPFQ[
           Expand[{a+1,c1,c2}-1],Expand[{a1,d1,d2}-1],arg] - 
    GeneralHypergeometricPFQ[
           Expand[{a,c1,c2}-1],Expand[{a1,d1,d2}-1],arg] ) /;a<a1 && 
 Length[Intersection[{c1,c2},{1}]] == 0

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] :=
   Module[ {a,a1,up,low},
     {a,up,low}  = SearchPar[0,uppar,lowpar];
     If[Not[FreeQ[a,HyperFail]], a1 = HyperFail,
        {a1,up,low} = SearchPar[a[[1]],up,low] ];
     EvalHyper[a,a1,up,low,arg] /;
   a1=!=HyperFail && ConditionHyp[a,a1,up,low,arg]
   ]/;Length[lowpar]>1

 EvalHyper[{a_,a1_},{b_,b1_},uppar_,lowpar_,arg_] :=
    b GeneralHypergeometricPFQ[Join[{a,b+1},uppar],
                        Join[{a1,b1},lowpar],arg ]/(b-a) -
    a GeneralHypergeometricPFQ[Join[{a+1,b},uppar],
                        Join[{a1,b1},lowpar],arg ]/(b-a)
  
 SearchPar[t_,{c1___,a_,c2___},{d1___,a1_,d2___}] := 
    {{a,a1},{c1,c2},{d1,d2}}/;a=!=t && IntegerQ[a1-a] && a1-a>0

 SearchPar[___] := {HyperFail,HyperFail,HyperFail}

 ConditionHyp[{a_,a1_},{b_,b1_},c_,d_,arg_] :=
   If[ Length[c] == 1 && Length[d] == 0 && N[arg]=!=1.||
       Length[c] == 2 && Length[d] == 1 ||
       Length[c] <= Length[d] ||
       N[Abs[arg]] < 1. || Not[NumberQ[N[arg]]] ||
       (Re[Expand[Apply[Plus,Join[c,-d,{a,b+1,-a1,-b1}] ]]] < 0 &&
        Re[Expand[Apply[Plus,Join[c,-d,{a+1,b,-a1,-b1}] ]]] < 0 &&
        Abs[arg] == 1) ||
       (Re[Expand[Apply[Plus,Join[c,-d,{a,b+1,-a1,-b1}] ]]] < 1 &&
        Re[Expand[Apply[Plus,Join[c,-d,{a+1,b,-a1,-b1}] ]]] < 1 &&
        Abs[arg] == 1 && arg=!=1),
     True,
     False]

(***************************************************************************
*                         Particular Cases
*
***************************************************************************)

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
   Module[ {answer = HypergeometricMid[ uppar,lowpar,arg ]},
     answer/;FreeQ[answer,HyperFail]
  ]/;
 Length[uppar] == Length[lowpar] + 1 

 HypergeometricMid[ {v1___,b_,v2___},u_,arg_ ] :=
   Module[ {answer = FormulaS[ Length[u],b,u[[1]]-1,arg]},
     answer/;FreeQ[answer,HyperFail]
   ]/;
   {v1,v2} === u-1 && Length[Union[u]] == 1

 FormulaS[ p_/;p=!=1,1,a_,1 ] := (-a)^p PolyGamma[p-1,a]/(p-1)!

 FormulaS[ p_,1,a_,-1 ] :=
      (-a/2)^p (PolyGamma[p-1,Expand[a/2]] - 
                PolyGamma[p-1,Expand[a/2+1/2]])/(p-1)! 
                                                                     
 FormulaS[ p_/;p=!=1,1,a_,z_ ] := a^p LerchPhi[z,p,a]

 FormulaS[ p_,b_,a_,1 ] :=
  Module[{r},
     -Gamma[1-b] (-a)^p/(p-1)! *
       (D[Gamma[r]/Gamma[1+r-b],{r,p-1}]//.{r->a}) ]/;
 Not[IntegerQ[b]]  

 FormulaS[ p_,b_Integer,a_,-1 ] :=
  Module[{r,k},
     -(-a)^p/(Gamma[p] Gamma[b])*
     (D[ Sum[Gamma[k] Pochhammer[1+k-r,b-1-k]/2^k,{k,1,b-1} ] +
         Gamma[b-r] (PolyGamma[0,Expand[r/2+1/2]]-
         PolyGamma[0,Expand[r/2]])/(2 Gamma[1-r]),{r,p-1} ]//.{r->a})]

 FormulaS[__] := HyperFail

(********************* Formula 9, p. 572, PBM *****************************)

 HypergeometricMid[ v_,{u1___,u2__},arg_ ] := 
   FormulaForPolylog[Length[v]-1,Length[{u2}],arg]/;
 Union[v]==={1} &&(Union[{u1}]==={2} && Union[{u2}]==={3}||
                   Length[{u1}]==0 && Union[{u2}]==={3} )

 FormulaForPolylog[ q_,n_,arg_ ] :=
    (-1)^q 2^n/arg Expand[
     Sum[Pochhammer[q,k]/k! PolyLog1[n-k,arg],{k,0,n-1}]/arg +
     Sum[Pochhammer[n,k] (-1)^(q-k)/k! PolyLog1[q-k,arg],{k,0,q-1}] -
     Binomial[n+q-1,q] ]/.PolyLog1->PolyLog


 HypergeometricMid[ __ ] := HyperFail

(**************************************************************************
*                      HyperBessel Functions
*
***************************************************************************)

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
   Module[ {answer =
              If[Not[Znak[arg]],
                 F03plus[lowpar,4 arg^(1/4)],
                 F03minus[lowpar,2 (-4 arg)^(1/4)]] },
     answer/;FreeQ[answer,HyperFail]
  ]/;
  Length[uppar]==0 && Length[lowpar]==3 

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
   Module[ {
     answer =
         F02p611[lowpar,
         PowerExpand[3 If[Znak[arg],-(-arg)^(1/3),arg^(1/3)]]] },
     answer/;FreeQ[answer,HyperFail]
   ]/;
  Length[uppar]==0 && Length[lowpar]==2 

(************************* Formula 7 & 8, p. 612, PBM *********************)

 F03plus[{v___,1/2,u___},z_] := HFormula7[{v,u},z]

 F03plus[{v___,3/2,u___},z_] := HFormula8[{v,u},z]

 F03plus[ __ ] := HyperFail

 HFormula7[ {a_,b_},z_ ] := 
    Gamma[Expand[2 a]] (z/4)^Expand[1-a 2] (BesselI[Expand[2 a-1],z] +
                   BesselJ[Expand[2 a-1],z])/(4^a)/;
 Expand[b-a]===1/2

 HFormula7[ {b_,a_},z_ ] := 
    Gamma[Expand[2 a]] (z/4)^Expand[1-a 2] (BesselI[Expand[2 a-1],z] +
                   BesselJ[Expand[2 a-1],z])/(4^a)/;
 Expand[b-a]===1/2

 HFormula8[ {a_,b_},z_ ] := 
    Gamma[Expand[2 a]] (BesselI[Expand[2 a-2],z] -
            BesselJ[Expand[2 a-2],z])/(2^(2 a+1) (z/4)^Expand[a 2])/;
 Expand[b-a]===1/2

 HFormula8[ {b_,a_},z_ ] := 
    Gamma[Expand[2 a]] (BesselI[Expand[2 a-2],z] -
            BesselJ[Expand[2 a-2],z])/(2^(2 a+1) (z/4)^Expand[a 2])/;
 Expand[b-a]===1/2

 HFormula7[ __ ] := HyperFail

 HFormula8[ __ ] := HyperFail

(************************* Formula 12, p. 612, PBM **************************)

 F03minus[ {v___,1/2,u___},z_ ] := HFormula1214[ {v,u},z ]

 F03minus[ {v___,3/2,u___},z_ ] := HFormula1619[ {v,u},z ]

 F03minus[ __ ] := HyperFail

 HFormula1214[ {1/4,3/4},z_ ] := Cosh[z] Cos[z]
 
 HFormula1214[ {3/4,1/4},z_ ] := Cosh[z] Cos[z]

 HFormula1214[ {5/4,3/4},z_ ] := (Sinh[z] Cos[z] + 
       Cosh[z] Sin[z])/(2 z)

 HFormula1214[ {3/4,5/4},z_ ] := (Sinh[z] Cos[z] + 
       Cosh[z] Sin[z])/(2 z)

 HFormula1214[ __ ] := HyperFail

 HFormula1619[ {5/4,3/4},z_ ] := Sinh[z] Sin[z] z^(-2)

 HFormula1619[ {3/4,5/4},z_ ] := Sinh[z] Sin[z] z^(-2)

 HFormula1619[ {5/4,7/4},z_ ] := (Cosh[z] Sin[z] - 
       Sinh[z] Cos[z]) z^(-3) 3/2

 HFormula1619[ {7/4,5/4},z_ ] := (Cosh[z] Sin[z] - 
       Sinh[z] Cos[z]) z^(-3) 3/2

 HFormula1619[ __ ] := HyperFail

(************************* Formula 1-3, p. 611, PBM *************************)

 F02p611[ {1/3,2/3},z_ ] := (E^z+2 E^(-z/2) Cos[z/2 Sqrt[3]])/3

 F02p611[ {2/3,4/3},z_ ] := (E^z-2 E^(-z/2) Cos[z/2 Sqrt[3]+Pi/3])/(3 z)

 F02p611[ {4/3,5/3},z_ ] := 2 (E^z-2 E^(-z/2) Cos[z/2 Sqrt[3]-Pi/3])/(3 z^2)

 F02p611[ __ ] := HyperFail
 
(**************************************************************************** 
*                         Hypergeometric3F2 
*
*****************************************************************************)

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
   Module[ { answer = F32[uppar,lowpar,arg] },
     answer/;FreeQ[answer,HyperFail]
   ]/;
  Length[uppar] == 3 && Length[lowpar] == 2

 F32[ {a___,1,b___},{v1_,v2_},1] := 
  (v1-1)(v2-1)(PolyGamma[0,v1-1] - PolyGamma[0,v2-1])/(v1-v2)/;
    v1=!=v2 && ({a,b}==={v1-1,v2-1} || {b,a}==={v1-1,v2-1})

 F32[ {1/2,1,1},{3/2,3/2},1 ] := 2 Catalan

 F32[ {1/2,1/2,1},{5/2,5/2},-1 ] := 9/8(4-Pi)

(************************* Formula 13, p. 498, PBM **************************)

 F32[ {v1___,a_,v2___,a1_,v3___},{u1___,a2_,u2___},z_ ] := 
         F32N1[{a,v1,v2,v3},{u1,u2},z]/;   
  Expand[a1-a-1/2]===0 && Expand[a2- 2 a]===0

 F32N1[ {a_,b_},{c_},z_ ] := ( (1+Sqrt[1-z])/2 )^(-b)*
         F21[b,2 a-b,b+1,1/2-Sqrt[1-z]/2]/;
  c-b===1
 
 F32N1[ __ ] := HyperFail

(************************* Formula 30, p. 499, PBM **************************)

 F32[ v_,{1/3,2/3},z_ ] := F32Formula30[ v,((Plus@@v)-1)/3, 
    If[ Znak[z],(-z)^(1/3) E^(Pi I/3),z^(1/3)] ]
      
 F32Formula30[ v_,midl_,arg_ ] := ( (1-arg)^(-3 midl) +
    (1-arg E^(2Pi I/3))^(-3 midl) + (1-arg E^(4Pi I/3))^(-3 midl) )/3
 Union[Expand[v - Sort[{midl,midl+1/3,midl+2/3}]]]==={0}
      
 F32Formula30[ __ ] := HyperFail

(************************* Formula 31, p. 499, PBM **************************)

 F32[ v_,{2/3,4/3},z_ ] := F32Formula31[ v,((Plus@@v)-1)/3, 
    If[ Znak[z],(-z)^(1/3) E^(Pi I/3),z^(1/3)] ]
      
 F32Formula31[ v_,midl_,arg_ ] := ( (1-arg)^(1-3 midl) +
    E^(-2Pi I/3)(1-arg E^(2Pi I/3))^(1-3 midl) + 
    E^(-4Pi I/3)(1-arg E^(4Pi I/3))^(1-3 midl) )/
    (3 arg (3 midl-1))
 Union[Expand[v - Sort[{midl,midl+1/3,midl+2/3}]]]==={0}
      
 F32Formula31[ __ ] := HyperFail

(************************* Formula 24, p. 535, PBM **************************)

 F32[ {v1___,a_,v2___,a1_,v3___},{u1_,u2_},1 ] := 
         F32N2[{a,v1,v2,v3},{u1,u2}]/;   
  Expand[a+a1-1]===0 

 F32N2[ {a_,b_},{c_,d_} ] := 2^(1-2 b) Pi Gamma[c] Gamma[d]/
      (Gamma[Expand[a/2+c/2]] Gamma[Expand[b+(1+a-c)/2]]*
       Gamma[Expand[(1+c-a)/2]] Gamma[Expand[1+b-(a+c)/2]])/;
   Expand[d-1 +c-2 b]===0

 F32N2[ {a_,b_},{d_,c_} ] := 2^(1-2 b) Pi Gamma[c] Gamma[d]/
      (Gamma[Expand[a/2+c/2]] Gamma[Expand[b+(1+a-c)/2]]*
       Gamma[Expand[(1+c-a)/2]] Gamma[Expand[1+b-(a+c)/2]])/;
   Expand[d-1 +c-2 b]===0
 
 F32N2[ __ ] := HyperFail

(************************* Formula 40-41, p. 537, PBM ***********************)
  
 F32[ {v1___,a_/;a=!=1,v2___},{u1___,b_,u2___},1 ] :=
   (b-1) (PolyGamma[0,b-1]-PolyGamma[0,b-a])/(a-1)/;
 {v1,v2}==={1,1} && {u1,u2}==={2}
  
 F32[ {1,1,1},{u1___,b_,u2___},1 ] :=
   (b-1) PolyGamma[1,b-1]/;
 {u1,u2}==={2}

(************************* Formula 42-44, p. 537, PBM ***********************)
  
 F32[ {v1___,a_/;a=!=1&&a=!=2,v2___},{u1___,b_,u2___},1 ] :=
   2 (b-1) (b-a) (PolyGamma[0,b-a+1]-PolyGamma[0,b-1])/((a-2)(a-1)) +
   2 (b-1)/(a-1)/;
 {v1,v2}==={1,1} && {u1,u2}==={3}
  
 F32[ {1,1,1},{u1___,b_,u2___},1 ] :=
   2 (2-b) + 2 (b-1)^2 PolyGamma[1,b]/;
 {u1,u2}==={3}
  
 F32[ {v1___,2,v2___},{u1___,b_,u2___},1 ] :=
   2 (b-1) (1 - (b-2) PolyGamma[1,b-1])/;
 {v1,v2}==={1,1} && {u1,u2}==={3}

(************************* Formula 50, 67, 70, p. 537, PBM ******************)
 
 F32[ {v1___,b_/;FreeQ[{1,2,3},b],v2___},u_,1 ] := 
    F32Formula50[Union[{v1,v2}][[1]],b]/;u===2+{v1,v2}

 F32Formula50[ a_,b_ ] := 
    Gamma[2+a] Gamma[1-b]/Gamma[2+a-b] ( a(a+1)(b-2a-1) (
    PolyGamma[0,2+a]-PolyGamma[0,2+a-b] ) +2a^2+(1-b)(1-2a^2) )
 
 F32[ {v1___,3,v2___},u_,1 ] := 
    F32Formula70[Union[{v1,v2}][[1]]]/;u===2+{v1,v2}

 F32Formula70[ a_ ] := a^2(a+1)^2/2 (3-2a+2(a-1)^2 PolyGamma[1,a])

 F32[ {v1___,2,v2___},u_,1 ] := 
    F32Formula67[Union[{v1,v2}][[1]]]/;u===2+{v1,v2}

 F32Formula67[ a_ ] := -a^2(a+1)^2 (-2+(2a-1) PolyGamma[1,a])

(************************* Formula 88, p. 539, PBM **************************)

 F32[ {v1___,n_/;Znak[n],v2___},{c_,d_},1 ] := F32Formula88[-n,{v1,v2},c]/;
 Expand[d-1+c-n-Plus@@{v1,v2}]===0

 F32[ {v1___,n_/;Znak[n],v2___},{d_,c_},1 ] := F32Formula88[-n,{v1,v2},c]/;
 Expand[d-1+c-n-Plus@@{v1,v2}]===0

 F32Formula88[ n_,{a_,b_},c_ ] := Pochhammer[c-a,n] Pochhammer[c-b,n]/(
    Pochhammer[c,n] Pochhammer[c-a-b,n])

(************************* Formula 97 & 99, p. 539, PBM ********************)

 F32[ {v1___,k_Integer?Negative a_Symbol,v2___},u_,1 ] := 
    F32Form99[ -k a,{v1,v2},u ]

 F32Form99[ n_,{a_,b_},{v1_,v2_} ] := 
   Pochhammer[1+a,n] Pochhammer[1+a/2-b,n]/(
   Pochhammer[1+a/2,n] Pochhammer[1+a-b,n] ) /; a+n=!=0 &&(
   v1===1+a+n && v2===1+a-b || v1===1+a-b && v2===1+a+n)

 F32Form99[ n_,{b_,a_},{v1_,v2_} ] := 
   Pochhammer[1+a,n] Pochhammer[1+a/2-b,n]/(
   Pochhammer[1+a/2,n] Pochhammer[1+a-b,n] ) /; a+n=!=0 &&(
   v1===1+a+n && v2===1+a-b || v1===1+a-b && v2===1+a+n)

 F32Form99[ n_,{a_,b_},{v1_,v2_} ] := (1+(-1)^n) *
   (-4)^(n/2) Pochhammer[1-a-b-n,n/2] Pochhammer[1/2,n/2]/(
   2 Pochhammer[1-a-n,n/2] Pochhammer[1-b-n,n/2]) /;
   v1===1-a-n && v2===1-b-n || v1===1-b-n && v2===1-a-n 

 F32Form99[ n_,{b_,a_},{v1_,v2_} ] := (1+(-1)^n) *
   (-4)^(n/2) Pochhammer[1-a-b-n,n/2] Pochhammer[1/2,n/2]/(
   2 Pochhammer[1-a-n,n/2] Pochhammer[1-b-n,n/2]) /;
   v1===1-a-n && v2===1-b-n || v1===1-b-n && v2===1-a-n 
 
 F32Form99[ __ ] := HyperFail

 F32[ __ ] := HyperFail

(****************************************************************************
*                         Hypergeometric2F1
*
*****************************************************************************)

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
   Hypergeometric2F1[ uppar[[1]],uppar[[2]],lowpar[[1]],arg] /;
 Length[uppar] == 2 && Length[lowpar] == 1 &&
 Accuracy[{uppar,lowpar,arg}] < Infinity

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
  Module[ {answer},
   answer = F21[uppar[[1]],uppar[[2]],lowpar[[1]],arg];
   answer = If[ FreeQ[answer,HyperFail],
                (answer/.LogTrig)/.Gamma[w_] :> Gamma[Expand[w]] ,
            If[ arg=!=1,
                (1-arg)^(lowpar[[1]]-uppar[[1]]-uppar[[2]]) *
                F21[ lowpar[[1]]-uppar[[1]],lowpar[[1]]-uppar[[2]],
                     lowpar[[1]],arg ],
                HyperFail]];
   If[ FreeQ[answer,HyperFail],
       (answer/.LogTrig)/.Gamma[w_] :> Gamma[Expand[w]] , 
       GaussFlag = False;   
       answer = Hypergeometric2F1[ uppar[[1]],uppar[[2]],
                                   lowpar[[1]],arg ]];
       GaussFlag = True;
       answer 
   ]/;
  Length[uppar] == 2 && Length[lowpar] == 1

 GaussFlag = True
(***************************** In Unit **************************************)

 F21[ a_,b_,c_,1 ] := Gamma[c] Gamma[c-a-b]/(Gamma[c-a] Gamma[c-b])/;
 Not[NumberQ[c-a-b]] || Re[c-a-b]>0

 F21[ 1,a_,b_,-1 ] := a PolyGamma[0,a/2+1/2]/2-a PolyGamma[0,a/2]/2/;
 b===a+1

 F21[ a_,1,b_,-1 ] := a PolyGamma[0,a/2+1/2]/2-a PolyGamma[0,a/2]/2/;
 b===a+1

 F21[ 1,a_,b_,-1 ] := a (a+1) (PolyGamma[0,(a+1)/2]-PolyGamma[0,a/2])-a-1/;
 b===a+2
 
 F21[ a_,1,b_,-1 ] := a (a+1) (PolyGamma[0,(a+1)/2]-PolyGamma[0,a/2])-a-1/;
 b===a+2
 
 F21[ 1,a_,b_,-1 ] := a (a+1) (a+2) (PolyGamma[0,(a+1)/2]-PolyGamma[0,a/2])-
                      a^2-7 a/2-3/;
 b===a+3
 
 F21[ a_,1,b_,-1 ] := a (a+1) (a+2) (PolyGamma[0,(a+1)/2]-PolyGamma[0,a/2])-
                      a^2-7 a/2-3/;
 b===a+3

 F21[ 2,a_,b_,-1 ] := (b-2) (b-1) ( 1-(b-1) (a+b-3) F21[1,a-1,b-1,-1] )/
                     (2 (1-a))/;
 a=!=1 

 F21[ a_,2,b_,-1 ] := (b-2) (b-1) ( 1-(b-1) (a+b-3) F21[1,a-1,b-1,-1] )/
                     (2 (1-a))/;
 a=!=1 

 F21[ 2,a_,b_,-1 ] := a (a-1) (PolyGamma[0,a/2]-PolyGamma[0,(a-1)/2])/2 -a/2 /;
 b===a+1

 F21[ a_,2,b_,-1 ] := a (a-1) (PolyGamma[0,a/2]-PolyGamma[0,(a-1)/2])/2 -a/2 /;
 b===a+1

 F21[ a_,b_,c_,-1 ] := 2^(-a) Sqrt[Pi] Gamma[1+a-b]/(Gamma[1/2+a/2]*
                  Gamma[1-b+a/2])/;
 Expand[c-1-a+b]===0

 F21[ b_,a_,c_,-1 ] := 2^(-a) Sqrt[Pi] Gamma[1+a-b]/(Gamma[1/2+a/2]*
                  Gamma[1-b+a/2])/;
 Expand[c-1-a+b]===0

 F21[ a_,b_,c_,-1 ] := 2^(-a) Sqrt[Pi] Gamma[c]*
   (1/(Gamma[a/2-b] Gamma[a/2+1/2])+a/(2 Gamma[a/2-b+1/2] Gamma[a/2+1]))/; 
 Expand[c-a+b]===0

 F21[ b_,a_,c_,-1 ] := 2^(-a) Sqrt[Pi] Gamma[c]*
   (1/(Gamma[a/2-b] Gamma[a/2+1/2])+a/(2 Gamma[a/2-b+1/2] Gamma[a/2+1]))/; 
 Expand[c-a+b]===0

 F21[ b_,a_,c_,2 ] := 2^(-2 a) (-a)!/(-2 a)! Pochhammer[1/2+b/2,-a] /; 
 Expand[c-2 a]===0 && ZnakSum[a] && ZnakSum[c]

 F21[ a_,b_,c_,2 ] := 2^(-2 a) (-a)!/(-2 a)! Pochhammer[1/2+b/2,-a] /; 
 Expand[c-2 a]===0 && ZnakSum[a] && ZnakSum[c]

(************************* Formula 106, p. 461, PBM **************************)

 F21[ a_,b_,1/2,z_/;Not[Znak[z]] ] := 
   ((1+Sqrt[z])^(-2 a)+(1-Sqrt[z])^(-2 a))/2/;Expand[b-a-1/2]===0

 F21[ b_,a_,1/2,z_/;Not[Znak[z]] ] := 
   ((1+Sqrt[z])^(-2 a)+(1-Sqrt[z])^(-2 a))/2/;Expand[b-a-1/2]===0

(************************* Formula 107, p. 461, PBM **************************)

 F21[ a_,b_,3/2,z_/;Not[Znak[z]] ] := 
   (-(1+Sqrt[z])^(1-2 a)+(1-Sqrt[z])^(1-2 a))/(2 (2 a-1) Sqrt[z])/;
 Expand[b-a-1/2]===0 && a=!=1/2

 F21[ b_,a_,3/2,z_/;Not[Znak[z]] ] := 
   (-(1+Sqrt[z])^(1-2 a)+(1-Sqrt[z])^(1-2 a))/(2 (2 a-1) Sqrt[z])/;
 Expand[b-a-1/2]===0 && a=!=1/2

(************************* Formula 132, p. 463, PBM **************************)

 F21[ 1,n_Integer,b_Rational?Positive,z_/;Abs[z]=!=1 ] := 
    Module[ {m},
     m = n-b+1/2;
      Pochhammer[1/2,m] F21[1,b-1/2,b,z]/((1-z)^m Pochhammer[b-1/2,m]) +
      (b-1) Sum[ Pochhammer[1/2-m,k]/(Pochhammer[2-n,k] (1-z)^(1+k)),
      {k,0,m-1}]/(n-1) ]/;
 b!=1/2 && Denominator[b]===2 && n-b-1/2 >= 0

 F21[ n_Integer,1,b_Rational?Positive,z_/;Abs[z]=!=1 ] := 
    Module[ {m},
     m = n-b+1/2;
      Pochhammer[1/2,m] F21[1,b-1/2,b,z]/((1-z)^m Pochhammer[b-1/2,m]) +
      (b-1) Sum[ Pochhammer[1/2-m,k]/(Pochhammer[2-n,k] (1-z)^(1+k)),
      {k,0,m-1}]/(n-1) ]/;
 b!=1/2 && Denominator[b]===2 && n-b-1/2 >= 0

(************************* Formula 133, p. 463, PBM **************************)

 F21[ 1,n_Integer,b_Rational?Positive,z_/;Abs[z]=!=1 ] := 
    Pochhammer[1/2,n] (2 Sqrt[z]*
    ArcSin[Sqrt[z]]/Sqrt[1-z] - Sum[(k-1)! z^k/Pochhammer[1/2,k],
        {k,1,n-1}])/(z^n (n-1)!)/;
 b-n===1/2 && Not[Znak[z]]

 F21[ n_Integer,1,b_Rational?Positive,z_/;Abs[z]=!=1 ] := 
    Pochhammer[1/2,n] (2 Sqrt[z]*
    ArcSin[Sqrt[z]]/Sqrt[1-z] - Sum[(k-1)! z^k/Pochhammer[1/2,k],
         {k,1,n-1}])/(z^n (n-1)!)/;
 b-n===1/2 && Not[Znak[z]]

(************************* Formula 134, p. 463, PBM **************************)

 F21[ 1,a_Rational?Positive,b_,z_/;Abs[z]=!=1 ] := 
   2 a z^(-a) ArcTanh[ Sqrt[z] ] -
   2 a z^(-a-1/2) Sum[z^k/(2 k-1),{k,1,a-1/2}]/;
 Denominator[a]==2 && b===a+1 && Not[Znak[z]]

 F21[ a_Rational?Positive,1,b_,z_/;Abs[z]=!=1 ] := 
   2 a z^(-a) ArcTanh[ Sqrt[z] ] -
   2 a z^(-a-1/2) Sum[z^k/(2 k-1),{k,1,a-1/2}]/;
 Denominator[a]==2 && b===a+1 && Not[Znak[z]]

(************************* Formula 135, p. 463, PBM **************************)

 F21[ 1,n_Integer,m_,z_/;Abs[z]=!=1 ] := 
    -n Log[1-z]/z^n - n z^(-n) Sum[z^k/k,{k,1,n-1}]/; m===n+1

 F21[ n_Integer,1,m_,z_/;Abs[z]=!=1 ] := 
    -n Log[1-z]/z^n - n z^(-n) Sum[z^k/k,{k,1,n-1}]/;m===n+1 

(************************* Formula 136, p. 463, PBM **************************)

 F21[ 1,1,m_Integer?Positive,z_/;z=!=1 ] := 
    (m-1) z/(z-1)^2( Sum[(z-1)^k z^(-k)/(m-k),{k,2,m-1}] -
                     ((z-1)/z)^m Log[1-z] ) 

(************************* Formula 137, p. 463, PBM **************************)

 F21[ v1_,v2_,a_,z_/;Abs[z]=!=1 && Not[Znak[z]] ] := 
    Pochhammer[1/2,a-1/2] (z-1)^(a-3/2) ( 2 Sqrt[z] *
    ArcTanh[ Sqrt[z] ] + Sum[(k-1)! z^k/(Pochhammer[1/2,k] *
    (z-1)^k),{k,1,a-3/2}])/((a-3/2)! z^(a-1/2))/;
 IntegerQ[a-1/2] && a-1/2 > 0 && (v1===1 && v2===1/2 || v2===1 && v1===1/2)

 F21[ v1_,v2_,a_,z_/;Abs[z]=!=1 && Znak[z] ] := 
    Pochhammer[1/2,a-1/2] (z-1)^(a-3/2) ( Sqrt[-z] *
    (-2) ArcTan[Sqrt[-z]] + Sum[(k-1)! z^k/(Pochhammer[1/2,k] *
    (z-1)^k),{k,1,a-3/2}])/((a-3/2)! z^(a-1/2))/;
 IntegerQ[a-1/2] && a-1/2 > 0 && (v1===1 && v2===1/2 || v2===1 && v1===1/2)

(************************* Formula 138, p. 464, PBM **************************)

 F21[ 1,1,b_,z_/;Abs[z]=!=1 && Not[Znak[z]] ] := 
    (2 b - 2) (z-1)^(b-3/2) z^(1/2-b) ( 
       z^(1/2) (1-z)^(-1/2) ArcSin[z^(1/2)] + 
       Sum[ z^k (z-1)^(-k) (2 k-1)^(-1),{k,1,b-3/2}])/;
 IntegerQ[b-1/2] && b>1  

 F21[ 1,1,b_,z_/;Abs[z]=!=1 && Znak[z] ] := 
    (2 b - 2) (z-1)^(b-3/2) z^(1/2-b) (
       -(-z)^(1/2) (1-z)^(-1/2) Log[Sqrt[-z]+Sqrt[1-z]] +
       Sum[ z^k (z-1)^(-k) (2 k-1)^(-1),{k,1,b-3/2}])/;
 IntegerQ[b-1/2] && b>1 

 F21[ 1,1,1/2,z_/;Abs[z]=!=1 && Not[Znak[z]] ] := 
    (1-z)^(-1) (1 + Sqrt[z] ArcSin[Sqrt[z]]/Sqrt[1-z])

(************************* Formula 1, p. 486, PBM **************************)

 F21[ a_,b_,1/2,z_/;Znak[z] ] := (1-z)^(-a) Cos[2 a ArcTan[Sqrt[-z]]]/;
 Expand[b-a-1/2]===0

 F21[ b_,a_,1/2,z_/;Znak[z] ] := (1-z)^(-a) Cos[2 a ArcTan[Sqrt[-z]]]/;
 Expand[b-a-1/2]===0

(************************* Formula 2, p. 486, PBM **************************)

 F21[ a_,b_,3/2,z_/;Znak[z] ] := 
   (1-z)^(1/2-a) Sin[Expand[(2 a-1) ArcTan[Sqrt[-z]]]]/((2 a-1) Sqrt[-z])/;
 Expand[b-a-1/2]===0 && a=!=1/2

 F21[ b_,a_,3/2,z_/;Znak[z] ] := 
   (1-z)^(1/2-a) Sin[Expand[(2 a-1) ArcTan[Sqrt[-z]]]]/((2 a-1) Sqrt[-z])/;
 Expand[b-a-1/2]===0 && a=!=1/2

(************************* Formula 9, p. 487, PBM **************************)

 F21[ 1,a_Rational?Positive,b_,z_/;Abs[z]=!=1 ] := 
    a Sum[ (-1)^(k-1) (-z)^(-k)/(a-k),{k,1,a-1/2}] +
    a (-z)^(-a) (2 Sin[a Pi] ArcTan[(-z)^(1/2)]-Cos[a Pi] Log[1+z])/;
 Denominator[a]==2 && b===a+1 && Znak[z]

 F21[ a_Rational?Positive,1,b_,z_/;Abs[z]=!=1 ] := 
    a Sum[ (-1)^(k-1) (-z)^(-k)/(a-k),{k,1,a-1/2}] +
    a (-z)^(-a) (2 Sin[a Pi] ArcTan[(-z)^(1/2)]-Cos[a Pi] Log[1+z])/;
 Denominator[a]==2 && b===a+1 && Znak[z]

(************************* Formula 131, p. 463, PBM ************************)
(*                              &&                                         *)
(************************* Formula 8, p. 486, PBM **************************)

 F21[ 1,a_Rational/;0<a<1,b_Rational,z_/;Abs[z]=!=1&&Not[Znak[z]] ] := 
    F21Formula131P463[Numerator[a],Denominator[a],z]/;b-a-1==0

 F21[ 1,a_Rational/;0<a<1,b_Rational,z_/;Abs[z]=!=1&&Znak[z] ] := 
    F21Formula8P486[Numerator[a],Denominator[a],-z]/;b-a-1==0

 F21[ 1,a_Rational/;a>1,b_Rational,z_/;Abs[z]=!=1 ] := 
    a/(z (a-1))(-1 + F21[1,a-1,a,z])/;b-a-1==0

 F21Formula131P463[ m_,n_,z_ ] := 
    Module[ {arg = z^(1/n), par1 = 2 Pi/n, par2 = 2 Pi m/n},
      -m z^(-m/n)/n ( Log[1-arg] + (-1)^m/2 (1+(-1)^n) Log[1+arg] +
      Sum[ Cos[par2 k] Log[1-2arg Cos[par1 k] + arg^2] -
           2Sin[par2 k] ArcTan[arg Sin[par1 k]/(1-arg Cos[par1 k])],
      {k,1,Floor[(n-1)/2]}])
    ]

 F21Formula8P486[ m_,n_,z_ ] := 
    Module[ {arg = z^(1/n), par1 = Pi/n, par2 = Pi m/n},
      -m z^(-m/n)/n ( (-1)^m (1-(-1)^n) Log[1+arg]/2 +
      Sum[ Cos[par2 (2k+1)] Log[1-2arg Cos[par1 (2k+1)] + arg^2] -
           2Sin[par2 (2k+1)] ArcTan[arg Sin[par1 (2k+1)]/
           (1-arg Cos[par1 (2k+1)])],
      {k,0,Floor[n/2]-1}]) 
    ]

(************************* Formula 5, p. 491, PBM ****************************)

 F21[ a_,b_,c_,1/2 ] := Sqrt[Pi] Gamma[c]/(Gamma[a/2+1/2] Gamma[b/2+1/2])/;
 Expand[2 c-a-b-1]===0

(************************* Polynomials ***************************************)

(************************* Formula 203, p. 467, PBM **************************)

 F21[ n_,m_,1/2,z_/;z=!=1 ] := ChebyshevT[-2n,Sqrt[1-z]]/;
 ZnakSum[n] && Expand[m+n]===0

(************************* Formula 207, p. 468, PBM **************************)

 F21[ a_,b_/;ZnakSum[b],c_,z_/;Not[Znak[z]] ] := 
    2^c/z^b Cos[2 b ArcCos[1/Sqrt[z]]]/;
 NotFreeSymbolicSum[] &&
 Expand[a-b-1/2]===0 && Expand[c-2 a]===0

 F21[ b_/;ZnakSum[b],a_,,c_,z_/;Not[Znak[z]] ] := 
    2^c/z^b Cos[2 b ArcCos[1/Sqrt[z]]]/;
 NotFreeSymbolicSum[] &&
 Expand[a-b-1/2]===0 && Expand[c-2 a]===0

(************************* Formula 203, p. 467, PBM **************************)

 F21[ n_,m_,c_,z_/;Not[Znak[z]] ] := (-1)^(-4n) (-2n)! *
    If[ Expand[c-2 n]===0, 1/(-2 n)!, Gamma[c]/Gamma[c-2n] ] *
    (z/4)^(-n) GegenbauerC[-2n,1-c+2n,1/Sqrt[z]]/;
 NotFreeSymbolicSum[] &&
 ZnakSum[n] && Expand[m-n-1/2]===0 && (!NumberQ[n] || !NumberQ[m])

 F21[ m_,n_,c_,z_/;Not[Znak[z]] ] := (-1)^(-4n) (-2n)! *
    If[ Expand[c-2 n]===0, 1/(-2 n)!, Gamma[c]/Gamma[c-2n] ] *
    (z/4)^(-n) GegenbauerC[-2n,1-c+2n,1/Sqrt[z]]/;
 NotFreeSymbolicSum[] &&
 ZnakSum[n] && Expand[m-n-1/2]===0 && (!NumberQ[n] || !NumberQ[m])

 NotFreeSymbolicSum[___] := If[ Names["SymbolicSum"] =!= {}, True, False ]

(************************** Elliptic *****************************************)

 F21[ 1/4,3/4,1,z_/;Not[Znak[z]] ] := 
        2 PowerExpand[EllipticK[2 Sqrt[z]/(1+Sqrt[z])]/
           (Pi Sqrt[1+Sqrt[z]])]

 F21[ 3/4,1/4,1,z_/;Not[Znak[z]] ] := 
        2 PowerExpand[EllipticK[2 Sqrt[z]/(1+Sqrt[z])]/
           (Pi Sqrt[1+Sqrt[z]])]

 F21[ 1/2,1/2,1,z_ ]  := 2 EllipticK[z]/Pi

 F21[ -1/2,1/2,1,z_ ] := 2 EllipticE[z]/Pi

 F21[ 1/2,-1/2,1,z_ ] := 2 EllipticE[z]/Pi

 F21[ 1/2,3/2,1,z_ ] := 2 EllipticE[z]/(Pi (1-z))

 F21[ 3/2,1/2,1,z_ ] := 2 EllipticE[z]/(Pi (1-z))

 F21[ __ ] := HyperFail


(**************************************************************************** 
*                         Hypergeometric1F1 
*
*****************************************************************************)

 GeneralHypergeometricPFQ[ {uppar_},{lowpar_},arg_ ] := 
   Hypergeometric1F1[ uppar,lowpar,arg] /;
 Accuracy[{uppar,lowpar,arg}] < Infinity

 GeneralHypergeometricPFQ[ {uppar_},{lowpar_},arg_ ] := 
  Module[ {r},
    r = F11[ uppar,lowpar,arg ];
    r = If[ FreeQ[r,HyperFail],r,
            E^arg F11[lowpar-uppar,lowpar,-arg] ];
    If[ FreeQ[r,HyperFail],r,KummerFlag[]:=False;
        r = Hypergeometric1F1[uppar,lowpar,arg]];
        KummerFlag[] := True;
        r
  ]
 
 KummerFlag[] := True

 F11[ 1,3/2,z_ ] := Sqrt[Pi] Erf[Sqrt[z]] E^z /(2 Sqrt[z])

 F11[ 1/2,3/2,z_ ] := Sqrt[Pi] Erf[Sqrt[-z]] /(2 Sqrt[-z])

(************************* Formula 11 , p. 579, PBM *************************)

 F11[ n_Integer?Positive,b_,z_/;Not[Znak[z]] ] := (b-1)/(n-1)!*
   Module[{var},D[var^(n-b) E^var Gamma[b-1,0,var],{var,n-1}]/.var->z]/;
 Not[NumberQ[b]] || Re[b]>1

(************************* Formula 5 , p. 579, PBM **************************)

 F11[ a_,b_,z_ ] := Gamma[a+1/2] z^(1/2-a) 2^(b-1) E^(z/2)*
    BesselI[a-1/2,z/2]/;
 Expand[b-2 a]===0

(************************* Formula 1 , p. 583, PBM **************************)

 F11[ a_,b_,z_/;Znak[z] ] := a (-z)^(-a) Gamma[a,0,-z]/;b-a-1===0

(************************* Formula 17 , p. 579, PBM *************************)

 F11[ a_/;Znak[a],b_,z_ ] := (-a)! LaguerreL[-a,b-1,z] Gamma[b]/Gamma[b-a]

 F11[ __ ] := HyperFail


(**************************************************************************** 
*                         Hypergeometric1F2 
*
*****************************************************************************)

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
   Module[ { answer = F12[ uppar[[1]],lowpar,arg ] },
     answer/;FreeQ[answer,HyperFail]
   ]/;
  Length[uppar] == 1 && Length[lowpar] == 2

 F12[ 1/2,{3/2,3/2},arg_/;Znak[arg] ] := 
   SinIntegral[Sqrt[-arg] 2]/(Sqrt[-arg] 2) 
 
 F12[ __ ] := HyperFail

(**************************************************************************** 
*                         Hypergeometric2F2 
*
*****************************************************************************)

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
   Module[ { answer = F22[ uppar,lowpar,arg ] },
     answer/;FreeQ[answer,HyperFail]
   ]/;
  Length[uppar] == 2 && Length[lowpar] == 2

 F22[ {1,1},{2,2},z_/;Not[Znak[z]]] := (ExpIntegralEi[z] - Log[z] -
	EulerGamma)/z

 F22[ {1,1},{2,2},z_/;Znak[z]] := (ExpIntegralEi[z] - Log[-z] -
	EulerGamma)/z

 F22[ __ ] := HyperFail

(**************************************************************************** 
*                         Hypergeometric2F0 
*
*****************************************************************************)

 GeneralHypergeometricPFQ[ uppar_,{},arg_ ] := 
   Module[ { answer = F20[ uppar,arg ] },
     answer/;FreeQ[answer,HyperFail]
   ]/;
  Length[uppar] == 2 

 F20[ {n_,m_},z_/;Znak[z] ] := 2^(2 n) (-z)^(-n) HermiteH[-2 n,1/Sqrt[-z]]/;
 ZnakSum[n] && Expand[m-n-1/2]===0

 F20[ {m_,n_},z_/;Znak[z] ] := 2^(2 n) (-z)^(-n) HermiteH[-2 n,1/Sqrt[-z]]/;
 ZnakSum[n] && Expand[m-n-1/2]===0

 F20[ __ ] := HyperFail

(**************************************************************************** 
*                         Hypergeometric2F3 
*
*****************************************************************************)

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
   Module[ {
     answer = F23[ uppar,lowpar,arg ] },
     answer/;FreeQ[answer,HyperFail]
   ]/;
  Length[uppar] == 2 && Length[lowpar] == 3

 F23[ uppar_,{v1___,1/2,v2___},arg_ ] := 
   F23ToBessel1[uppar,Sort[{v1,v2}],arg]

 F23[ uppar_,{v1___,3/2,v2___},arg_ ] := 
   F23ToBessel2[uppar,Sort[{v1,v2}],arg]

 F23[ __ ] := HyperFail

(************************* Formula 6 - 7 , p. 609, PBM **********************)

 F23ToBessel1[ {a_,a1_},{b_,b1_},arg_/;Not[Znak[arg]] ] :=
   Cosh[Sqrt[arg]] BesselI[b1-1,Sqrt[arg]] Gamma[b1] arg^(1/2-b1/2)/;
 Expand[a1-a-1/2]===0 && Expand[b1-b-1/2]===0 && Expand[2 a-b]===0

 F23ToBessel1[ {a_,a1_},{b_,b1_},arg_/;Znak[arg] ] :=
   Cos[Sqrt[-arg]] BesselJ[b1-1,Sqrt[-arg]] Gamma[b1] (-arg)^(1/2-b1/2)/;
 Expand[a1-a-1/2]===0 && Expand[b1-b-1/2]===0 && Expand[2 a-b]===0

 F23ToBessel1[ __ ] := HyperFail

(************************* Formula 8 - 9 , p. 609, PBM **********************)

 F23ToBessel2[ {1,1},{2,2},arg_/;Znak[arg] ] :=
   -EulerGamma/arg - Log[2 Sqrt[-arg]]/arg + CosIntegral[2 Sqrt[-arg]]/arg

 F23ToBessel2[ {a_,a1_},{b_,b1_},arg_/;Not[Znak[arg]] ] :=
   Sinh[Sqrt[arg]] BesselI[b-1,Sqrt[arg]] Gamma[b]/arg^(b/2)/;
 Expand[a1-a-1/2]===0 && Expand[b1-b-1/2]===0 && Expand[2 a-b1]===0

 F23ToBessel2[ {a_,a1_},{b_,b1_},arg_/;Znak[arg] ] :=
   Sin[Sqrt[-arg]] BesselJ[b-1,Sqrt[-arg]] Gamma[b]/(-arg)^(b/2)/;
 Expand[a1-a-1/2]===0 && Expand[b1-b-1/2]===0 && Expand[2 a-b1]===0

 F23ToBessel2[ __ ] := HyperFail

(*========================================================================*)

 GeneralHypergeometricPFQ[ uppar_,lowpar_,arg_ ] := 
     HypergeometricPFQ[ Sort[uppar],Sort[lowpar],arg ]
 
 KeiperFlag = True

 Hypergeometric2F1[ a_,b_,c_,arg_ ] :=
   Module[{r},
     KeiperFlag = False;
     r = Hypergeometric2F1[a,b,c,arg];
     KeiperFlag = True;
     r /; Head[r]=!=Hypergeometric2F1 
   ] /; KeiperFlag

 Hypergeometric2F1[ a_,b_,c_,arg_ ] :=
   Block[ { r, answer, Hypergeometric2F1 },
        answer = HypergeometricPFQ[{a,b},{c},arg];
        If[ Length[r=Expand[answer]]<10,
            Simplify[r],
            answer ] /; Head[answer]=!=Hypergeometric2F1
   ] /; Accuracy[{a,b,c,arg}] === Infinity && GaussFlag &&
 KeiperFlag

 Hypergeometric1F1[ a_,b_,arg_ ] :=
   Block[ { answer,Hypergeometric1F1 },
            answer = HypergeometricPFQ[{a},{b},arg];
     answer/;Head[answer]=!=Hypergeometric1F1
   ] /; Accuracy[{a,b,arg}] === Infinity && KummerFlag[]

(***************************************************************************
*                     HypergeometricPFQRegularized
*
****************************************************************************)
HypergeometricPFQRegularized[{ }, { }, z_] := Exp[z];
HypergeometricPFQRegularized[{a_}, { }, z_] := (1 - z)^(-a);
HypergeometricPFQRegularized[{ }, {c_}, z_] :=
        Hypergeometric0F1Regularized[c, z];
HypergeometricPFQRegularized[{a_}, {c_}, z_] :=
        Hypergeometric1F1Regularized[a, c, z];
HypergeometricPFQRegularized[{a_, b_}, {c_}, z_] :=
        Hypergeometric2F1Regularized[a, b, c, z];

(* Sorting *)
HypergeometricPFQRegularized[a_List, b_List, z_] :=
        HypergeometricPFQRegularized[Sort[a], Sort[b], z] /;
                (!OrderedQ[a] || !OrderedQ[b])
 
(* Cancelation *)
HypergeometricPFQRegularized[{a0___, a1_, a2___}, 
                             {b0___, b1_, b2___}, z_] :=
    HypergeometricPFQRegularized[{a0, a2}, {b0, b2}, z]/Gamma[b1] /; 
a1 == b1
 
(* Negative Integers in the Denominator *)
HypergeometricPFQRegularized[a_List, {b0___, c_, b2___}, z_] :=
    Module[{b = {b0, b2}, c1 = 1-c},
        Apply[Times, Map[Pochhammer[#, c1]&, a]] z^c1 *
        HypergeometricPFQRegularized[a+c1, Append[b, 1]+c1, z] /
        Apply[Times, Map[Pochhammer[#, c1]&, b]]
    ] /; (c == Round[c] && Round[c] <= 0)
 
(* Approximate Numerical Evaluation *)
HypergeometricPFQRegularized[a_List, b_List, z_] :=
    Module[{pfq = HypergeometricPFQ[a, b, z]},
        pfq/N[Apply[Times, Map[Gamma, b]], Precision[pfq]] /; NumberQ[pfq]]
 
(* Derivatives wrt z *)
Derivative[0, 0, n_Integer][HypergeometricPFQRegularized] ^:=
    (Apply[Times, Map[Pochhammer[#, n]&, #1]] *
        HypergeometricPFQRegularized[#1 + n, #2 + n, #3] /
        Apply[Times, Map[Pochhammer[#, n]&, #2]])& /; n > 0

(*========================================================================*)

 Znak[n_ a_] := True/;NumberQ[n] && Im[n]===0 && n<0
 
 Znak[Complex[0,n_] a_.] := True/;n<0
 
 Znak[n_] := True/;NumberQ[n] && Im[n]===0 &&n<0

 Znak[a_] := False

 ZnakSum[n_ a_] := True/;NumberQ[n] && Im[n]===0 && n<0
 
 ZnakSum[Complex[0,n_] a_.] := True/;n<0
 
 ZnakSum[n_] := True/;NumberQ[n] && Im[n]===0 &&n<0
 
 ZnakSum[n_ + c_] := True/;ZnakSum[n] && ZnakSum[c]
 
 ZnakSum[n_ + c_] := True/;NumberQ[n] && Znak[c]

 ZnakSum[a_] := False

 LogTrig = {
  Cos[ArcSin[w_]]^m_. :> (1-w^2)^(m/2),
  Cos[ArcTan[w_]]^m_. :> 1/(1+w^2)^(m/2),
  Sec[ArcCos[w_]]^m_. :> w^(-m),
  Sec[ArcSin[w_]]^m_. :> (1-w^2)^(-m/2),
  Sec[ArcTan[w_]]^m_. :> (1+w^2)^(m/2),
  Sin[ArcCos[w_]]^m_. :> (1-w^2)^(m/2),
  Sin[ArcTan[w_]]^m_. :> w^m/(1+w^2)^(m/2),
  Csc[ArcSin[w_]]^m_. :> w^(-m),
  Csc[ArcCos[w_]]^m_. :> (1-w^2)^(-m/2),
  Csc[ArcTan[w_]]^m_. :> (1+w^2)^(m/2)/w^m,
  Tan[ArcSin[w_]]^m_. :> w^m/(1-w^2)^(m/2),
  Tan[ArcCos[w_]]^m_. :> (1-w^2)^(m/2)/w^m,
  Cot[ArcCos[w_]]^m_. :> w^m/(1-w^2)^(m/2),
  Cot[ArcSin[w_]]^m_. :> (1-w^2)^(m/2)/w^m,
  ArcTan[Tan[w_]]^m_. :> w^m
        }

 MultPochham[a_,k_] := Times@@(Pochhammer[ #,k]&/@a)

 SimpPolyGammaH[ expr_Plus ] := 
    Module[ {inter},
      inter = Plus@@Cases[expr,a_/;FreeQ[a,PolyGamma]];
      If[ Length[inter] < 7,
          Factor[inter],
          inter ]
      + SimpPolyGammaH1[ expr - inter ]
    ]

 SimpPolyGammaH1[ expr_Plus ] :=
    Module[ {listP,argP,exprN},
      listP = Union[ List@@(#&/@expr) /.
              { a_. PolyGamma[n_,w_] :> {n,w} }];
      argP  = FindArgPolyGamma[ listP ]; 
      If[ Length[argP] != 0,
          SimpPolyGammaH[ Expand[expr//.
             BuildRule[(PolyGamma@@#)&/@argP,
                       ((PolyGamma@@{#[[1]],#[[2]]-1})&/@argP) + 
                       ((Together[#[[2]]-1]^(-1-#[[1]]) (#[[1]])!*
                        (-1)^(#[[1]]) )&/@argP)]]],
          exprN = 
             CollectPolyGamma[ expr,listP ]/.
             {a_ PolyGamma[w__] :> Factor[a] PolyGamma[w]};
          If[ Length[exprN] < 6 && 
              Length[FindArg2PolyGamma[ listP ]] != 0,
              exprN//.SimpPolyGammaSum,
              exprN ]
        ]
    ]

 SimpPolyGammaH1[ expr_ ] := expr 

 SimpPolyGammaH[ expr_ ] := expr 

 BuildRule[{l_,lR___},{r_,rR___}] :=
   Join[{l :> r},BuildRule[{lR},{rR}]]

 BuildRule[{},{}] := {}
        
 FindArgPolyGamma[ {v1___,{nw_,w_},v2___,{nv_,v_},v3___} ] :=
   Join[{{nv,v}},FindArgPolyGamma[ {v1,v2,v3} ]] /; 
 IntegerQ[v-w] && Positive[v-w] && nw==nv

 FindArgPolyGamma[ {___} ] := {}

 FindArg2PolyGamma[ {v1___,{nw_,w_},v2___,{nv_,v_},v3___} ] :=
   Join[{{nv,v}},FindArg2PolyGamma[ {v1,v2,v3} ]] /; 
 nw==nv && (v+w===0 || v+w===1)  

 FindArg2PolyGamma[ {___} ] := {}

 CollectPolyGamma[ expr_,{ {n_,arg_},w___ } ] :=
   CollectPolyGamma[ Collect[expr,PolyGamma@@{n,arg} ], {w} ]

 CollectPolyGamma[ expr_,{} ] := expr

 SimpGammaH[f_Plus] := Map[ SimpGammaH[#]&,f ]

 SimpGammaH[ expr_/;Length[expr]>1 ] :=
     (expr/.{Gamma[_]:>1}) *
     SimpGamma1[ Times@@Cases[expr,Gamma[a_]^n_.]] 

 SimpGamma1[ Times[v1___,Gamma[w1_]^n_.,v2___,Gamma[w2_]^m_.] ] := 
     If[ (w2-w1)>0,
           SimpGamma1[v1 v2 ]/Factor[Pochhammer[w1,w2-w1]^n],
           SimpGamma1[v1 v2 ] Factor[Pochhammer[w2,w1-w2]^n] 
       ] /; IntegerQ[w2-w1] && IntegerQ[n] && n>0 && n+m == 0

 SimpGamma1[Times[v1___,Gamma[u_]^n_.,v2___,Gamma[v_]^m_.]] := 
     If[ SameQ[Expand[v-u],1/2],(2^Expand[1-2 u] Pi^(1/2))^Sign[m] *
         SimpGamma1[ Gamma[u]^(n-Sign[m])*
                     Gamma[v]^(m-Sign[m]) v1 v2 *
                     Gamma[Expand[2 u]]^Sign[m] ],
         (2^Expand[1-2 v] Pi^(1/2))^Sign[n] *
         SimpGamma1[ Gamma[u]^(n-Sign[n])*
                     Gamma[v]^(m-Sign[n]) v1 v2 *
                     Gamma[Expand[2 v]]^Sign[n] ]
      ]/;
  Abs[Expand[u-v]]===1/2 

 SimpGamma1[Times[v1___,Gamma[u_]^n_.,v2___,Gamma[v_]^m_.]] := 
     SimpGamma1[ v1 v2 Gamma[u]^(n-Sign[n]) Gamma[v]^(m-Sign[m]) ] *
     (Pi/Sin[Expand[Pi u]])^Sign[n] /;
   Expand[u+v]===1 && m n > 0

 SimpGamma1[Times[v1___,Gamma[1+u_]^n_.,v2___,Gamma[1+v_]^m_.]] := 
    u^Sign[n] (Pi/Sin[Expand[Pi u]])^Sign[n]*
    SimpGamma1[ v1 v2 Gamma[1+u]^(n-Sign[n]) Gamma[1+v]^(m-Sign[m]) ] /;
   Expand[u+v]===0 && m n > 0

 SimpGamma1[v__] := v

 SimpGammaH[v__] := v

 SimpPolyGammaSum = {
  a_. PolyGamma[k_,z_] + b_ PolyGamma[k_,x_] +c_.:>
    c + a Pi Module[{var},D[Cot[Pi var],{var,k}]/.var->x]/;EvenQ[k] && 
    Znak[b] && Expand[a+b]===0 && Expand[z+x-1]===0,
  a_. PolyGamma[k_,z_] + b_ PolyGamma[k_,x_] +c_.:>
    c - a Pi Module[{var},D[Cot[Pi var],{var,k}]/.var->x]/;
    OddQ[k] && Expand[a-b]===0 && Expand[z+x-1]===0,
  a_. PolyGamma[k_,z_] + b_ PolyGamma[k_,x_] +c_.:>
    c + a k! x^(-k-1) + a Pi Module[{var},
    D[Cot[Pi var],{var,k}]/.var->x]/;
    EvenQ[k] && Znak[b] && Expand[a+b]===0 && Expand[z+x]===0,
  a_. PolyGamma[k_,z_] + b_ PolyGamma[k_,x_] +c_.:>
    c + a k! x^(-k-1) - a Pi Module[{var},
    D[Cot[Pi var],{var,k}]/.var->x]/;
    OddQ[k] && Expand[a-b]===0 && Expand[z+x]===0,
  Cot[a_] :> Cot[Expand[a]],
  Csc[a_] :> Csc[Expand[a]]
      }

 SimpTrigSum = {
  a_. Tan[x_] + b_. Cot[x_] +c_.:> c + 2 a Csc[2 x]/;a===b,
  a_. Tanh[x_] + n_?Negative b_. Coth[x_] +c_.:> c-2 a Csch[2 x]/;
      a+b n===0&&Not[Znak[a]],
  n_?Negative a_. Tanh[x_] + b_. Coth[x_] +c_.:> c+2 b Csch[2 x]/;
      a n+b===0&&Not[Znak[b]]
  }

 PolyGammaRule = {
  PolyGamma[v_]     :> PolyGamma[0,v],
  PolyGamma[k_Integer,n_Integer + v_] :> 
    PolyGamma[k,n+v-1] + (-1)^k k!/(n+v-1)^(k+1)/;n>0,
  PolyGamma[k_Integer,n_] :> PolyGamma[k,n-1] +(-1)^k k! (n-1)^(-k-1)/;
    NumberQ[n] && Re[n]>1,
  PolyGamma[n_,v_]  :> PolyGamma[n,Expand[v]]
                 }
(*========================================================================*)

 Unprotect[ PolyLog ]

 PolyLog[2,1/2] := Pi^2/12 - Log[2]^2/2

 PolyLog[2,I] := -Pi^2/48 + I Catalan

 PolyLog[2,-I] := -Pi^2/48 - I Catalan

(*========================================================================*)

 Unprotect[ Zeta ]

 Zeta[z_,1/2] := (2^z-1) Zeta[z]

 Zeta[2,1/4] := Pi^2 + 8 Catalan

 Zeta[2,3/4] := Pi^2 - 8 Catalan

 Zeta[n_Integer,v_Rational] := Zeta[n,v-1] - (v-1)^(-n)/;v>1

(*========================================================================*)

 Unprotect[ PolyGamma ]

 PolyGamma[1,1/4] := Pi^2 + 8 Catalan

 PolyGamma[1,3/4] := Pi^2 - 8 Catalan
 
 PolyGamma[0,n_Rational] := -EulerGamma - Log[2 Denominator[n]] -
    Pi Cos[Pi n]/(2 Sin[Pi n]) + 2 Sum[Cos[2 n i Pi] *
            Log[Sin[Pi i/Denominator[n]]],
            {i,1,Floor[(Denominator[n]-1)/2]}] /; 
 n>0 && n<1

(*========================================================================*)

(* Protect[ Hypergeometric2F1, Hypergeometric1F1, Zeta, PolyLog, PolyGamma ] *)

 End[ ] (* "HypergeometricPFQ`Private` *)

 SetAttributes[HypergeometricPFQ, { ReadProtected, Protected } ];
 SetAttributes[HypergeometricPFQRegularized, ReadProtected];

 End[ ] (* "System`" *)

