home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / PRELOAD.PAK / POLYGAMM.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  10.2 KB  |  328 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: PolyGamma-functions Simplifications *)
  15.  
  16. (* :Author:
  17.     Victor S. Adamchik, October 1991.
  18. *)
  19.  
  20. (* :Context: Simplify`PolyGamma` *)
  21.  
  22. (* :Keywords: polygamma *)
  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.     PolyGamma-functions using a sequence of 10 rules which
  33.     apply to the expression. 
  34. *)
  35.  
  36. (*========================================================================*)
  37.  
  38.  Begin["System`"]
  39.  
  40.  SimplifyPolyGamma::usage =
  41. "SimplifyPolyGamma[expr] performs a sequence of polygamma-function 
  42.    transformations on expr, and returns the simplest form it finds."
  43.  
  44.  Begin["Simplify`PolyGamma`Private`"]
  45.  
  46. (*========================================================================*)
  47.  
  48.  Unprotect[ SimplifyPolyGamma]
  49.  
  50.  SimplifyPolyGamma[ expr_ ] := SimpPolyGamma[expr//Simplify]
  51.  
  52.  SimpPolyGamma[f_Times] := SimpPolyGamma[#]&/@f
  53.  
  54.  SimpPolyGamma[ expr_/;!FreeQ[Denominator[expr],PolyGamma[___]] ] := 
  55.    SimpPolyGamma[Numerator[expr]] / SimpPolyGamma[Denominator[expr]]  
  56.  
  57.  SimpPolyGamma[ expr_ ] := 
  58.    SimpPolyGammaG[ Expand[expr/.PolyGamma[v_] :> PolyGamma[0,v]] ]
  59.  
  60.  SimpPolyGammaG[ expr_/;FreeQ[expr,PolyGamma[___]] ] := expr
  61.  
  62.  SimpPolyGammaG[ PolyGamma[a_,b_]^n_.] := SimpPolyGammaN[{{{a,b},1}}]^n
  63.  
  64.  SimpPolyGammaG[ expr_Plus  ] := 
  65.    Block[ { PolyGamma, answer,exprnew,
  66.         rr = Cases[expr,s_. PolyGamma[a__]^n_], rrr},
  67.      exprnew = expr-Plus@@rr-
  68.       (rrr=Plus@@Cases[expr,a_/;FreeQ[a,PolyGamma[__]]]);
  69.      answer =  Plus@@(SimpPolyGamma[#]&/@rr) + rrr;
  70.      answer +=
  71.          If[ FreeQ[ exprnew,PolyGamma[__]],
  72.          exprnew,
  73.          rrr = If[ MatchQ[exprnew, s_. PolyGamma[a__]],
  74.                {exprnew},
  75.                Cases[exprnew,s_. PolyGamma[a__]]
  76.            ]; 
  77.          rr = exprnew-Plus@@rrr;
  78.          If[ rr === exprnew, SimpPolyGamma[exprnew],
  79.              SimpPolyGamma[exprnew-Plus@@rrr] + 
  80.              SimpPolyGammaN[ rrr/.
  81.             PolyGamma[a_,b_] n_. :> If[ FreeQ[n,PolyGamma],
  82.             {{a,Expand[b]},n}, FailS] ]
  83.         ]
  84.         ]; 
  85.     If[ FreeQ[answer, FailS], answer, expr ] 
  86.   ]
  87.  
  88.  SimpPolyGammaG[ f_[expr___] ] := f@@(SimpPolyGamma[#]&/@{expr})
  89.  
  90.  SimpPolyGammaN[ expr_ ] := FailS/; !FreeQ[expr,FailS]
  91.  
  92.  SimpPolyGammaN[ expr_ ] := 
  93.    Module[ { rr = expr, i = 0, answer = 0, max },
  94.      max = Max[(#/.{{a_,_},_} :> a)&/@expr];
  95.      While[ i <= max,
  96.     answer += 
  97.       Multiplication[new = Cases[rr, {{i,_},_}]];
  98.       rr = Complement[rr, new]; i+=1
  99.      ];
  100.      answer
  101.   ]
  102.  
  103.  Multiplication[ expr_ ] := 
  104.     If[ Length[expr]<2, 
  105.     SimpPolyGamma1[expr]
  106.     ,
  107.     Block[ { SimpPolyGamma1 },
  108.      (MultiplicationVar[Cases[expr,{{_,a_/;!NumberQ[a]},_}],False] //.
  109.      SimpPolyGamma1[a_] + SimpPolyGamma1[b_] :>SimpPolyGamma1[Join[a,b]])
  110.      + 
  111.      (MultiplicationVar[Cases[expr,{{_,a_/;NumberQ[a]},_}],True] //.
  112.      SimpPolyGamma1[a_] + SimpPolyGamma1[b_] :>SimpPolyGamma1[Join[a,b]])
  113.     ]
  114.    ]
  115.  
  116.  MultiplicationVar[ expr_, const_ ] := 
  117.     If[ Length[expr]<2, 
  118.     SimpPolyGamma1[expr]
  119.     ,
  120.     Module[ { r, rr },
  121.        If[ const,
  122.         r = (#/.{{_,n_Rational},d_}:> 
  123.             {0,d})&/@(rr=Cases[expr, {{_,n_Rational},_}]),
  124.         r = (#/.{{_,n_Rational+s},d_}:> 
  125.             {s,d})&/@(rr=Cases[expr, {{_,n_Rational+s},_}])
  126.        ];
  127.        Plus@@(MultiplicationVar1[Sort[#]]&/@(
  128.        Cases[rr, {{_,n_Rational + #[[1]]},#[[2]]}]&/@Union[r])) +
  129.        SimpPolyGamma1[Complement[expr,rr]]
  130.     ]
  131.    ]
  132.  
  133.  MultiplicationVar1[ expr_ ] := 
  134.     If[ Length[expr]<2, 
  135.     SimpPolyGamma1[expr]
  136.     ,
  137.     Module[ { var, flag, max, univ, p=expr[[1,1,1]], r=expr[[1,2]] },
  138.        var = expr[[1,1,2]]/.{n_Rational+s_. :> s+Floor[n]};
  139.        flag = If[Positive[expr[[1,1,2]]-var],1,-1]; 
  140.        max = Max[(#/.{{_,s_},_}:>Denominator[s-var])&/@expr];
  141.        univ = Complement[Table[ 
  142.         {{p,i/max flag+var},r},{i,max-1}],expr];
  143.        If[ max - Length[univ] > Floor[max/2] + 1,
  144.         If[IntegerQ[var] && Negative[var],
  145.            r ( If[ p==0, (1-max) EulerGamma,
  146.                (max^(p+1)-1) (-1)^(p+1) p! Zeta[p+1] 
  147.             ] +
  148.                max^(p+1) p! Sum[1/i^(p+1),{i,-max var}] -
  149.                p! Sum[1/i^(p+1),{i,-var}] - 
  150.                If[p==0, max Log[max], 0] )
  151.            ,
  152.            If[ flag == 1,
  153.             SimpPolyGamma1[{{{p,Expand[max var]}, max^(p+1) r}, 
  154.                     {{p,var},-r} }] -
  155.             If[p==0, r max Log[max], 0],
  156.             SimpPolyGamma1[{{{p,Expand[max var-max]}, max^(p+1) r}, 
  157.                     {{p,var-1},-r} }] -
  158.             If[p==0, r max Log[max], 0]
  159.            ]
  160.         ]  +
  161.         SimpPolyGamma1[ Join[Complement[ Join[expr,univ],
  162.                 Table[{{p,i/max flag+var},r}, {i,max-1}]
  163.                  ], (#/.{{a_,b_},c_}:>{{a,b},-c})&/@univ ]]
  164.         ,
  165.         rr = Cases[Union[expr],{{_,a_/;Denominator[a-var]==max},_}];
  166.          SimpPolyGamma1[rr] +
  167.         MultiplicationVar1[ 
  168.             Delete[expr,Flatten[Position[expr,#]&/@rr,1]]]
  169.        ]
  170.     ]
  171.    ]
  172.  
  173.  SimpPolyGamma1[ {} ] := 0
  174.  
  175.  SimpPolyGamma1[ expr_ ] :=
  176.    Module[ { new, exprnew=expr, i=1, answer=0 },
  177.      While[ i<=3 && Length[exprnew]!=0,
  178.     answer +=
  179.     If[ i==3, Plus@@(function[#[[1]],i]&/@(new=FindListCond[exprnew, i])),
  180.         Plus@@(Fold[ function[#1, #2, i]&,
  181.               #[[1]], Rest[#] ]&/@(new=FindListCond[exprnew, i])) ];
  182.     If[ !FreeQ[answer,polygamma]&&i!=3,True,i+=1];
  183.     exprnew = Join[ ComplementNotSet[exprnew,Flatten[new,1]],
  184.             If[MatchQ[answer,s_. polygamma[a__]],
  185.                {answer},
  186.                Cases[answer,s_. polygamma[a__]]
  187.             ]/.s_. polygamma[w1_,w2_] :>{{w1,w2},s}
  188.           ];
  189.     answer = answer/.polygamma[__]:>0
  190.      ];
  191.      SimpLog[answer] + Plus@@(Times[PolyGamma@@#[[1]], #[[2]]]&/@exprnew) 
  192.    ]
  193.  
  194.  FindListCond[ {v1___, {{p_,u_},n_},v2___ }, 3 ] :=  
  195.    Join[ {{{{p,u},n}}}, FindListCond[{v1,v2}, 3] ] /; 
  196.  IntegerQ[u] || RationalQ[u]
  197.  
  198.  FindListCond[ {v1___, {{p_,u_},n_}, v2___, {{p_,v_},m_}, v3___},c_/;c!=3] := 
  199.    Join[ {{{{p,u},n},{{p,v},m}}}, FindListCond[{v1,v2,v3}, c] ] /; 
  200.  cond[{{p,u},n},{{p,v},m}, c]
  201.  
  202.  FindListCond[___] := {}
  203.  
  204.  cond[{{p_,u_/;!IntegerQ[u]},n_},{{p_,v_},m_}, 1] := 
  205.     IntegerQ[u+v] && If[ EvenQ[p], Znak[m n], !Znak[m n]]
  206.  
  207.  cond[{{p_,u_},n_},{{p_,v_},m_}, 2] := IntegerQ[u-v] && Negative[m n]
  208.  
  209.  function[{{p_,u_},n_},{{p_,v_},m_}, 1] := 
  210.     If[ n+m===0, 
  211.     If[ Positive[u+v],
  212.         -m p! Sum[ 1/(u-k)^(p+1), {k, 0, u+v-1} ],
  213.         m p! Sum[ 1/(u+k)^(p+1), {k, 1, -u-v} ] 
  214.     ] -
  215.     n Pi Module[ { y }, Pi^p D[Cot[y],{y,p}]/.y ->Expand[Pi u] ] -
  216.     (-1)^p n p!/u^(p+1)
  217.     ,
  218.     If[ FreeQ[ {n,m},Complex] && NumberQ[n] && NumberQ[m],
  219.         If[ Positive[u+v],
  220.             -m p! Sum[ 1/(u-k)^(p+1), {k, 0, u+v-1} ],
  221.             m p! Sum[ 1/(u+k)^(p+1), {k, 1, -u-v} ] 
  222.         ] -
  223.         Sign[n] Min[Abs[n], Abs[m]]*
  224.         ( Pi Module[ { y }, Pi^p D[Cot[y], {y,p}]/.y->Expand[Pi u] ] +
  225.           (-1)^p p!/u^(p+1)) + 
  226.           (n-Sign[n] Min[Abs[n], Abs[m]]) polygamma[p,u] +
  227.           Sign[m] (Abs[m]-Min[Abs[n], Abs[m]]) polygamma[p,-u],
  228.         n PolyGamma[p,u] + m PolyGamma[p,v]
  229.     ]
  230.     ]
  231.  
  232.  function[{{1,1/4},n_}, 3] := n(Pi^2 + 8 Catalan)
  233.  function[{{1,3/4},n_}, 3] := n(Pi^2 - 8 Catalan) 
  234.  function[{{k_,1},n_}, 3]  := (-1)^(k+1) k! n Zeta[k+1]
  235.  function[{{k_,1/2},n_},3] := (-1)^(k+1) k! n Zeta[k+1] (2^(k+1)-1)
  236.  function[{{0,k_/;0<k<1},n_},3]   :=  
  237.     n(-EulerGamma - Log[2 Denominator[k]] -
  238.     Pi Cos[Pi k]/(2 Sin[Pi k]) + 
  239.     2 Sum[Cos[2 k i Pi] Log[Sin[Pi i/Denominator[k]]],
  240.              {i,1,Floor[(Denominator[k]-1)/2]}])
  241.  function[{{p_,k_Integer},n_},3] :=  
  242.    (-1)^p p! n (-Zeta[p+1]+Sum[i^(-p-1),{i,k-1}])
  243.  function[{{p_,k_/;k>1 || k<0},n_},3]   :=  
  244.    If[ k>1, n(-1)^p p!(k-1)^(-p-1) + function[{{p,k-1},n},3],
  245.         -n(-1)^p p!k^(-p-1) + function[{{p,k+1},n},3]
  246.    ]/; p==0 || Denominator[k]==2 || Denominator[k]==4
  247.  function[{{p_,k_},n_},3]   :=  n polygamma[p,k]
  248.  function[{{p_,u_},n_},{{p_,v_},m_}, 2] := 
  249.     If[ Positive[u-v],
  250.     n p! (-1)^p Sum[ 1/(v+k)^(p+1), {k, 0, u-v-1} ] +
  251.     (n+m) polygamma[p,v],
  252.     m p! (-1)^p Sum[ 1/(u+k)^(p+1), {k, 0, v-u-1} ] +
  253.     (n+m) polygamma[p,u]
  254.     ]
  255.  
  256.  ComplementNotSet[ list1_, list2_ ] :=
  257.     If[ list2==={}, list1,
  258.     Join[ Complement[list1,list2],
  259.           ComplementNotSet[SameElem[list1],SameElem[list2]]
  260.     ]
  261.    ]
  262.  
  263.  SameElem[ {v1___,e1_,v2___,e2_,v3___} ] := 
  264.    Join[{e1},SameElem[{v1,v2,v3,e2}]] /; e1===e2
  265.  
  266.  SameElem[ ___ ] := {}
  267.  
  268.  RationalQ[ s_ ] := If[ Head[s]===Rational, True, False ]
  269.  
  270.  Znak[ s_?Negative ] := True  
  271.  Znak[ s_?Negative w_ ] := True  
  272.  Znak[ _ ] := False
  273.  
  274.  SimpLog[ expr_/; !FreeQ[expr,Log[_]] ] :=
  275.    Module[ {exprN = Expand[expr],list,div,pos = {},i=0 },
  276.      exprN = exprN//.{
  277.     Log[a_Power] :> a[[2]] Log[a[[1]]] /;NumberQ[N[a]],
  278.     Log[a_ s_./;NumberQ[N[a]]] :> 
  279.         Log[Numerator[a]]-Log[Denominator[a]]+Log[s]
  280.     };
  281.      list = Union[Cases[exprN,w_. Log[n_Integer]/;FreeQ[w,Log[_]]]/.
  282.                       w_. Log[n_] :> n];
  283.      div = GCD@@list;
  284.      If[ Length[list]>1 && div!=1, 
  285.          list = Complement[list,{div}];
  286.          SimpLog[exprN//.BuildRule[
  287.          Log[#]&/@list ,Log[div]+Log[#]&/@(list/div) ]]
  288.     ,
  289.          If[ Length[list]<=2,
  290.              exprN
  291.          ,
  292.              While[Length[pos]==0 && Length[list]>1, 
  293.               i = i+ 1;
  294.               listN = {list[[1]],#}&/@Rest[list];
  295.               list = Rest[list];
  296.               pos = Position[ (GCD@@#)&/@listN,a_/; a != 1 ]
  297.              ];
  298.              If[ Length[pos]!=0,
  299.                  list = Take[listN,pos[[1]]][[1]];
  300.                  div = GCD@@list; 
  301.                  SimpLog[exprN//.BuildRule[Log[#]&/@list ,      
  302.                          Log[div]+Log[#]&/@(list/div) ]]
  303.          ,
  304.                  exprN
  305.                ]
  306.           ]
  307.      ]
  308.    ]    
  309.  
  310.  SimpLog[ expr_ ] := expr
  311.  
  312.  BuildRule[{l_,lR___},{r_,rR___}] :=
  313.   Join[{l :> r},BuildRule[{lR},{rR}]]
  314.  
  315.  BuildRule[{},{}] := {}
  316.  
  317.  
  318. (*========================================================================*)
  319.  
  320.  End[]  (* Simplify`PolyGamma`Private` *)
  321.  
  322.  SetAttributes[SimplifyPolyGamma, { ReadProtected, Protected } ];
  323.  
  324.  End[]   (* System` *)
  325.  
  326. (*========================================================================*)
  327.  
  328.