V_GAMMALNS Log of Gamma(x) for positive or negative real x [y,s]=(x) Inputs: x real matrix of input values Outputs: y log(gamma(x)) or, if s is present, log(abs(gamma(x))) s sign(gamma(x)) If s is not specified then y will b complex if gamma(x) is negative (i.e. if floor(x) is negative and odd). Uses gamma(x) = pi/gamma(1-z)/sin(pi*z) which is 5.5.3 from [1] [1] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8.
0001 function [y,s]=v_gammalns(x) 0002 %V_GAMMALNS Log of Gamma(x) for positive or negative real x [y,s]=(x) 0003 % 0004 % Inputs: x real matrix of input values 0005 % 0006 % Outputs: y log(gamma(x)) or, if s is present, log(abs(gamma(x))) 0007 % s sign(gamma(x)) 0008 % 0009 % If s is not specified then y will b complex if gamma(x) is negative (i.e. if floor(x) is negative and odd). 0010 % Uses gamma(x) = pi/gamma(1-z)/sin(pi*z) which is 5.5.3 from [1] 0011 % 0012 % [1] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. 0013 % NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8. 0014 0015 % Copyright (C) Mike Brookes 2017 0016 % Version: $Id: v_gammalns.m 10865 2018-09-21 17:22:45Z dmb $ 0017 % 0018 % VOICEBOX is a MATLAB toolbox for speech processing. 0019 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0020 % 0021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0022 % This program is free software; you can redistribute it and/or modify 0023 % it under the terms of the GNU General Public License as published by 0024 % the Free Software Foundation; either version 2 of the License, or 0025 % (at your option) any later version. 0026 % 0027 % This program is distributed in the hope that it will be useful, 0028 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0030 % GNU General Public License for more details. 0031 % 0032 % You can obtain a copy of the GNU General Public License from 0033 % http://www.gnu.org/copyleft/gpl.html or by writing to 0034 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0036 persistent l 0037 m=x<=0; 0038 y=zeros(size(x)); 0039 s=ones(size(x)); 0040 if ~all(m(:)) % non-negative x values 0041 y(~m)=gammaln(x(~m)); 0042 end 0043 f=m & (x==fix(x)); % find and non-positive integers 0044 if any(f(:)) 0045 y(f)=Inf; % set output to infinity 0046 m=m & ~f; % do not consider these values furhter 0047 end 0048 if any(m(:)) % negative x values 0049 if isempty(l) 0050 l=log(pi); 0051 end 0052 t=sin(pi*x(m)); 0053 if nargout>1 0054 p=t<0; 0055 s(m)=1-2*(p); 0056 y(m)=l-gammaln(1-x(m))-log(abs(t)); 0057 else 0058 y(m)=l-gammaln(1-x(m))-log(t); 0059 end 0060 end 0061