# Implement the Polygamma function

The Polygamma function of order $$\m\$$, $$\\psi^{(m)}(z)\$$, is the $$\(m + 1)\$$th derivative of the logarithm of the gamma function, which is also the $$\m\$$th derivative of the digamma function. Your task is to take an integer $$\m\$$ and a positive real number $$\z\$$ and output $$\\psi^{(m)}(z)\$$

## Definitions

For those unfamiliar with the functions above (Gamma, Digamma and Polygamma), here are a few different definitions for each:

### $$\\Gamma(z)\$$

• The gamma function is an extension of the factorial ($$\x! = 1\cdot2\cdot3\cdots(x-1)\cdot(x)\$$) to real numbers
• $$\\Gamma(z) = \int_{0}^{\infty}x^{z-1}e^{-x}dx\$$
• $$\\Gamma(n) = (n - 1)! \:,\:\: n \in \mathbb{N}\$$
• $$\\Gamma(n+1) = n\Gamma(n) \:,\:\: n \in \mathbb{N}\$$

### $$\\psi(z)\$$

• The digamma function is the logarithmic derivative of the gamma function
• $$\\psi(z) = \frac{d}{dz}\ln(\Gamma(z))\$$
• $$\\psi(z) = \frac{\Gamma'(z)}{\Gamma(z)}\$$
• $$\\psi(z + 1) = \psi(z) + \frac{1}{z}\$$

### $$\\psi^{(m)}(z)\$$

• The polygamma function of order $$\m\$$ is the $$\m\$$th derivative of the digamma function
• $$\\psi^{(m)}(z) = \frac{d^m}{dz^m}\psi(z)\$$
• $$\\psi^{(m)}(z) = \frac{d^{m+1}}{dz^{m+1}}\ln(\Gamma(z))\$$
• $$\\psi^{(m)}(z+1)= \psi^{(m)}(z) + (-1)^m\frac{m!}{z^{m+1}}\$$

You are to take two inputs, a natural number $$\m\$$ and a positive real number $$\z\$$, and output $$\\psi^{(m)}(z)\$$. The inputs and outputs will always fit within the number bounds of your language, but your algorithm must work theoretically for any and all inputs.

As the output is usually going to be a real number, rather than an integer, the output should be correct to at least 10 significant figures. Trailing zeros may be omitted for exact values. For example, if the output is an integer, trailing decimal 0s are not required, but are allowed if you want.

This is so the shortest code in bytes wins.

## Test cases

Results may differ due to floating point inaccuracies, Python's scipy library was used to generate the values. Values are rounded to 15d.p., unless otherwise stated.

 m,                  z -> ψ⁽ᵐ⁾(z)
17,                  2 -> 1357763223.715975761413574
5,                 40 -> 0.0000002493894351
9,           53.59375 -> 0.00000000001201026493
35,                  9 -> 469354.958166260155849
46,                  5 -> -7745723758939047727202304.000000000000000
7, 1.2222222222222222 -> 1021.084176496877490
28,               6.25 -> -2567975.924144014250487
2,               7.85 -> -0.018426049840992


This table has the values of $$\\psi^{(m)}(z)\$$ for $$\0 \le m \le 9\$$ and $$\1 \le z \le 20\$$:


+---+------------------------+---------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+
|   |           1            |          2          |         3          |         4          |         5          |         6          |         7          |         8          |         9          |         10         |         11         |         12         |         13         |         14         |         15         |         16         |         17         |         18         |         19         |         20         |
+---+------------------------+---------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+
| 0 |   -0.577215664901533   |  0.422784335098467  | 0.922784335098467  | 1.256117668431800  | 1.506117668431800  | 1.706117668431800  | 1.872784335098467  | 2.015641477955610  | 2.140641477955610  | 2.251752589066721  | 2.351752589066721  | 2.442661679975812  | 2.525995013309145  | 2.602918090232222  | 2.674346661660794  | 2.741013328327460  | 2.803513328327460  | 2.862336857739225  | 2.917892413294781  | 2.970523992242149  |
| 1 |   1.644934066848227    |  0.644934066848227  | 0.394934066848226  | 0.283822955737115  | 0.221322955737115  | 0.181322955737115  | 0.153545177959338  | 0.133137014694031  | 0.117512014694031  | 0.105166335681686  | 0.095166335681686  | 0.086901872871768  | 0.079957428427324  | 0.074040268664010  | 0.068938227847684  | 0.064493783403239  | 0.060587533403239  | 0.057127325790783  | 0.054040906037696  | 0.051270822935203  |
| 2 |   -2.404113806319188   |  -0.404113806319189 | -0.154113806319189 | -0.080039732245115 | -0.048789732245114 | -0.032789732245115 | -0.023530472985855 | -0.017699569195768 | -0.013793319195768 | -0.011049834970802 | -0.009049834970802 | -0.007547205368999 | -0.006389797961592 | -0.005479465690312 | -0.004750602716552 | -0.004158010123959 | -0.003669728873959 | -0.003262645625435 | -0.002919710097314 | -0.002628122402315 |
| 3 |   6.493939402266829    |  0.493939402266829  | 0.118939402266829  | 0.044865328192755  | 0.021427828192755  | 0.011827828192755  | 0.007198198563125  | 0.004699239795945  | 0.003234396045945  | 0.002319901304290  | 0.001719901304290  | 0.001310093231071  | 0.001020741379219  | 0.000810664701232  | 0.000654479778283  | 0.000535961259764  | 0.000444408525389  | 0.000372570305061  | 0.000315414383708  | 0.000269374221340  |
| 4 |  -24.886266123440890   |  -0.886266123440879 | -0.136266123440878 | -0.037500691342113 | -0.014063191342113 | -0.006383191342113 | -0.003296771589026 | -0.001868795150638 | -0.001136373275638 | -0.000729931168235 | -0.000489931168235 | -0.000340910050701 | -0.000244459433417 | -0.000179820455575 | -0.000135196191875 | -0.000103591253604 | -0.000080703070010 | -0.000063799959344 | -0.000051098643488 | -0.000041405977726 |
| 5 |  122.081167438133861   |  2.081167438133896  | 0.206167438133897  | 0.041558384635954  | 0.012261509635954  | 0.004581509635954  | 0.002009493175049  | 0.000989510004771  | 0.000531746332896  | 0.000305945162117  | 0.000185945162117  | 0.000118208290511  | 0.000078020533309  | 0.000053159387985  | 0.000037222150950  | 0.000026687171526  | 0.000019534614153  | 0.000014563111016  | 0.000011034967722  | 0.000008484266206  |
| 6 |  -726.011479714984489  |  -6.011479714984437 | -0.386479714984435 | -0.057261607988551 | -0.013316295488551 | -0.004100295488551 | -0.001528279027645 | -0.000654007738836 | -0.000310684984930 | -0.000160150871077 | -0.000088150871077 | -0.000051203486564 | -0.000031109607963 | -0.000019635233198 | -0.000012804988755 | -0.000008590996985 | -0.000005908787970 | -0.000004154139804 | -0.000002978092040 | -0.000002172607350 |
| 7 |  5060.549875237640663  |  20.549875237639476 | 0.862375237639470  | 0.094199654649073  | 0.017295357774073  | 0.004392957774073  | 0.001392271903016  | 0.000518000614207  | 0.000217593204539  | 0.000100511115987  | 0.000050111115987  | 0.000026599144024  | 0.000014877714841  | 0.000008699205352  | 0.000005284083130  | 0.000003317553637  | 0.000002144087193  | 0.000001421585007  | 0.000000964233099  | 0.000000667475582  |
| 8 | -40400.978398747647589 | -80.978398747634884 | -2.228398747634885 | -0.179930526327158 | -0.026121932577158 | -0.005478092577158 | -0.001477178082416 | -0.000478010895205 | -0.000177603485537 | -0.000073530517936 | -0.000033210517936 | -0.000016110901963 | -0.000008296615840 | -0.000004494456155 | -0.000002542957742 | -0.000001494142013 | -0.000000907408791 | -0.000000567407762 | -0.000000364140247 | -0.000000239189714 |
| 9 | 363240.911422382690944 | 360.911422382626938 | 6.536422382626807  | 0.391017718703625  | 0.044948382766125  | 0.007789470766125  | 0.001788099024012  | 0.000503455497598  | 0.000165497161722  | 0.000061424194120  | 0.000025136194120  | 0.000011145599233  | 0.000005284884641  | 0.000002652620244  | 0.000001398085550  | 0.000000768796112  | 0.000000438758675  | 0.000000258758130  | 0.000000157124373  | 0.000000097937278  |
+---+------------------------+---------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+

• Sandbox (deleted). Please don't hesitate to comment any suggestions. – caird coinheringaahing Sep 1 '20 at 13:44
• "An extension of the factorial to real numbers" is not a definition, as there are infinitely many ways to do that. (and the gamma function isn't defined for all real numbers either) – Carla Cvekla Sep 3 '20 at 17:36
• @CarlaCvekla I’m trying to make it understandable what the various functions are, rather than provide a rigorous mathematical definition – caird coinheringaahing Sep 3 '20 at 17:37

# JavaScript (ES7),  68 66 61  59 bytes

Expects (m)(z).

(m,n=m)=>g=z=>n?-n--*g(z):eval("for(k=5e6;k--;)n-=z++**~m")


Try it online!

This is based on the following series representation (from Wikipedia):

$$\psi^{(m)}(z)=(-1)^{m+1}m!\sum_{k=0}^{\infty}\dfrac{1}{(z+k)^{m+1}}$$

### Commented

(m, n = m) =>                // outer function taking m and saving a copy in n
g = z =>                     // inner recursive function taking z
n ?                        // if n is not equal to 0:
-n--                     //   yield -n to invert the sign; decrement n afterwards
* g(z)                   //   multiply by the result of a recursive call
:                          // else:
eval(                    //   evaluate as JS code:
"for(k = 5e6; k--;)" + //     repeat 5 million times:
"n -= z++ ** ~m"     //       subtract z ** -(m+1) from n; increment z
)                        //   end of eval(), which returns the final value of n


# Mathematica, 9 bytes

PolyGamma


Of course Mathematica has a builtin for this

• Here's an answer in Mathematica that doesn't use any gamma-related builtins and is much smarter than this. Go upvote that! – user Sep 2 '20 at 13:48

# R, 8 bytes

psigamma


Try it online!

Takes inputs z,m (test harness stolen from Dominic's answer).

R has a builtin as part of its Special Functions of Mathematics including various forms of the gamma function.

• You should have refreshed your MATLAB / Octave for more golfiness :-P – Luis Mendo Sep 1 '20 at 17:32
• @LuisMendo I did feel particularly smug having beaten Mathematica, too!! – Giuseppe Sep 1 '20 at 17:57

# Mathematica, 32 bytes (30 characters)

Without any Gamma related builtin, uses Bubbler's formula

Sum[#!/(-#2-x)^(#+1),{x,0,∞}]&


Try it online!

# Mathematica, 29 bytes

Without PolyGamma[z] or PolyGamma[n, z]

Log@Gamma@x~D~{x,#+1}/.x->#2&


Try it online!

# Mathematica, 27 bytes

With PolyGamma[z] (this is the equivalent of the digamma function, or $$\\large\psi^0(z)\$$)

PolyGamma@x~D~{x,#}/.x->#2&


Try it online!

A few Mathematica programs that don't use the builtin PolyGamma[n, z].

# 05AB1E, 16 15 bytes

₄nÝ+I±mOI!IÉ·<P


-1 byte thanks to @ovs.

First input is $$\z\$$, second input is $$\m\$$.

Explanation:

Uses the same algorithm as in @Arnauld's JavaScript answer, so make sure to upvote him.
Or to be more precise, it uses the algorithm:

$$\psi^{(m)}(z)=(m\text{%}2\times2-1)\times m!\times\sum_{k=0}^{1000^2}{(z+k)^{\sim m}}$$

₄         # Push 1000
n        # Square it to 1000000
Ý       # Pop and push a list in the range [0,1000000]
+      # Add the first (implicit) input-integer z to each value
I     # Push the second input m
±    # Take it's bitwise-NOT: -m-1
m   # Take each value to the power this -m-1
O  # Sum all values in the list together
I!        # Push the second input m again, and take its factorial
IÉ        # Push the second input m again, and check if it's odd
# (1 if truthy; 0 if falsey)
·       # Double that
<      # And decrease it by 1
P         # And finally take the product of all three values on the stack
# (after which it is output implicitly as result)


NOTE: If there are any very minor inaccuracies in the decimals, the ₄n ($$\1\text{,}000\text{,}000\$$) could be replaced with žm ($$\9\text{,}876\text{,}543\text{,}210\$$), although it would be too slow to run on TIO in that case.

• I>mz can be one byte shorter as I±m. – ovs Sep 1 '20 at 15:19
• @ovs Thanks! I now notice Arnauld had something similar in his answer. – Kevin Cruijssen Sep 1 '20 at 15:29

# APL (Dyalog Unicode) 18.0, 20 bytes

+/!⍤⊣÷(-(⍳!9)+⊢)*1+⊣


Try it online!

-2 bytes thanks to Adám and ngn.

# APL (Dyalog Unicode), 22 bytes

{+/(!⍺)÷(-⍵+⍳1e6)*1+⍺}


Try it online!

Left argument is $$\m\$$, right arg is $$\z\$$.

Uses a slight modification of the formula used by other answers:

\begin{aligned} \psi^{(m)}(z)&=(-1)^{m+1}m!\sum_{k=0}^{\infty}\dfrac{1}{(z+k)^{m+1}} \\ &\approx\sum_{k=0}^{10^6-1}\dfrac{m!}{(-z-k)^{m+1}} \end{aligned}

### How it works

{+/(!⍺)÷(-⍵+⍳1e6)*1+⍺}  ⍝ ⍺←m, ⍵←z
-⍵+⍳1e6        ⍝ vector of -(z+0..999999)
(       )*1+⍺   ⍝ raise each to the power of 1+m
+/(!⍺)÷                ⍝ divide m! by each of above and sum them

• shorter as a train: +/!⍤⊣÷(-(⍳1e6)+⊢)*1+⊣ (but crashes on tio) – ngn Sep 2 '20 at 7:41
• 1e6 can be !9 – Adám Sep 3 '20 at 11:28

# Octave / MATLAB, 4 bytes

@psi


Try it online!

# R, 525145 44 bytes

Edit: -1+1 bytes thanks to Giuseppe, who also pointed out that there is already a built-in R function, psigamma, that solves the task for only 8 bytes

Edit2: ...and -6 more bytes thanks to Robin Ryder

function(m,z)gamma(M<-m+1)*sum((-z:-1e4)^-M)


Try it online!

Uses the same formula as Arnauld's answer.
Series representations like this are very well-suited to R as a natively-vectorized language.

Change the 1e4 to higher values (up to 9e9 without increasing the byte count) for progressively higher-accuracy and slower runtime.

• There is also of course the psigamma function. – Giuseppe Sep 1 '20 at 16:20
• And function(m,z)gamma(M<-m+1)*(-1)^M*sum((z+0:1e4)^-M) is a byte shorter – Giuseppe Sep 1 '20 at 16:26
• You can even increase the bound to 9e9 without adding bytes. :-) – Robin Ryder Sep 1 '20 at 20:19
• Since (-1)^M == (-1)^-M, you can push into the sum for -6 bytes: Try it online! – Robin Ryder Sep 1 '20 at 21:33
• 44 bytes -- it'll lose precision slightly since it'll sum ceiling(z) fewer terms than the way you have it implemented. – Giuseppe Sep 2 '20 at 16:42

# Python 3, 38 bytes

from scipy.special import*
f=polygamma


Try it online!

• You can remove 2 bytes by not naming it, even though it does make testing it slightly weird: Try it online! – caird coinheringaahing Sep 1 '20 at 23:07
• @cairdcoinheringaahing I know that's the general consensus, but I've decided it's just too stange (you could take it to the max and not have any bytes because it's just there somewhere) - thanks anyway! :-) – Noodle9 Sep 1 '20 at 23:09

# Java, 168148 102 bytes

(m,z)->{double p=1-m%2*2,f=0;long i=m;for(;i>0;)p*=i--;for(;i<1e7;)f-=p*Math.pow(z+i++,~m);return f;};


## Explanation

For convenience, here the auto-formatted version:

(m, z) -> {
double p = 1 - m % 2 * 2, f = 0;
long i = m;
for (; i > 0; ) p *= i--;
for (; i < 1e7; ) f -= p * Math.pow(z + i++, ~m);
return f;
};


So typical Java code: quite verbose. At least my version.

Edit: could save 20 bytes thanks to @user

Edit: saved even more bytes thanks to @ceilingcat

Try it online

• Try double p(long m,double z){long i,x=i=1;for(;i<m;)x*=++i;double pow=Math.pow(-1,m+1)*x;double f=0;for(int k=0;k<10000000;++k)f+=pow/Math.pow(z+k,m+1);return f;} (159 bytes) – user Sep 3 '20 at 13:38

# Pyth, 24 bytes

**^_1JhhQ*FhQsm^+deQ_JCG


Try it online! (link points to slightly different code that sums 1e5 terms instead of 1.56e62 terms to make code runnable and avoid overflow errors)

## Explanation

Uses the same algorithm as in @Arnauld's JavaScript answer, so make sure to upvote him.

**^_1JhhQ*FhQsm^+deQ_JCG
JhhQ                  : Set J to first input + 1
^_1J                     : -1 ^ J
*       *FhQ              : times factorial of first input
*            s             : times sum of
m            : mapping
^+deQ_J     :   F(d): (d + (second input)) ^ -J
CG   : on range(1.56e62)


# Scala, 68 bytes

Saved 3 bytes and fixed my answer thanks to Arnauld

Uses the algorithm from Arnauld's answer

m=>z=>(0 to 1<<20 map(z+_ pow ~m)sum)*(m%2*2-1)*(1.0/:(1 to m))(_*_)


Try it online!

• I don't know Scala and there might be a better trick, but you can save at least one byte by using 1<<20 instead of 999999. – Arnauld Sep 1 '20 at 16:18