6月25日(日)にShibuya.lispのもくもく会がありました。14:30にオンラインで集合して、自己紹介と質問タイム（勉強になった！）があり、15時頃～17時頃にもくもくと作業をしていきました。（もくもく会の後の雑談は予定を1時間近くオーバーして盛り上がりました。色々と知らない話を聞くことができて楽しかったです。）
さて2023年2月に、Matthew F. Causley, “The Gamma Function via Interpolation,” Numerical Algorithms, vol.90, pp.687–707 (2022). DOI: 10.1007/s11075-021-01204-8に掲載されているガンマ関数のMatlabコードをCommon Lispに移植しました。
このときは従来法であるスターリング法を用いたgamma_B関数を移植したのですが、新しいアルゴリズムによって精度が向上したうえに高速だというgamma_aaaについては実装できていませんでした。ガンマ関数はベータ分布やt分布の計算に使われる統計分析の基本関数ですので、精度向上も高速化もうれしいことです。今回のもくもく会を利用して、その実装を移植することにしました。
Matlab版では入力として行列が渡されても対応できるような書き方になっていましたが、Common Lisp版では引数として値ひとつだけにしました（必要があればmapcarでリストに適用すればいいですし）。
;;; Gamma Function ;;; ;;; Reference: ;;; Causley, M.F. The gamma function via interpolation. ;;; Numer Algor 90, 687-707 (2022). ;;; https://doi.org/10.1007/s11075-021-01204-8 ;;; preprint: https://arxiv.org/abs/2104.00697 ;;; (defun gamma_aaa (z) "gamma function using a scaled barycentric rational interpolant." (let* ((r 5.51L0) (tt '(#C(0.5d0 0.0L0) #C(0.5d0 -1.0L0) #C(0.5d0 1.0L0) #C(0.5d0 40.0L0) #C(0.5d0 -5.0L0) #C(0.5d0 4.0L0) #C(0.5d0 -3.0L0))) (w '(#C(-0.058033315398988594147056119254557L0 0.0L0) #C(-0.12329392903700113481857414399201L0 -0.05023735799303798155168720995789L0) #C(-0.072017314427899076223482666136988L0 0.029346047538194301729230772934898L0) #C(-0.73570545082472338371815112623153L0 0.35269523425582927078636430451297L0) #C( 0.39424018689617629229715589644911L0 -0.046173606361601587932952384107921L0) #C(-0.10309397777341289259567247427185L0 0.04351009147705412610784847515788L0) #C(-0.17024770255373244953744915619609L0 -0.32884604768510888872512509806256L0))) (f '(#C( 722.24538019924227683077333495021L0 0.0L0) #C( -47.561377245304413463600212708116L0 245.59392283177459148646448738873L0) #C( -47.561377245304413463600212708116L0 -245.59392283177459148646448738873L0) #C( 2.3652595366167963319981026870664L0 -1.1292734670349124925792239082512L0) #C( -7.7668988926260489336073078447953L0 10.095560385519366519702089135535L0) #C( -14.060483019799770332269872596953L0 -14.194015555290931729359726887196L0) #C( -27.239490936407644738892486202531L0 24.743535230939201596811471972615L0))) (j (<= (realpart z) 0.5L0))) (when j (setf z (- 1L0 z))) (let* ((c (mapcar #'(lambda (x) (/ 1 (- z x))) tt)) (g (/ (apply #'+ (mapcar #'(lambda (x y z) (* x y z)) c w f)) (apply #'+ (mapcar #'* c w)))) (gg (* g (exp (- (* (- z 0.5L0) (log (+ z r))) z r))))) (when (realp z) (setf gg (realpart gg))) (when j (multiple-value-bind (m zt) (truncate (realpart z)) (/ (* (expt -1L0 m) pi) (* gg (sin (* pi zt)))))) gg)))
ガンマ関数と階乗には Γ(x) = (x-1)! という関係があるので、整数入力に対しては近似値ではなく整数値が返ってくるほうが望ましい気がします。そこで、以下のように呼び出すようにしました。
(defun gamma (z) (if (and (>= z 0) (integerp z)) (factorial (1- z)) (gamma_aaa (coerce z 'long-float))))
Matlabのgammaと比較してみましょう。
>> gamma(4) ans = 6 >> gamma(5.5) ans = 52.342777784553526 >> gamma(18) ans = 3.556874280960000e+14
CL-USER> (gamma 4) 6 CL-USER> (gamma 5.5) 52.342777784553526d0 CL-USER> (gamma 18) 355687428096000
同じような計算結果になっているので、ヨシ！
