home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / PRELOAD.PAK / GAMMA.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  9.0 KB  |  278 lines
  1.  
  2. (* :Copyright: Copyright 1991  Wolfram Research, Inc.
  3.         Permission is hereby granted to modify and/or make copies of
  4.         this file for any purpose other than direct profit, or as part
  5.         of a commercial product, provided this copyright notice is left
  6.         intact.  Sale, other than for the cost of media, is prohibited.
  7.  
  8.         Permission is hereby granted to reproduce part or all of
  9.         this file, provided that the source is acknowledged.
  10. *)
  11.  
  12. (* :Version: Mathematica 2.1 *)
  13.  
  14. (* :Title: Gamma-functions Simplifications *)
  15.  
  16. (* :Author:
  17.     Victor S. Adamchik, October 1991.
  18. *)
  19.  
  20. (* :Context: Simplify`Gamma` *)
  21.  
  22. (* :Keywords: gamma *)
  23.  
  24. (* :Requirements: none. *)
  25.  
  26. (* :Warnings: none. *)
  27.  
  28. (* :Limitations: none known. *)
  29.  
  30. (* :Summary:
  31.     This package allows to simplify any expressions with 
  32.     Gamma-functions using a sequence of 8 rules which
  33.     apply to the expression. 
  34. *)
  35.  
  36. (*========================================================================*)
  37.  
  38.  Begin["System`"]
  39.  
  40.  SimplifyGamma::usage =
  41. "SimplifyGamma[expr] performs a sequence of gamma-function transformations 
  42.    on expr, and returns the simplest form it finds."
  43.  
  44.  Begin["Simplify`Gamma`Private`"]
  45.  
  46.  
  47. (*========================================================================*)
  48.  
  49.  Unprotect[ SimplifyGamma]
  50.  
  51.  SimplifyGamma[ expr_ ] := SimpGamma[Expand[expr]]
  52.  
  53.  SimpGamma[ expr_Plus/;Not[FreeQ[expr,Gamma[u__]/;Length[{u}]>1]] ] :=
  54.     Module[ {answer,r},
  55.       answer = SimpIncompleteGamma[expr,
  56.             Union[(r = Cases[ expr,a_. Gamma[u_,n_,m_] ])/.
  57.                 b_. Gamma[v_,n_,m_] :> v] ];
  58.       answer = SimpIncompleteGamma1[Expand[answer],
  59.             Union[r/.b_. Gamma[v_,n_,m_] :> Gamma[v,n,m]] ];
  60.       answer /; FreeQ[answer,FailInt]
  61.     ]
  62.  
  63.  SimpGamma[f_Plus] := (SimpGamma[#]&/@f)/;!FreeQ[f,Gamma[_]]
  64.  
  65.  SimpGamma[ expr_ ] := 
  66.    SimpGamma[Numerator[expr]] / SimpGamma[Denominator[expr]] /;
  67.  Head[Denominator[expr]]===Plus && !FreeQ[Denominator[expr],Gamma[_]]  
  68.  
  69.  SimpGamma[ expr_ ] := If[ FreeQ[expr,Gamma[_]], expr, SimpGammaG[expr] ]
  70.  
  71.  SimpIncompleteGamma[ expr_,args_/;Length[args] >1 ] :=
  72.     Module[ {el,min,var},
  73.       el = Union[ args/.n_. + v_. :> v/;NumberQ[n] ];
  74.       If[ Length[el] >1, FailInt, 
  75.           argn = args/.el[[1]]->var;
  76.           min =  Min[argn/.n_. + v_. var :> n/;NumberQ[n]];
  77.           SimpGamma[ (expr/.el[[1]]->var)//.Gamma[k_. + v_. var,n_,m_] :>
  78.           (k+v var-1) Gamma[k+v var-1,n,m] - m^(k+v var-1) E^(-m)/;
  79.           NumberQ[k] && k>min ]/.var->el[[1]] ]
  80.      ]
  81.  
  82.  SimpIncompleteGamma[ expr_,args_ ] := FailInt
  83.  
  84.  SimpIncompleteGamma1[ expr_,{arg1_,arg___} ] :=
  85.     Module[ {r,rnew},
  86.       rnew = 
  87.       Factor[Plus@@((r=Cases[expr,a_. arg1])/.a_. arg1 :> a)] arg1; 
  88.       SimpIncompleteGamma1[ (expr - Plus@@r)+rnew,{arg} ]
  89.     ]
  90.  
  91.  SimpIncompleteGamma1[ expr_, _ ] := expr
  92.  
  93.  SimpGammaG[ Gamma[a_]^n_.] := Gamma[a]^n
  94.  
  95.  SimpGammaG[ expr_Times  ] := 
  96.    Module[ { rr = Cases[expr,Gamma[a_]^n_./;!IntegerQ[n]],
  97.         exprnew, rrr },
  98.      exprnew = expr/(rrr=Times@@rr *Times@@Cases[expr,a_/;FreeQ[a,Gamma[_]]]);
  99.      rrr *
  100.      If[ FreeQ[ exprnew,Gamma[_]] || MatchQ[exprnew,Gamma[a_]^n_.],
  101.      exprnew,
  102.      rrr = Cases[exprnew,Gamma[a_]^n_.];
  103.      rr = exprnew/Times@@rrr;
  104.      If[ rr === exprnew, SimpGamma[#]&/@ exprnew,
  105.          SimpGamma[rr] SimpGammaN[ rrr/.Gamma[a_]^n_.:>{Expand[a],n} ]
  106.      ]
  107.     ]
  108.   ]
  109.  
  110.  SimpGammaG[ f_[expr___] ] := f@@(SimpGamma[#]&/@{expr})
  111.  
  112.  SimpGammaN[ expr_ ] := 
  113.     If[ Length[expr]<2, 
  114.     SimpGamma1[ expr ],
  115.     Block[ { rr =  Cases[expr,{n_Rational+s_.,a_}], SimpGamma1 },
  116.       rr = 
  117.         SimpGamma1[Complement[expr,rr]] * 
  118.         MultiplicationVar[Select[rr,#[[2]]>0&]] *
  119.         MultiplicationVar[#/.{a_,b_}:>{a,-b}&/@Select[
  120.             rr,#[[2]]<0&]]^(-1)//.{
  121.         gamma[a_]^n_ :> gamma[#/.{c_,d_}:>{c,d n}&/@a],
  122.         gamma[a_] gamma[b_] :> gamma[Join[a,b]] };
  123.       rr/.gamma->SimpGamma1/.
  124.         SimpGamma1[a_] SimpGamma1[b_] :> SimpGamma1[Join[a,b]]
  125.     ]
  126.     ] 
  127.  
  128.  SimpGamma1[ {} ] := 1
  129.  
  130.  SimpGamma1[ {v1___, {1,n_}, v2___} ] := SimpGamma1[{v1,v2}]
  131.  
  132.  SimpGamma1[ {v1___, {u_,n_/;Abs[n]>1}, v2___} ] := 
  133.     SimpGamma1[ Join[{v1,v2},Table[{u,Sign[n]}, {i,Abs[n]}]] ]
  134.  
  135.  SimpGamma1[ {v1___, {u1_,n1_}, v2___, {u2_,n2_},v3___ } ] := 
  136.     SimpGamma1[ {v1,v2,v3} ]/;n1+n2==0 && u1===u2
  137.  
  138.  SimpGamma1[ expr_ ] :=
  139.    If[ Length[expr]==1, Gamma[expr[[1,1]]]^expr[[1,2]] ,
  140.    Module[ { new, exprnew=expr, i=1, answer=1 },
  141.      While[ i<=7 && Length[exprnew]!=0,
  142.     answer *=
  143.     Times@@(Fold[ function[#1[[1]], #1[[2]], #2[[1]], #2[[2]], i]&,
  144.               #[[1]], Rest[#] ]&/@(new=FindListCond[exprnew, i]));
  145.     If[ !FreeQ[answer,gamma],True,i+=1];
  146.     exprnew = Join[ ComplementNotSet[exprnew,Flatten[new,1]],
  147.             Flatten[Cases[answer,gamma[a_]^n_.]/.
  148.               gamma[a_]^n_.:>Table[{a,Sign[n]},{i,Abs[n]}],1]
  149.           ];
  150.     answer = answer/.gamma[_]:>1
  151.      ];
  152.      answer Times@@(Power[Gamma[ #[[1]] ],#[[2]]]&/@exprnew) 
  153.    ]
  154.   ]
  155.  
  156.  MultiplicationVar[ expr_ ]:=
  157.     If[ Length[expr]<2, 
  158.     gamma[expr]
  159.     ,
  160.     Module[ { r, rr },
  161.       r = (rr=Cases[expr, {n_Rational + s_,_}])/.
  162.             {n_Rational + s_,_}:> s;
  163.        Times@@(Multiplication[Sort[#]]&/@(
  164.        Cases[rr, {n_Rational + #,_}]&/@Union[r])) *
  165.        Multiplication[ComplementNotSet[expr,rr]]
  166.     ]
  167.    ]
  168.  
  169.  Multiplication[ {} ] := 1
  170.  
  171.  Multiplication[ expr_ ] :=
  172.   If[ Length[expr]<2,
  173.       gamma[expr],
  174.       Module[ {rr = expr/.{a_,b_/;b>1} :> {a,1}, 
  175.            rest, univ, max, flag, var},
  176.     rest = Join[ SameElem[rr], 
  177.         #/.{a_,b_} :>{a,b-1}&/@Select[expr, #[[2]]>1&]];
  178.     rr = Union[rr];
  179.     var = expr[[1,1]]/.{n_Rational+s_. :> s+Floor[n]};
  180.     flag = If[Positive[expr[[1,1]]-var],1,-1]; 
  181.     max = Max[(#/.{s_,_}:>Denominator[s-var])&/@expr];
  182.     univ = Complement[Table[ {i/max flag+var,1}, {i,max-1}], rr];
  183.     If[ max - Length[univ] > Floor[max/2] + 1,
  184.         MultiplicationResult[var,max,rr,flag,rest,univ]
  185.         ,
  186.         If[ GCD@@((#/.{s_,_}:>Denominator[s-var])&/@expr) != 1, 
  187.         var = expr[[1,1]]; 
  188.             max = Max[(#/.{s_,_}:>Denominator[s-var])&/@expr];
  189.         univ = Complement[Table[ {i/max flag+var,1}, {i,max-1}], rr];
  190.         If[ max!=1 && max - Length[univ] > Floor[max/2] + 1,
  191.             MultiplicationResult[var,max,rr,flag,rest,univ]
  192.                 ,
  193.             rr = Cases[expr,{a_/;Denominator[a-var]==max,_}];
  194.             If[rr=!=expr,Multiplication[rr],gamma[rr] ] *
  195.             Multiplication[ 
  196.             Delete[expr,Flatten[Position[expr,#]&/@rr,1]] ]
  197.         ]
  198.             ,
  199.         rr = Cases[expr,{a_/;Denominator[a-var]==max,_}];
  200.         gamma[rr] *
  201.         Multiplication[ Delete[expr,Flatten[Position[expr,#]&/@rr,1]]]
  202.        ]
  203.     ]
  204.      ]
  205.   ]
  206.  
  207.  MultiplicationResult[ var_,max_,rr_,flag_,rest_,univ_ ] :=
  208.     If[ var===0 || (IntegerQ[var] && Negative[var]),
  209.         max^(-var max+1/2) (2 Pi)^((max-1)/2) *
  210.         (-1)^(-var-var max) (-var)!/(max (-max var)!)
  211.         , 
  212.         max^Expand[-var max+1/2] (2 Pi)^((max-1)/2) *
  213.         gamma[{{Expand[max var],1},{var,-1}}]
  214.     ] *
  215.     Multiplication[Join[
  216.         Complement[Join[rr,univ],Table[ {i/max flag+var,1},{i,max-1}]],
  217.              rest]] *
  218.     gamma[univ/.{a_,1} :> {a,-1}]
  219.  
  220.  FindListCond[ {v1___, {u_,n_}, v2___, {v_,m_}, v3___}, c_ ] := 
  221.    Join[ {{{u,n},{v,m}}}, FindListCond[{v1,v2,v3},c] ]/; cond[u,n,v,m,c]
  222.  
  223.  FindListCond[__] := {}
  224.  
  225.  cond[u_, n_, v_, m_, 1] :=  m n > 0 && IntegerQ[u+v] && !IntegerQ[u]
  226.  function[u_, n_, v_, m_, 1] := 
  227.     If[ Positive[u+v],
  228.     Pochhammer[-v,u+v]^Sign[n],
  229.     Pochhammer[u,Abs[u+v]]^Sign[-n] ]*
  230.     v^Sign[-m] (-Pi/Sin[Expand[Pi v]])^Sign[m]
  231.  cond[u_, n_, v_, m_, 2] := m n < 0 && Expand[v/2-u+1/2]===0  
  232.  function[u_, n_, v_, m_, 2] := (2^(v-1) Pi^(-1/2))^Sign[m] *
  233.                 gamma[Expand[u-1/2]]^Sign[m]
  234.  cond[u_, n_, v_, m_, 3] := m n < 0 && Expand[u/2-v+1/2]===0 
  235.  function[u_, n_, v_, m_, 3] := (2^(u-1) Pi^(-1/2))^Sign[n] *
  236.                 gamma[Expand[v-1/2]]^Sign[n]
  237.  cond[u_, n_, v_, m_, 4] := m n < 0 && Expand[v/2-u]===0  
  238.  function[u_, n_, v_, m_, 4] := (2^(v-1) Pi^(-1/2))^Sign[m] *
  239.                 gamma[Expand[u+1/2]]^Sign[m]
  240.  cond[u_, n_, v_, m_, 5] := m n < 0 && Expand[u/2-v]===0 
  241.  function[u_, n_, v_, m_, 5] := (2^(u-1) Pi^(-1/2))^Sign[n] *
  242.                 gamma[Expand[v+1/2]]^Sign[n]
  243.  cond[u_, n_, v_, m_, 6] :=  m n > 0 && IntegerQ[2(u-v)] && !IntegerQ[u-v]
  244.  function[u_, n_, v_, m_, 6] := 
  245.     If[ Positive[u-v],
  246.     Pochhammer[v+1/2,Floor[u-v]]^Sign[n] *
  247.     (2^Expand[1-2 v] Pi^(1/2))^Sign[n] gamma[Expand[2 v]]^Sign[n],
  248.     Pochhammer[u+1/2,Floor[v-u]]^Sign[n] *
  249.     (2^Expand[1-2 u] Pi^(1/2))^Sign[n] gamma[Expand[2 u]]^Sign[n]
  250.     ]
  251.  cond[u_, n_, v_, m_, 7] := m n < 0 && IntegerQ[u-v]
  252.  function[u_, n_, v_, m_, 7] := 
  253.     If[ Negative[u-v],
  254.     Factor[Pochhammer[u,v-u]]^If[n>0,-1,1],
  255.     Factor[Pochhammer[v,u-v]]^If[n>0,1,-1] ]
  256.  
  257.  SameElem[ {v1___,e1_,v2___,e2_,v3___} ] := 
  258.    Join[{e1},SameElem[{v1,v2,v3,e2}]] /; e1===e2
  259.  
  260.  SameElem[ ___ ] := {}
  261.  
  262.  ComplementNotSet[ list1_, list2_ ] :=
  263.     If[ list2==={}, list1,
  264.     Join[ Complement[list1,list2],
  265.           ComplementNotSet[SameElem[list1],SameElem[list2]]
  266.     ]
  267.    ]
  268.  
  269. (*========================================================================*)
  270.  
  271.  End[]  (* Simplify`Gamma`Private` *)
  272.  
  273.  SetAttributes[SimplifyGamma, { ReadProtected, Protected } ];
  274.  
  275.  End[]   (* System` *)
  276.  
  277. (*========================================================================*)
  278.