lambdaway :: gamma

5 | list | login |load

| |

_img https://wikimedia.org/api/rest_v1/media/math/render/svg/ad82e51acce72b2ae91dc4b9dfd6bc7544beb74b

_img https://upload.wikimedia.org/wikipedia/commons/thumb/6/6b/Gamma_plus_sin_pi_z.svg/2560px-Gamma_plus_sin_pi_z.svg.png  
{center The gamma function, Γ(z) in blue, plotted along with Γ(z) + sin(πz) in green. Notice the intersection at positive integers. Both are valid analytic continuations of the factorials to the non-integers.}

_h1 gamma function

_p From [[wikipedia|https://en.wikipedia.org/wiki/Gamma_function]] the {b Γ function} is one commonly used extension of the factorial function to complex numbers. The gamma function is defined for all complex numbers except the non-positive integers. For every positive integer n, {b Γ(n)=(n-1)!}.

_p The task is to implement one algorithm to compute the Gamma function (in the real field only). The following lambdatalk code is derived from [[rosettacode.org|https://rosettacode.org/wiki/Gamma_function]], with [[Lanczos|https://en.wikipedia.org/wiki/Lanczos_approximation]] approximation.

{pre
'{def gamma

{def gamma.p
  {A.new     0.99999999999980993
           676.5203681218851
         -1259.1392167224028   
           771.32342877765313
          -176.61502916214059
            12.507343278686905   
            -0.13857109526572012
             9.9843695780195716e-6
             1.5056327351493116e-7
 }}

{def gamma.rec
  {lambda {:x :a :i}
   {if {< :i {A.length {gamma.p}}}
    then {gamma.rec :x 
                    {+ :a {/ {A.get :i {gamma.p}} {+ :x :i}} } 
                    {+ :i 1}}
    else :a
 }}}

{lambda {:x}
  {if {< :x 0.5}
   then {/ {PI}
           {* {sin {* {PI} :x}}
              {gamma {- 1 :x}}}}
   else {let { {:x {- :x 1}}
               {:t {+ {- :x 1} 7 0.5}}
             } {* {sqrt {* 2 {PI}}}
                  {pow :t {+ :x 0.5}}
                  {exp -:t}
                  {gamma.rec :x {A.first {gamma.p}} 1}} 
}}}}
-> {def gamma

'{S.map {lambda {:i}
               {div} Γ(:i) = {gamma :i}} 
       {S.serie -5.5 5.5 0.5}}
-> 
{S.map {lambda {:i}
               {div} Γ(:i) = {gamma :i}} 
       {S.serie -5.5 5.5 0.5}}
}

_h2 Notes
_ul for positive integer values {b Γ(n) = (n-1)!}
_ul for n=0.5 {b Γ(0.5) = √π = {sqrt {PI}}}
_ul ...