Mercurial > octave-antonio

comparison libinterp/corefcn/psi.cc @ 20161:65e22ba879f0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
psi: add support to compute the polygamma function (kth-derivative). * libinterp/corefcn/psi.cc: previously, only the digamma function, k == 0, was being computed. Add support for polygamma function, add tests, and improve documentation. * liboctave/cruft/slatec-fn/dpsifn.f, liboctave/cruft/slatec-fn/psifn.f: the two functions that actually compute the the polygamma functions, copied verbatim from SLATEC, and under public domain. * liboctave/cruft/slatec-fn/module.mk: add dpsifn.f and psifn.f to the build system. * liboctave/numeric/lo-specfun.cc: add new signature for function psi to compute polygamma function that wraps the Fortran DPSIFN and PSIFN functions. * liboctave/numeric/lo-specfun.h: declare new function and document all psi() with doxygen.
author Carnë Draug <carandraug@octave.org>
date Sun, 03 May 2015 22:52:07 +0100
parents bd565f3e0ecb
children
comparison
equal deleted inserted replaced
20157:e410d62ae2c8 20161:65e22ba879f0
32 32
33 #include "lo-specfun.h" 33 #include "lo-specfun.h"
34 34
35 DEFUN (psi, args, , 35 DEFUN (psi, args, ,
36 "-*- texinfo -*-\n\ 36 "-*- texinfo -*-\n\
37 @deftypefn {Function File} {} psi (@var{z})\n\ 37 @deftypefn {Function File} {} psi (@var{z})\n\
38 Compute the psi (digamma) function.\n\ 38 @deftypefnx {Function File} {} psi (@var{k}, @var{z})\n\
39 Compute the psi (polygamma) function.\n\
40 \n\
41 The polygamma functions are the @var{k}th derivative of the logarithm\n\
42 of the gamma function. If unspecified, @var{k} defaults to zero. A value\n\
43 of zero computes the digamma function, a value of 1, the trigamma function,\n\
44 and so on.\n\
45 \n\
46 The digamma function is defined:\n\
39 \n\ 47 \n\
40 @tex\n\ 48 @tex\n\
41 $$\n\ 49 $$\n\
42 \\Psi (z) = {d (log (\\Gamma (z))) \\over dx}\n\ 50 \\Psi (z) = {d (log (\\Gamma (z))) \\over dx}\n\
43 $$\n\ 51 $$\n\
48 psi (z) = d (log (gamma (z))) / dx\n\ 56 psi (z) = d (log (gamma (z))) / dx\n\
49 @end group\n\ 57 @end group\n\
50 @end example\n\ 58 @end example\n\
51 @end ifnottex\n\ 59 @end ifnottex\n\
52 \n\ 60 \n\
61 When computing the digamma function (when @var{k} equals zero), @var{z}\n\
62 can have any value real or complex value. However, for polygamma functions\n\
63 (@var{k} higher than 0), @var{z} must be real and non-negative.\n\
64 \n\
53 @seealso{gamma, gammainc, gammaln}\n\ 65 @seealso{gamma, gammainc, gammaln}\n\
54 @end deftypefn") 66 @end deftypefn")
55 { 67 {
56 octave_value retval; 68 octave_value retval;
57 69
58 const octave_idx_type nargin = args.length (); 70 const octave_idx_type nargin = args.length ();
59 71 if (nargin < 1 || nargin > 2)
60 if (nargin != 1)
61 { 72 {
62 print_usage (); 73 print_usage ();
63 return retval; 74 return retval;
64 } 75 }
65 76
77 const octave_value oct_z = (nargin == 1) ? args(0) : args(1);
78 const octave_idx_type k = (nargin == 1) ? 0 : args(0).idx_type_value ();
79 if (error_state || k < 0)
80 {
81 error ("psi: K must be a non-negative integer");
82 return retval;
83 }
84 else if (k == 0)
85 {
66 #define FLOAT_BRANCH(T, A, M, E) \ 86 #define FLOAT_BRANCH(T, A, M, E) \
67 if (args(0).is_ ## T ##_type ()) \ 87 if (oct_z.is_ ## T ##_type ()) \
68 { \ 88 { \
69 const A ## NDArray z = args(0).M ## array_value (); \ 89 const A ## NDArray z = oct_z.M ## array_value (); \
70 A ## NDArray psi_z (z.dims ()); \ 90 A ## NDArray psi_z (z.dims ()); \
71 \ 91 \
72 const E* zv = z.data (); \ 92 const E* zv = z.data (); \
73 E* psi_zv = psi_z.fortran_vec (); \ 93 E* psi_zv = psi_z.fortran_vec (); \
74 const octave_idx_type n = z.numel (); \ 94 const octave_idx_type n = z.numel (); \
75 for (octave_idx_type i = 0; i < n; i++) \ 95 for (octave_idx_type i = 0; i < n; i++) \
76 psi_zv[i] = psi (zv[i]); \ 96 *psi_zv++ = psi (*zv++); \
77 \ 97 \
78 retval = psi_z; \ 98 retval = psi_z; \
79 } 99 }
80 100
81 if (args(0).is_complex_type ()) 101 if (oct_z.is_complex_type ())
82 { 102 {
83 FLOAT_BRANCH(double, Complex, complex_, Complex) 103 FLOAT_BRANCH(double, Complex, complex_, Complex)
84 else FLOAT_BRANCH(single, FloatComplex, float_complex_, FloatComplex) 104 else FLOAT_BRANCH(single, FloatComplex, float_complex_, FloatComplex)
105 else
106 {
107 error ("psi: Z must be a floating point");
108 }
109 }
85 else 110 else
86 { 111 {
87 error ("psi: Z must be a floating point"); 112 FLOAT_BRANCH(double, , , double)
88 } 113 else FLOAT_BRANCH(single, Float, float_, float)
114 else
115 {
116 error ("psi: Z must be a floating point");
117 }
118 }
119
120 #undef FLOAT_BRANCH
89 } 121 }
90 else 122 else
91 { 123 {
124 if (! oct_z.is_real_type ())
125 {
126 error ("psi: Z must be real value for polygamma (K > 0)");
127 return retval;
128 }
129
130 #define FLOAT_BRANCH(T, A, M, E) \
131 if (oct_z.is_ ## T ##_type ()) \
132 { \
133 const A ## NDArray z = oct_z.M ## array_value (); \
134 A ## NDArray psi_z (z.dims ()); \
135 \
136 const E* zv = z.data (); \
137 E* psi_zv = psi_z.fortran_vec (); \
138 const octave_idx_type n = z.numel (); \
139 for (octave_idx_type i = 0; i < n; i++) \
140 { \
141 if (*zv < 0) \
142 { \
143 error ("psi: Z must be non-negative for polygamma (K > 0)"); \
144 return retval; \
145 } \
146 *psi_zv++ = psi (k, *zv++); \
147 } \
148 retval = psi_z; \
149 }
150
92 FLOAT_BRANCH(double, , , double) 151 FLOAT_BRANCH(double, , , double)
93 else FLOAT_BRANCH(single, Float, float_, float) 152 else FLOAT_BRANCH(single, Float, float_, float)
94 else 153 else
95 { 154 {
96 error ("psi: Z must be a floating point"); 155 error ("psi: Z must be a floating point for polygamma (K > 0)");
97 } 156 }
98 }
99 157
100 #undef FLOAT_BRANCH 158 #undef FLOAT_BRANCH
159 }
101 160
102 return retval; 161 return retval;
103 } 162 }
104 163
105 /* 164 /*
153 %!assert (imag (psi (1/2 + i*z)), 1/2 * pi * tanh (pi * z), eps) 212 %!assert (imag (psi (1/2 + i*z)), 1/2 * pi * tanh (pi * z), eps)
154 213
155 ## Abramowitz and Stegun, page 259 eq 6.3.13 214 ## Abramowitz and Stegun, page 259 eq 6.3.13
156 %!assert (imag (psi (1 + i*z)), - 1./(2*z) + 1/2 * pi * coth (pi * z), eps*10) 215 %!assert (imag (psi (1 + i*z)), - 1./(2*z) + 1/2 * pi * coth (pi * z), eps*10)
157 216
217 ## Abramowitz and Stegun, page 260 eq 6.4.5
218 %!test
219 %! for z = 0:20
220 %! assert (psi (1, z + 0.5), 0.5 * (pi^2) - 4 * sum ((2*(1:z) -1) .^(-2)), eps*10)
221 %! endfor
222
223 ## Abramowitz and Stegun, page 260 eq 6.4.6
224 %!test
225 %! z = 0.1:0.1:20;
226 %! for n = 0:8
227 %! ## our precision goes down really quick when computing n is too high,
228 %! assert (psi (n, z+1), psi (n, z) + ((-1)^n) * factorial (n) * (z.^(-n-1)), 0.1)
229 %! endfor
230
231 ## Test input validation
232 %!error psi ()
233 %!error psi (1, 2, 3)
234 %!error <Z must be> psi ("non numeric")
235 %!error <K must be a non-negative integer> psi (-5, 1)
236 %!error <Z must be non-negative for polygamma> psi (5, -1)
237 %!error <Z must be a floating point> psi (5, uint8 (-1))
238 %!error <Z must be real value for polygamma> psi (5, 5i)
239
158 */ 240 */