17
\$\begingroup\$

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}}\$

Task

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 -> ψ(m)(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 |
+---+------------------------+---------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+
asked Sep 1, 2020 at 13:43
\$\endgroup\$
3
  • \$\begingroup\$ Sandbox (deleted). Please don't hesitate to comment any suggestions. \$\endgroup\$ Commented Sep 1, 2020 at 13:44
  • \$\begingroup\$ "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) \$\endgroup\$ Commented Sep 3, 2020 at 17:36
  • \$\begingroup\$ @CarlaCvekla I’m trying to make it understandable what the various functions are, rather than provide a rigorous mathematical definition \$\endgroup\$ Commented Sep 3, 2020 at 17:37

12 Answers 12

16
\$\begingroup\$

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
answered Sep 1, 2020 at 14:17
\$\endgroup\$
11
\$\begingroup\$

Mathematica, 9 bytes

PolyGamma

Of course Mathematica has a builtin for this

answered Sep 1, 2020 at 16:21
\$\endgroup\$
1
  • 3
    \$\begingroup\$ Here's an answer in Mathematica that doesn't use any gamma-related builtins and is much smarter than this. Go upvote that! \$\endgroup\$ Commented Sep 2, 2020 at 13:48
9
\$\begingroup\$

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].

answered Sep 1, 2020 at 19:56
\$\endgroup\$
8
\$\begingroup\$

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.

answered Sep 1, 2020 at 17:01
\$\endgroup\$
2
  • 1
    \$\begingroup\$ You should have refreshed your MATLAB / Octave for more golfiness :-P \$\endgroup\$ Commented Sep 1, 2020 at 17:32
  • 2
    \$\begingroup\$ @LuisMendo I did feel particularly smug having beaten Mathematica, too!! \$\endgroup\$ Commented Sep 1, 2020 at 17:57
7
\$\begingroup\$

05AB1E, (削除) 16 (削除ここまで) 15 bytes 

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

-1 byte thanks to @ovs.

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

Try it online or verify all test cases.

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}}$$

4 # 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 4n (\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.

answered Sep 1, 2020 at 15:04
\$\endgroup\$
2
  • 1
    \$\begingroup\$ I>mz can be one byte shorter as I±m. \$\endgroup\$ Commented Sep 1, 2020 at 15:19
  • \$\begingroup\$ @ovs Thanks! I now notice Arnauld had something similar in his answer. \$\endgroup\$ Commented Sep 1, 2020 at 15:29
6
\$\begingroup\$

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
answered Sep 2, 2020 at 1:13
\$\endgroup\$
2
  • \$\begingroup\$ shorter as a train: +/!⍤⊣÷(-(⍳1e6)+⊢)*1+⊣ (but crashes on tio) \$\endgroup\$ Commented Sep 2, 2020 at 7:41
  • \$\begingroup\$ 1e6 can be !9 \$\endgroup\$ Commented Sep 3, 2020 at 11:28
5
\$\begingroup\$

Octave / MATLAB, 4 bytes

@psi

Try it online!

answered Sep 1, 2020 at 17:25
\$\endgroup\$
4
\$\begingroup\$

R, (削除) 52 (削除ここまで) (削除) 51 (削除ここまで) (削除) 45 (削除ここまで) 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.

answered Sep 1, 2020 at 16:03
\$\endgroup\$
1
  • 1
    \$\begingroup\$ You can even increase the bound to 9e9 without adding bytes. :-) \$\endgroup\$ Commented Sep 1, 2020 at 20:19
2
\$\begingroup\$

Python 3, 38 bytes

from scipy.special import*
f=polygamma

Try it online!

answered Sep 1, 2020 at 22:59
\$\endgroup\$
2
  • 1
    \$\begingroup\$ You can remove 2 bytes by not naming it, even though it does make testing it slightly weird: Try it online! \$\endgroup\$ Commented Sep 1, 2020 at 23:07
  • 2
    \$\begingroup\$ @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! :-) \$\endgroup\$ Commented Sep 1, 2020 at 23:09
2
\$\begingroup\$

Java, (削除) 168 (削除ここまで) (削除) 148 (削除ここまで) 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

I used the same algorithm as in @Arnauld's JavaScript answer. Please vote his answer up.

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

answered Sep 3, 2020 at 10:27
\$\endgroup\$
1
  • 1
    \$\begingroup\$ 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) \$\endgroup\$ Commented Sep 3, 2020 at 13:38
1
\$\begingroup\$

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)
answered Sep 1, 2020 at 15:56
\$\endgroup\$
1
\$\begingroup\$

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!

answered Sep 1, 2020 at 16:03
\$\endgroup\$
1
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ Commented Sep 1, 2020 at 16:18

Your Answer

Draft saved
Draft discarded

Sign up or log in

Sign up using Google
Sign up using Email and Password

Post as a guest

Required, but never shown

Post as a guest

Required, but never shown

By clicking "Post Your Answer", you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.