MathGroup Archive 2003

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: New version, new bugs

  • To: mathgroup at smc.vnet.net
  • Subject: [mg43110] Re: New version, new bugs
  • From: Maxim <dontsendhere@.>
  • Date: Tue, 12 Aug 2003 04:43:22 -0400 (EDT)
  • References: <bgq9q4$d50$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Mathematica gets the limit of the hypergeometric function wrong when
parameters approach negative integer values (also of the beta function and
probably some others). Here's the code to do this task for the
hypergeometric and beta functions, maybe someone will find it useful:

ClearAll[lim]
Options[lim] = {Direction -> Automatic};
lim[$expr_, n_Symbol -> n0_, opt___?OptionQ] := Module[
    {ser, expr=$expr, ord = 1, ans},

    ser[La_, Lb_, z_] := Module[
        {ck0, ck, k0, k, zz, ans},
        If[Length[La] > Length[Lb] + 1,
          Return[Indeterminate] ];
        k0 = Max[-Select[Limit[#, n -> n0, opt], IntegerQ] & /@ {La, Lb},
0];
        ck = Times @@ Pochhammer[La, k]/Times @@ Pochhammer[Lb, k];
        ans = Sum[ck*z^k/k!, {k, 0, k0 + 1}];
        ck0 = ck /. k -> k0 + 1;
        ck = Times @@ Pochhammer[La + k0 + 1, k - k0 - 1]/
            Times @@ Pochhammer[Lb + k0 + 1, k - k0 - 1];
        ck += O[n, n0]^ord;
        If[Head[ck] =!= SeriesData,
          Return[Indeterminate] ];
        ans += ck0*MapAt[Sum[#*zz^k/k!, {k, k0 + 2, Infinity}] & /@ # &, ck,

        3] /. zz -> z
        ];

    expr = expr /.
        Beta[z1_, z2_, a_, b_] -> Beta[z2, a, b] - Beta[z1, a, b] /.
      Beta[z_:1, a_, b_] -> z^a/a*Hypergeometric2F1[a, 1 - b, a + 1, z];
    While[
      ans = expr /.
          {Hypergeometric0F1[a_, z_] :> ser[{}, {a}, z],
            Hypergeometric0F1Regularized[a_, z_] :>
              ser[{}, {a}, z]/Gamma[a],
            Hypergeometric1F1[a_, b_, z_] :> ser[{a}, {b}, z],
            Hypergeometric1F1Regularized[a_, b_, z_] :>
              ser[{a}, {b}, z]/Gamma[b],
            Hypergeometric2F1[a_, b_, c_, z_] :> ser[{a, b}, {c}, z],
            Hypergeometric2F1Regularized[a_, b_, c_, z_] :>
              ser[{a, b}, {c}, z]/Gamma[c],
            HypergeometricPFQ[a_, b_, z_] :> ser[a, b, z],
             HypergeometricPFQRegularized[a_, b_, z_] :>
               ser[a, b, z]/Times @@ Gamma[b]};
      Head[ans] === SeriesData && ans[[4]] == ans[[5]] == 0,
      ord++
      ];
    ans = Limit[Normal[ans], n -> n0, opt];
    ans = ans /. Gamma[a_ /; IntegerQ[a] && Positive[a], z_] :>
          (a - 1)!*E^-z*Sum[z^k/k!, {k, 0, a - 1}];
    ans = Collect[ans, E^_ | _[_], Apart];
    ans /; ans =!= Indeterminate
    ]


Example:

In[4]:=
lim[Hypergeometric2F1[n, 2*n, (n + 1)^2 - 3, x], n -> -1]

Out[4]=
1 - (8*x)/3 + x^2/3 + (-2 + (4*x)/3)*Log[1 - x]


Now here's the list of errors I encountered while writing these 50 lines of
code:

1) when writing out sums of rising factorials:

In[1]:=
Module[
  {c = Pochhammer[-2, -1 + k]},
  Sum[c,{k,1,Infinity}]
]

Out[1]=
Infinity*Pochhammer[-2, -1 + k]

in version 4.2 this sum is left unevaluated -- hardly an improvement;

2) when testing cases where the result involves sums of digamma functions:

In[1]:=
Sum[PolyGamma[0, k]/(2^k*k!), {k, 1, Infinity}]

Out[1]=
(1/2)*(2*EulerGamma - 2*Sqrt[E]*EulerGamma -
3*E^(1/4)*Sqrt[2*Pi]*BesselI[-(1/2), 1/4]*Log[2] -
   E^(1/4)*Sqrt[2*Pi]*Derivative[1, 0][BesselI][-(1/2), 1/4])

correct value is 2*EulerGamma - ExpIntegralEi[1/2] + Sqrt[E]*Gamma[0, 1/2] -
Log[2] - Sqrt[E]*Log[2]; this sum was also left unevaluated in 4.2;

3) at first I was scanning the lists of coefficients in a more procedural
way, and noticed a strange behaviour or Increment; here's an example to
demonstrate it more clearly:

In[1]:=
Compile[
    {},
    Module[
      {L={10,20,30},a=1},
      L[[a++]]++;
      Append[L,a]
      ]
    ][]

Out[1]=
{10,31,30,4}

Imagine using a C compiler with this sort of glitches. Once again the bug
appears in a function that doesn't have a precise definition; all this
certainly would have been of little importance for an educational software,
but not for a serious development tool.

Also I once again ran into the fact that named and unnamed (Slot[k])
function parameters behave differently; consider the following definition:

In[1]:=
f = #1 /. {x_, y_} :> Refine[Sign[{x, y}], #2] &;
f[{a - b, b - a}, a > b]

Out[2]=
{1, -1}

The function f simply returns the signs of the pair of elements under given
assumptions, but not in this case:

In[3]:=
f[{y, x}, y > 0 && x < 0]

Out[3]=
{-1, 1}

Ironically, I wrote about this issue myself before, and I still keep falling
for it; it's such an intuitive language. This brings up an interesting
consideration: are such quirks worth retaining merely for the sake of
compatibility?

Maxim Rytin
m.r at prontomail.com



  • Prev by Date: Re: Palette Question
  • Next by Date: [newbie] can't get the answer
  • Previous by thread: Re: Re: New version, new bugs
  • Next by thread: Re: Re: New version, new bugs