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 -> ψ⁽ᵐ⁾(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  |
+---+------------------------+---------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+
Improve this question
\$\endgroup\$
3
  • \$\begingroup\$ Sandbox (deleted). Please don't hesitate to comment any suggestions. \$\endgroup\$
    – caird coinheringaahin g ♦
    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\$
    – Carla Cvekla
    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\$
    – caird coinheringaahin g ♦
    Sep 3, 2020 at 17:37

12 Answers 12

Reset to default
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
Improve this answer
\$\endgroup\$
11
\$\begingroup\$

Mathematica, 9 bytes

PolyGamma

Of course Mathematica has a builtin for this

Improve this answer
\$\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\$
    – user
    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].

Improve this answer
\$\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.

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

05AB1E, 16 15 bytes

₄nÝ+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}}$$

₄         # 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.

Improve this answer
\$\endgroup\$
2
  • 1
    \$\begingroup\$ I>mz can be one byte shorter as I±m. \$\endgroup\$
    – ovs
    Sep 1, 2020 at 15:19
  • \$\begingroup\$ @ovs Thanks! I now notice Arnauld had something similar in his answer. \$\endgroup\$
    – Kevin Cruijssen
    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
Improve this answer
\$\endgroup\$
2
  • \$\begingroup\$ shorter as a train: +/!⍤⊣÷(-(⍳1e6)+⊢)*1+⊣ (but crashes on tio) \$\endgroup\$
    – ngn
    Sep 2, 2020 at 7:41
  • \$\begingroup\$ 1e6 can be !9 \$\endgroup\$
    – Adám
    Sep 3, 2020 at 11:28
5
\$\begingroup\$

Octave / MATLAB, 4 bytes

@psi

Try it online!

Improve this answer
\$\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.

Improve this answer
\$\endgroup\$
1
  • 1
    \$\begingroup\$ You can even increase the bound to 9e9 without adding bytes. :-) \$\endgroup\$
    – Robin Ryder
    Sep 1, 2020 at 20:19
2
\$\begingroup\$

Python 3, 38 bytes

from scipy.special import*
f=polygamma

Try it online!

Improve this answer
\$\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\$
    – caird coinheringaahin g ♦
    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\$
    – Noodle9
    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

Improve this answer
\$\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\$
    – user
    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)
Improve this answer
\$\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!

Improve this answer
\$\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\$
    – Arnauld
    Sep 1, 2020 at 16:18

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.