!--- ! Copyright(C) 1996 Takuya OOURA ! 2011 Youhei SASAKI ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: ! 1. Redistributions of source code must retain the above copyright ! notice, this list of conditions and the following disclaimer. ! 2. Redistributions in binary form must reproduce the above copyright ! notice, this list of conditions and the following disclaimer in the ! documentation and/or other materials provided with the distribution. !! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE HOLDERS OR ! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, ! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, ! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR ! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF ! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING ! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! !++ ! !HISTORY 2009/02/06 S.Takehiro: ! downloaded from http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.tar.gz ! 2011/02/18 Y.SASAKI ! rewrite fortran >= 90 ! !+++ != TITLE Gamma function produced by Dr. Ooura, RIMS, Kyoto Univ. ! !* log Gamma function in double precision !* upstream source url: http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.tar.gz ! * rewrite for Fortran 90 by Youhei SASAKI ! real(8) function dlgamma(x) implicit none real(8), intent(in) :: x real(8) :: a(0 : 21) real(8) :: b(0 : 97) real(8) :: c(0 : 64) real(8) :: d(0 : 6) real(8) :: y, t, v, w integer :: i, k real(8), parameter :: pi = 3.141592653589793238d0 data (a(i), i = 0, 10) / & & 0.00009967270908702825d0, -0.00019831672170162227d0, & & -0.00117085315349625822d0, 0.00722012810948319552d0, & & -0.00962213009367802970d0, -0.04219772092994235254d0, & & 0.16653861065243609743d0, -0.04200263501129018037d0, & & -0.65587807152061930091d0, 0.57721566490153514421d0, & & 0.99999999999999999764d0 / data (a(i), i = 11, 21) / & & 0.00004672097259011420d0, -0.00006812300803992063d0, & & -0.00132531159076610073d0, 0.00733521178107202770d0, & & -0.00968095666383935949d0, -0.04217642811873540280d0, & & 0.16653313644244428256d0, -0.04200165481709274859d0, & & -0.65587818792782740945d0, 0.57721567315209190522d0, & & 0.99999999973565236061d0 / data (b(i), i = 0, 13) / & & -0.00000000004587497028d0, 0.00000000019023633960d0, & & 0.00000000086377323367d0, 0.00000001155136788610d0, & & -0.00000002556403058605d0, -0.00000015236723372486d0, & & -0.00000316805106385740d0, 0.00000122903704923381d0, & & 0.00002334372474572637d0, 0.00111544038088797696d0, & & 0.00344717051723468982d0, 0.03198287045148788384d0, & & -0.32705333652955399526d0, 0.40120442440953927615d0 / data (b(i), i = 14, 27) / & & -0.00000000005184290387d0, -0.00000000083355121068d0, & & -0.00000000256167239813d0, 0.00000001455875381397d0, & & 0.00000013512178394703d0, 0.00000029898826810905d0, & & -0.00000358107254612779d0, -0.00002445260816156224d0, & & -0.00004417127762011821d0, 0.00112859455189416567d0, & & 0.00804694454346728197d0, 0.04919775747126691372d0, & & -0.24818372840948854178d0, 0.11071780856646862561d0 / data (b(i), i = 28, 41) / & & 0.00000000030279161576d0, 0.00000000160742167357d0, & & -0.00000000405596009522d0, -0.00000005089259920266d0, & & -0.00000002029496209743d0, 0.00000135130272477793d0, & & 0.00000391430041115376d0, -0.00002871505678061895d0, & & -0.00023052137536922035d0, 0.00045534656385400747d0, & & 0.01153444585593040046d0, 0.07924014651650476036d0, & & -0.12152192626936502982d0, -0.07916438300260539592d0 / data (b(i), i = 42, 55) / & & -0.00000000050919149580d0, -0.00000000115274986907d0, & & 0.00000001237873512188d0, 0.00000002937383549209d0, & & -0.00000030621450667958d0, -0.00000077409414949954d0, & & 0.00000816753874325579d0, 0.00002412433382517375d0, & & -0.00026061217606063700d0, -0.00091000087658659231d0, & & 0.01068093850598380797d0, 0.11395654404408482305d0, & & 0.07209569059984075595d0, -0.10971041451764266684d0 / data (b(i), i = 56, 69) / & & 0.00000000040119897187d0, -0.00000000013224526679d0, & & -0.00000001002723190355d0, 0.00000002569249716518d0, & & 0.00000020336011868466d0, -0.00000118097682726060d0, & & -0.00000300660303810663d0, 0.00004402212897757763d0, & & -0.00001462405876235375d0, -0.00164873795596001280d0, & & 0.00513927520866443706d0, 0.13843580753590579416d0, & & 0.32730190978254056722d0, 0.08588339725978624973d0 / data (b(i), i = 70, 83) / & & -0.00000000015413428348d0, 0.00000000064905779353d0, & & 0.00000000160702811151d0, -0.00000002655645793815d0, & & 0.00000007619544277956d0, 0.00000047604380765353d0, & & -0.00000490748870866195d0, 0.00000821513040821212d0, & & 0.00014804944070262948d0, -0.00122152255762163238d0, & & -0.00087425289205498532d0, 0.14438703699657968310d0, & & 0.61315889733595543766d0, 0.55513708159976477557d0 / data (b(i), i = 84, 97) / & & 0.00000000001049740243d0, -0.00000000025832017855d0, & & 0.00000000139591845075d0, -0.00000000021177278325d0, & & -0.00000005082950464905d0, 0.00000037801785193343d0, & & -0.00000073982266659145d0, -0.00001088918441519888d0, & & 0.00012491810452478905d0, -0.00049171790705139895d0, & & -0.00425707089448266460d0, 0.13595080378472757216d0, & & 0.89518356003149514744d0, 1.31073912535196238583d0 / data (c(i), i = 0, 12) / & & 0.0000000116333640008d0, -0.0000000833156123568d0, & & 0.0000003832869977018d0, -0.0000015814047847688d0, & & 0.0000065010672324100d0, -0.0000274514060128677d0, & & 0.0001209015360925566d0, -0.0005666333178228163d0, & & 0.0029294103665559733d0, -0.0180340086069185819d0, & & 0.1651788780501166204d0, 1.1031566406452431944d0, & & 1.2009736023470742248d0 / data (c(i), i = 13, 25) / & & 0.0000000013842760642d0, -0.0000000069417501176d0, & & 0.0000000342976459827d0, -0.0000001785317236779d0, & & 0.0000009525947257118d0, -0.0000052483007560905d0, & & 0.0000302364659535708d0, -0.0001858396115473822d0, & & 0.0012634378559425382d0, -0.0102594702201954322d0, & & 0.1243625515195050218d0, 1.3888709263595291174d0, & & 2.4537365708424422209d0 / data (c(i), i = 26, 38) / & & 0.0000000001298977078d0, -0.0000000008029574890d0, & & 0.0000000049454846150d0, -0.0000000317563534834d0, & & 0.0000002092136698089d0, -0.0000014252023958462d0, & & 0.0000101652510114008d0, -0.0000774550502862323d0, & & 0.0006537746948291078d0, -0.0066014912535521830d0, & & 0.0996711934948138193d0, 1.6110931485817511402d0, & & 3.9578139676187162939d0 / data (c(i), i = 39, 51) / & & 0.0000000000183995642d0, -0.0000000001353537034d0, & & 0.0000000009984676809d0, -0.0000000076346363974d0, & & 0.0000000599311464148d0, -0.0000004868554120177d0, & & 0.0000041441957716669d0, -0.0000377160856623282d0, & & 0.0003805693126824884d0, -0.0045979851178130194d0, & & 0.0831422678749791178d0, 1.7929113303999329439d0, & & 5.6625620598571415285d0 / data (c(i), i = 52, 64) / & & 0.0000000000034858778d0, -0.0000000000297587783d0, & & 0.0000000002557677575d0, -0.0000000022705728282d0, & & 0.0000000207024992450d0, -0.0000001954426390917d0, & & 0.0000019343161886722d0, -0.0000204790249102570d0, & & 0.0002405181940241215d0, -0.0033842087561074799d0, & & 0.0713079483483518997d0, 1.9467574842460867884d0, & & 7.5343642367587329552d0 / data (d(i), i = 0, 6) / & & -0.00163312359200500807d0, 0.00083644533703385956d0, & & -0.00059518947575728181d0, 0.00079365057505415415d0, & & -0.00277777777735463043d0, 0.08333333333333309869d0, & & 0.91893853320467274178d0 / !! main loop w = x if (x .lt. 0.0d0) w = 1 - x if (w .lt. 0.5d0) then k = 0 if (w .ge. 0.25d0) k = 11 y = ((((((((((a(k) * w + a(k + 1)) * w + & & a(k + 2)) * w + a(k + 3)) * w + a(k + 4)) * w + & & a(k + 5)) * w + a(k + 6)) * w + a(k + 7)) * w + & & a(k + 8)) * w + a(k + 9)) * w + a(k + 10)) * w y = -dlog(y) else if (w .lt. 3.5d0) then t = w - 4.5d0 / (w + 0.5d0) k = idint(t + 4) t = t - (k - 3.5d0) k = k * 14 y = ((((((((((((b(k) * t + b(k + 1)) * t + & & b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & & b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & & b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & & b(k + 11)) * t + b(k + 12)) * t + b(k + 13) else if (w .lt. 8) then k = (idint(w)) - 3 t = w - (k + 3.5d0) k = k * 13 y = (((((((((((c(k) * t + c(k + 1)) * t + & & c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + & & c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + & & c(k + 8)) * t + c(k + 9)) * t + c(k + 10)) * t + & & c(k + 11)) * t + c(k + 12) else v = 1 / w t = v * v y = (((((d(0) * t + d(1)) * t + d(2)) * t + & & d(3)) * t + d(4)) * t + d(5)) * v + d(6) y = y + ((w - 0.5d0) * dlog(w) - w) end if if (x .lt. 0) y = dlog(pi / dsin(pi * x)) - y dlgamma = y end function dlgamma