Do you want it to work for complex or only for real argument?
I've been looking for a numeric Gamma method that gives 34 digits of precision suitable for quad-precision floating-point [2] although I would settle for 31 digits suitable for "double-double".
By the way, on the theoretical front, one function that I did not study, but probably deserves study, is g(z) := Gamma(z+1/2)^(1/z). does this have a convergent continued fraction expansion? It has a divergent asymptotic series .36787944117144232159552377016146086744581113103177*z+ .33805859406623990237027945096151887426851375834020+ .13999922423542209328692750196611434611925705150275/z+ .33493044087723619416023793118971174508783220448207e-1/z^2+ .56720131162423401908596776888208182906669281170266e-2/z^3+ .11415650363382196952422045360200835828013297526599e-2/z^4+ .40066835176315930489725300869291932605805052785655e-4/z^5- .18459730484462595896459837456234342865776092168523e-3/z^6+ .12180790117916786446998042861431093723893158454545e-3/z^7+ .17568169709328468660924962236969118264708016900776e-3/z^8- .23036114978935780189878504071470392725192413835113e-3/z^9- .26501793497134862116045528462277497265836478867190e-3/z^10+ .59068400115606570603275441175404496195447459946218e-3/z^11+ .62050016149812141546059352772177092271212715711371e-3/z^12- .20943016044158043700802700877762421178065590234503e-2/z^13- .21036447470556705667818040016973832680344857168705e-2/z^14+ .99842712425733569346391610042789127327983738573824e-2/z^15+ .97772628167677602737106584758217900245031666550502e-2/z^16- .61986382824479867397792991698201564463192370984501e-1/z^17- .59747962332351720667423814714163223752965683269949e-1/z^18+... but perhaps it becomes convergent when converted to a continued fraction (aka "Pade summation"). This would probably be hard to prove, if true.