Implementar a função Polygamma

Question

A função de ordem do Polygamma \$m\$, \$\psi^{(m)}(z)\$, é o \$(m + 1)\$a derivada do logaritmo da função gama , que também é a \$m\$derivada da função digamma . Sua tarefa é pegar um inteiro \$m\$e um número real positivo \$z\$e saída \$\psi^{(m)}(z)\$

Definições

Para aqueles que não estão familiarizados com as funções acima (Gamma, Digamma e Polygamma), aqui estão algumas definições diferentes para cada um:

\$\Gamma(z)\$

• A função gama é uma extensão do fatorial ( \$x! = 1\cdot2\cdot3\cdots(x-1)\cdot(x)\$) para números reais
• \$\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)\$

• A função digamma é a derivada logarítmica da função gama
• \$\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)\$

• A função poligama da ordem \$m\$é o \$m\$derivada da função digamma
• \$\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}}\$

Tarefa

Você deve pegar duas entradas, um número natural \$m\$e um número real positivo \$z\$, e saída \$\psi^{(m)}(z)\$. As entradas e saídas sempre caberão dentro dos limites numéricos de sua linguagem, mas seu algoritmo deve funcionar teoricamente para todas e quaisquer entradas .

Como a saída geralmente será um número real, em vez de um inteiro, a saída deve ser correta para pelo menos 10 algarismos significativos. Zeros à direita podem ser omitidos para valores exatos. Por exemplo, se a saída for um inteiro, os decimais finais 0não são necessários, mas são permitidos se você quiser.

Este é o código-golfe, então o código mais curto em bytes vence.

Casos de teste

Os resultados podem ser diferentes devido a imprecisões de ponto flutuante, a biblioteca scipy do Python foi usada para gerar os valores. Os valores são arredondados para 15d.p., salvo indicação em contrário.

 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

Esta tabela possui os valores de \$\psi^{(m)}(z)\$para \$0 \le m \le 9\$e \$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  |
+---+------------------------+---------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+--------------------+

Answer 1

JavaScript (ES7),  68 66 61  59 bytes

Espera (m)(z).

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

Experimente online!

Isso é baseado na seguinte representação de série ( da Wikipedia ):

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

Comentado

(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

Answer 2

Mathematica, 9 bytes

PolyGamma

Claro que o Mathematica tem um builtin para isso

Answer 3

R , 8 bytes

psigamma

Experimente online!

Recebe entradas z,m(equipamento de teste roubado da resposta de Dominic ).

R tem um embutido como parte de suas Funções Especiais de Matemática, incluindo várias formas da função gama.

Answer 4

Mathematica, 32 bytes (30 caracteres)

Sem qualquer Gammaembutido relacionado, usa a fórmula de Bubbler

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

Experimente online!

Mathematica, 29 bytes

Sem PolyGamma[z]ouPolyGamma[n, z]

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

Experimente online!

Mathematica, 27 bytes

Com PolyGamma[z](isto é o equivalente da função digamma, ou \$\large\psi^0(z)\$)

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

Experimente online!

Alguns programas do Mathematica que não usam o builtin PolyGamma[n, z].

Answer 5

05AB1E , 16 15 bytes

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

-1 byte graças a @ovs .

A primeira entrada é \$z\$, a segunda entrada é \$m\$.

Experimente online ou verifique todos os casos de teste .

Explicação:

Usa o mesmo algoritmo da resposta JavaScript de @Arnauld , portanto, certifique-se de votar nele.
Ou, para ser mais preciso, usa o algoritmo:

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

NOTA: Se houver alguma imprecisão muito pequena nos decimais, o ₄n( \$1\text{,}000\text{,}000\$) pode ser substituído por žm( \$9\text{,}876\text{,}543\text{,}210\$), embora seja muito lento para executar no TIO nesse caso.

Answer 6

Octave / MATLAB, 4 bytes

@psi

Experimente online!

Answer 7

APL (Dyalog Unicode) 18,0, 20 bytes

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

Experimente online!

-2 bytes graças a Adám e ngn.

---

APL (Dyalog Unicode) , 22 bytes

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

Experimente online!

O argumento esquerdo é \$m\$, o argumento certo é \$z\$.

Usa uma ligeira modificação da fórmula usada por outras respostas:

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

Como funciona

{+/(!⍺)÷(-⍵+⍳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

Answer 8

R , 52 51 45 44 bytes

Editar: -1 + 1 bytes graças a Giuseppe, que também apontou que já existe uma função R embutida psigamma, que resolve a tarefa para apenas 8 bytes

Edit2: ... e -6 mais bytes graças a Robin Ryder

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

Experimente online!

Usa a mesma fórmula da resposta de Arnauld.
Representações de série como essa são muito adequadas para R como uma linguagem nativamente vetorizada.

Altere 1e4para valores mais altos (até 9e9sem aumentar a contagem de bytes) para obter uma precisão cada vez maior e um tempo de execução mais lento.

Answer 9

Python 3 , 38 bytes

from scipy.special import*
f=polygamma

Experimente online!

Answer 10

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

Explicação

Usei o mesmo algoritmo da resposta JavaScript de @Arnauld . Por favor, vote na resposta dele.

Por conveniência, aqui a versão formatada automaticamente:

(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;
    };

Um código Java típico: bastante detalhado. Pelo menos minha versão.

Editar: pode economizar 20 bytes graças a @user

Editar: economizou ainda mais bytes graças ao @ceilingcat

Experimente online

Answer 11

Pyth , 24 bytes

**^_1JhhQ*FhQsm^+deQ_JCG

Experimente online! (o link aponta para um código ligeiramente diferente que soma os 1e5termos em vez de 1.56e62termos para tornar o código executável e evitar erros de estouro)

Explicação

Usa o mesmo algoritmo da resposta de JavaScript de @Arnauld , portanto, certifique-se de votar nele.

**^_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)

Answer 12

Scala , 68 bytes

Economizei 3 bytes e consertei minha resposta graças a Arnauld

Usa o algoritmo da resposta de Arnauld

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

Experimente online!