# Expanding PolyGamma function error

Bug introduced in 6.0 or earlier and fixed in 10.3.0

For large half-integer arguments of the PolyGamma[] function, FunctionExpand[] or FullSimplify[] return strange, error-like, expressions

FunctionExpand[PolyGamma[0, 1051/2]]
(* => ... + SystemPolyGammaDumpres$14430 *)  or even (* Hold[RuleCondition[$ConditionHold[$ConditionHold[ SystemPolyGammaDumpres$5680]],
Head[SystemPolyGammaDumpres$5680] =!= DirectedInfinity]] *)  These, combined further with simple algebraic expressions like +, give FunctionExpand[PolyGamma[0, 1051/2] + 1] (* => ... + Hold[Fail] *)  (why Hold[Fail]? and not some message?) This "error" does not appear when I decrease the argument of PolyGamma[]: FunctionExpand[PolyGamma[0, 21/2]] (* 62075752/14549535 - EulerGamma - Log[4] *)  Is this a bug of my Mathematica version (10.0.2.0 for Mac OS X)? Recurrence identities for polygamma should be implemented in Mathematica (compare output of PolyGamma[0, 10/2 + 1] - PolyGamma[0, 10/2] and PolyGamma[0, 2001/2 + 1] - PolyGamma[0, 2001/2]); how can one trace that bug and possibly overcome it? • Gamma[] is just evaluated, for HarmonicNumber[] I do get the same as for PolyGamma[]. – mmal Sep 2, 2015 at 14:05 • A clear bug, I think; somebody botched the recursive relation behind the scenes. Can you try Block[{$RecursionLimit = Infinity}, FunctionExpand[PolyGamma[1051/2]]]? Sep 2, 2015 at 14:08
• This seems to go back a while. I've filed a bug report with the developers. Sep 2, 2015 at 15:00
• Sorry, I haven't mentioned that I'm interested in exact values and a way to overcome this bug not in approximate (numerical) results.
– mmal
Sep 2, 2015 at 15:16
• Fixed in the development version today. Sep 2, 2015 at 23:21

As noted by ilian, this has been fixed in the newest version. For those stuck with older versions, a workaround can be devised based on this expansion formula, like so:

a = 1051/2;
n = IntegerPart[a]; f = FractionalPart[a];
PolyGamma[f] + Sum[1/(f + k), {k, 0, n - 1}, Method -> "Procedural"]
// FunctionExpand


This works even for the larger a = 4601/2 case, unlike the suggestion I gave in the comments.

As previously mentioned in the comments, this bug has been fixed as of Mathematica 10.3.

\$Version

(* "10.3.0 for Mac OS X x86 (64-bit) (October 9, 2015)" *)

FunctionExpand[PolyGamma[0, 1051/2]]

(* 1285419799936525483590440356485913118557238953850917571655841\
5741182979315173079330077832682124967619286005471667823256760099546705\
2891030373156621089354809760128773005810025097953873373875179741234456\
9187486859122888279655863758285716796306889323506583432571529367704503\
6137729023370199592517121547190898489561519446980549309597633794551282\
8024779606014776869981600684213070729430041879996083178849511286173189\
91345788279487197852589070073000030156838138/
15624578896737929535308911953771213336402116274975379996559013353031\
1680190991120771924495781500410869490340481605874727401745947540728422\
2556840660803953943256878960592583459598665600326021616974229885783431\
2790170314522959586741781679470734395597442825041097450304608994258061\
5845822338517884531761963978132838586622754530461858355206268627504641\
8438386432184286804056434659183179146119622607050759848035815087097489\
678659193191774142898224573450093125 - EulerGamma - Log[4] *)

FunctionExpand[HarmonicNumber[1051/2 - 1] - EulerGamma]

(* 12854197999365254835904403564859131185572389538509175716558415\
7411829793151730793300778326821249676192860054716678232567600995467052\
8910303731566210893548097601287730058100250979538733738751797412344569\
1874868591228882796558637582857167963068893235065834325715293677045036\
1377290233701995925171215471908984895615194469805493095976337945512828\
0247796060147768699816006842130707294300418799960831788495112861731899\
1345788279487197852589070073000030156838138/
15624578896737929535308911953771213336402116274975379996559013353031\
1680190991120771924495781500410869490340481605874727401745947540728422\
2556840660803953943256878960592583459598665600326021616974229885783431\
2790170314522959586741781679470734395597442825041097450304608994258061\
5845822338517884531761963978132838586622754530461858355206268627504641\
8438386432184286804056434659183179146119622607050759848035815087097489\
678659193191774142898224573450093125 - EulerGamma - Log[4] *)

FunctionExpand[PolyGamma[0, 9999/2]] // N

(* 8.51699 *)

• Is HarmonicNumber[] also covered by this fix? Oct 17, 2015 at 3:02
• @J.M. Yes, example added. Oct 17, 2015 at 4:45