
(* :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.
*)

(* :Version: Mathematica 2.1 *)

(* :Title: Gamma-functions Simplifications *)

(* :Author:
	Victor S. Adamchik, October 1991.
*)

(* :Context: Simplify`Gamma` *)

(* :Keywords: gamma *)

(* :Requirements: none. *)

(* :Warnings: none. *)

(* :Limitations: none known. *)

(* :Summary:
	This package allows to simplify any expressions with 
	Gamma-functions using a sequence of 8 rules which
	apply to the expression. 
*)

(*========================================================================*)
 
 Begin["System`"]

 SimplifyGamma::usage =
"SimplifyGamma[expr] performs a sequence of gamma-function transformations 
   on expr, and returns the simplest form it finds."

 Begin["Simplify`Gamma`Private`"]


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

 Unprotect[ SimplifyGamma]

 SimplifyGamma[ expr_ ] := SimpGamma[Expand[expr]]

 SimpGamma[ expr_Plus/;Not[FreeQ[expr,Gamma[u__]/;Length[{u}]>1]] ] :=
    Module[ {answer,r},
      answer = SimpIncompleteGamma[expr,
			Union[(r = Cases[ expr,a_. Gamma[u_,n_,m_] ])/.
				b_. Gamma[v_,n_,m_] :> v] ];
      answer = SimpIncompleteGamma1[Expand[answer],
            Union[r/.b_. Gamma[v_,n_,m_] :> Gamma[v,n,m]] ];
      answer /; FreeQ[answer,FailInt]
    ]

 SimpGamma[f_Plus] := (SimpGamma[#]&/@f)/;!FreeQ[f,Gamma[_]]

 SimpGamma[ expr_ ] := 
   SimpGamma[Numerator[expr]] / SimpGamma[Denominator[expr]] /;
 Head[Denominator[expr]]===Plus && !FreeQ[Denominator[expr],Gamma[_]]  

 SimpGamma[ expr_ ] := If[ FreeQ[expr,Gamma[_]], expr, SimpGammaG[expr] ]

 SimpIncompleteGamma[ expr_,args_/;Length[args] >1 ] :=
    Module[ {el,min,var},
      el = Union[ args/.n_. + v_. :> v/;NumberQ[n] ];
      If[ Length[el] >1, FailInt, 
          argn = args/.el[[1]]->var;
          min =  Min[argn/.n_. + v_. var :> n/;NumberQ[n]];
          SimpGamma[ (expr/.el[[1]]->var)//.Gamma[k_. + v_. var,n_,m_] :>
          (k+v var-1) Gamma[k+v var-1,n,m] - m^(k+v var-1) E^(-m)/;
          NumberQ[k] && k>min ]/.var->el[[1]] ]
     ]

 SimpIncompleteGamma[ expr_,args_ ] := FailInt

 SimpIncompleteGamma1[ expr_,{arg1_,arg___} ] :=
    Module[ {r,rnew},
      rnew = 
      Factor[Plus@@((r=Cases[expr,a_. arg1])/.a_. arg1 :> a)] arg1; 
      SimpIncompleteGamma1[ (expr - Plus@@r)+rnew,{arg} ]
    ]

 SimpIncompleteGamma1[ expr_, _ ] := expr

 SimpGammaG[ Gamma[a_]^n_.] := Gamma[a]^n

 SimpGammaG[ expr_Times  ] := 
   Module[ { rr = Cases[expr,Gamma[a_]^n_./;!IntegerQ[n]],
		exprnew, rrr },
     exprnew = expr/(rrr=Times@@rr *Times@@Cases[expr,a_/;FreeQ[a,Gamma[_]]]);
     rrr *
     If[ FreeQ[ exprnew,Gamma[_]] || MatchQ[exprnew,Gamma[a_]^n_.],
	 exprnew,
	 rrr = Cases[exprnew,Gamma[a_]^n_.];
	 rr = exprnew/Times@@rrr;
	 If[ rr === exprnew, SimpGamma[#]&/@ exprnew,
	     SimpGamma[rr] SimpGammaN[ rrr/.Gamma[a_]^n_.:>{Expand[a],n} ]
	 ]
    ]
  ]

 SimpGammaG[ f_[expr___] ] := f@@(SimpGamma[#]&/@{expr})

 SimpGammaN[ expr_ ] := 
    If[ Length[expr]<2, 
	SimpGamma1[ expr ],
	Block[ { rr =  Cases[expr,{n_Rational+s_.,a_}], SimpGamma1 },
	  rr = 
	    SimpGamma1[Complement[expr,rr]] * 
	    MultiplicationVar[Select[rr,#[[2]]>0&]] *
	    MultiplicationVar[#/.{a_,b_}:>{a,-b}&/@Select[
			rr,#[[2]]<0&]]^(-1)//.{
		gamma[a_]^n_ :> gamma[#/.{c_,d_}:>{c,d n}&/@a],
		gamma[a_] gamma[b_] :> gamma[Join[a,b]] };
	  rr/.gamma->SimpGamma1/.
		SimpGamma1[a_] SimpGamma1[b_] :> SimpGamma1[Join[a,b]]
	]
    ] 
 
 SimpGamma1[ {} ] := 1

 SimpGamma1[ {v1___, {1,n_}, v2___} ] := SimpGamma1[{v1,v2}]

 SimpGamma1[ {v1___, {u_,n_/;Abs[n]>1}, v2___} ] := 
	SimpGamma1[ Join[{v1,v2},Table[{u,Sign[n]}, {i,Abs[n]}]] ]

 SimpGamma1[ {v1___, {u1_,n1_}, v2___, {u2_,n2_},v3___ } ] := 
	SimpGamma1[ {v1,v2,v3} ]/;n1+n2==0 && u1===u2

 SimpGamma1[ expr_ ] :=
   If[ Length[expr]==1, Gamma[expr[[1,1]]]^expr[[1,2]] ,
   Module[ { new, exprnew=expr, i=1, answer=1 },
     While[ i<=7 && Length[exprnew]!=0,
	answer *=
	Times@@(Fold[ function[#1[[1]], #1[[2]], #2[[1]], #2[[2]], i]&,
		      #[[1]], Rest[#] ]&/@(new=FindListCond[exprnew, i]));
	If[ !FreeQ[answer,gamma],True,i+=1];
	exprnew = Join[ ComplementNotSet[exprnew,Flatten[new,1]],
			Flatten[Cases[answer,gamma[a_]^n_.]/.
			  gamma[a_]^n_.:>Table[{a,Sign[n]},{i,Abs[n]}],1]
		  ];
	answer = answer/.gamma[_]:>1
     ];
     answer Times@@(Power[Gamma[ #[[1]] ],#[[2]]]&/@exprnew) 
   ]
  ]

 MultiplicationVar[ expr_ ]:=
    If[ Length[expr]<2, 
	gamma[expr]
	,
	Module[ { r, rr },
	  r = (rr=Cases[expr, {n_Rational + s_,_}])/.
			{n_Rational + s_,_}:> s;
	   Times@@(Multiplication[Sort[#]]&/@(
	   Cases[rr, {n_Rational + #,_}]&/@Union[r])) *
	   Multiplication[ComplementNotSet[expr,rr]]
	]
   ]

 Multiplication[ {} ] := 1
 
 Multiplication[ expr_ ] :=
  If[ Length[expr]<2,
      gamma[expr],
      Module[ {rr = expr/.{a_,b_/;b>1} :> {a,1}, 
	       rest, univ, max, flag, var},
	rest = Join[ SameElem[rr], 
		#/.{a_,b_} :>{a,b-1}&/@Select[expr, #[[2]]>1&]];
	rr = Union[rr];
	var = expr[[1,1]]/.{n_Rational+s_. :> s+Floor[n]};
	flag = If[Positive[expr[[1,1]]-var],1,-1]; 
	max = Max[(#/.{s_,_}:>Denominator[s-var])&/@expr];
	univ = Complement[Table[ {i/max flag+var,1}, {i,max-1}], rr];
	If[ max - Length[univ] > Floor[max/2] + 1,
	    MultiplicationResult[var,max,rr,flag,rest,univ]
	    ,
	    If[ GCD@@((#/.{s_,_}:>Denominator[s-var])&/@expr) != 1, 
		var = expr[[1,1]]; 
	        max = Max[(#/.{s_,_}:>Denominator[s-var])&/@expr];
		univ = Complement[Table[ {i/max flag+var,1}, {i,max-1}], rr];
		If[ max!=1 && max - Length[univ] > Floor[max/2] + 1,
		    MultiplicationResult[var,max,rr,flag,rest,univ]
	            ,
		    rr = Cases[expr,{a_/;Denominator[a-var]==max,_}];
		    If[rr=!=expr,Multiplication[rr],gamma[rr] ] *
		    Multiplication[ 
			Delete[expr,Flatten[Position[expr,#]&/@rr,1]] ]
		]
	        ,
		rr = Cases[expr,{a_/;Denominator[a-var]==max,_}];
		gamma[rr] *
		Multiplication[ Delete[expr,Flatten[Position[expr,#]&/@rr,1]]]
	   ]
	]
     ]
  ]

 MultiplicationResult[ var_,max_,rr_,flag_,rest_,univ_ ] :=
	If[ var===0 || (IntegerQ[var] && Negative[var]),
	    max^(-var max+1/2) (2 Pi)^((max-1)/2) *
	    (-1)^(-var-var max) (-var)!/(max (-max var)!)
	    , 
	    max^Expand[-var max+1/2] (2 Pi)^((max-1)/2) *
	    gamma[{{Expand[max var],1},{var,-1}}]
	] *
	Multiplication[Join[
	    Complement[Join[rr,univ],Table[ {i/max flag+var,1},{i,max-1}]],
		     rest]] *
	gamma[univ/.{a_,1} :> {a,-1}]

 FindListCond[ {v1___, {u_,n_}, v2___, {v_,m_}, v3___}, c_ ] := 
   Join[ {{{u,n},{v,m}}}, FindListCond[{v1,v2,v3},c] ]/; cond[u,n,v,m,c]

 FindListCond[__] := {}

 cond[u_, n_, v_, m_, 1] :=  m n > 0 && IntegerQ[u+v] && !IntegerQ[u]
 function[u_, n_, v_, m_, 1] := 
    If[ Positive[u+v],
	Pochhammer[-v,u+v]^Sign[n],
	Pochhammer[u,Abs[u+v]]^Sign[-n] ]*
    v^Sign[-m] (-Pi/Sin[Expand[Pi v]])^Sign[m]
 cond[u_, n_, v_, m_, 2] := m n < 0 && Expand[v/2-u+1/2]===0  
 function[u_, n_, v_, m_, 2] := (2^(v-1) Pi^(-1/2))^Sign[m] *
				gamma[Expand[u-1/2]]^Sign[m]
 cond[u_, n_, v_, m_, 3] := m n < 0 && Expand[u/2-v+1/2]===0 
 function[u_, n_, v_, m_, 3] := (2^(u-1) Pi^(-1/2))^Sign[n] *
				gamma[Expand[v-1/2]]^Sign[n]
 cond[u_, n_, v_, m_, 4] := m n < 0 && Expand[v/2-u]===0  
 function[u_, n_, v_, m_, 4] := (2^(v-1) Pi^(-1/2))^Sign[m] *
				gamma[Expand[u+1/2]]^Sign[m]
 cond[u_, n_, v_, m_, 5] := m n < 0 && Expand[u/2-v]===0 
 function[u_, n_, v_, m_, 5] := (2^(u-1) Pi^(-1/2))^Sign[n] *
				gamma[Expand[v+1/2]]^Sign[n]
 cond[u_, n_, v_, m_, 6] :=  m n > 0 && IntegerQ[2(u-v)] && !IntegerQ[u-v]
 function[u_, n_, v_, m_, 6] := 
    If[ Positive[u-v],
	Pochhammer[v+1/2,Floor[u-v]]^Sign[n] *
	(2^Expand[1-2 v] Pi^(1/2))^Sign[n] gamma[Expand[2 v]]^Sign[n],
	Pochhammer[u+1/2,Floor[v-u]]^Sign[n] *
	(2^Expand[1-2 u] Pi^(1/2))^Sign[n] gamma[Expand[2 u]]^Sign[n]
    ]
 cond[u_, n_, v_, m_, 7] := m n < 0 && IntegerQ[u-v]
 function[u_, n_, v_, m_, 7] := 
    If[ Negative[u-v],
	Factor[Pochhammer[u,v-u]]^If[n>0,-1,1],
	Factor[Pochhammer[v,u-v]]^If[n>0,1,-1] ]

 SameElem[ {v1___,e1_,v2___,e2_,v3___} ] := 
   Join[{e1},SameElem[{v1,v2,v3,e2}]] /; e1===e2

 SameElem[ ___ ] := {}

 ComplementNotSet[ list1_, list2_ ] :=
    If[ list2==={}, list1,
	Join[ Complement[list1,list2],
	      ComplementNotSet[SameElem[list1],SameElem[list2]]
	]
   ]

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

 End[]  (* Simplify`Gamma`Private` *)

 SetAttributes[SimplifyGamma, { ReadProtected, Protected } ];

 End[]   (* System` *)

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