GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
psi.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2016-2017 CarnĂ« Draug
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include "ov.h"
28 #include "defun.h"
29 #include "error.h"
30 #include "dNDArray.h"
31 #include "fNDArray.h"
32 
33 #include "lo-specfun.h"
34 
35 DEFUN (psi, args, ,
36  doc: /* -*- texinfo -*-
37 @deftypefn {} {} psi (@var{z})
38 @deftypefnx {} {} psi (@var{k}, @var{z})
39 Compute the psi (polygamma) function.
40 
41 The polygamma functions are the @var{k}th derivative of the logarithm
42 of the gamma function. If unspecified, @var{k} defaults to zero. A value
43 of zero computes the digamma function, a value of 1, the trigamma function,
44 and so on.
45 
46 The digamma function is defined:
47 
48 @tex
49 $$
50 \Psi (z) = {d (log (\Gamma (z))) \over dx}
51 $$
52 @end tex
53 @ifnottex
54 
55 @example
56 @group
57 psi (z) = d (log (gamma (z))) / dx
58 @end group
59 @end example
60 
61 @end ifnottex
62 
63 When computing the digamma function (when @var{k} equals zero), @var{z}
64 can have any value real or complex value. However, for polygamma functions
65 (@var{k} higher than 0), @var{z} must be real and non-negative.
66 
67 @seealso{gamma, gammainc, gammaln}
68 @end deftypefn */)
69 {
70  int nargin = args.length ();
71 
72  if (nargin < 1 || nargin > 2)
73  print_usage ();
74 
75  const octave_value oct_z = (nargin == 1) ? args(0) : args(1);
76  const octave_idx_type k = (nargin == 1) ? 0 : args(0).idx_type_value ("psi: K must be an integer");
77  if (k < 0)
78  error ("psi: K must be non-negative");
79 
81 
82  if (k == 0)
83  {
84 #define FLOAT_BRANCH(T, A, M, E) \
85  if (oct_z.is_ ## T ##_type ()) \
86  { \
87  const A ## NDArray z = oct_z.M ## array_value (); \
88  A ## NDArray psi_z (z.dims ()); \
89  \
90  const E* zv = z.data (); \
91  E* psi_zv = psi_z.fortran_vec (); \
92  const octave_idx_type n = z.numel (); \
93  for (octave_idx_type i = 0; i < n; i++) \
94  *psi_zv++ = octave::math::psi (*zv++); \
95  \
96  retval = psi_z; \
97  }
98 
99  if (oct_z.is_complex_type ())
100  {
101  FLOAT_BRANCH(double, Complex, complex_, Complex)
102  else FLOAT_BRANCH(single, FloatComplex, float_complex_, FloatComplex)
103  else
104  error ("psi: Z must be a floating point");
105  }
106  else
107  {
108  FLOAT_BRANCH(double, , , double)
109  else FLOAT_BRANCH(single, Float, float_, float)
110  else
111  error ("psi: Z must be a floating point");
112  }
113 
114 #undef FLOAT_BRANCH
115  }
116  else
117  {
118  if (! oct_z.is_real_type ())
119  error ("psi: Z must be real value for polygamma (K > 0)");
120 
121 #define FLOAT_BRANCH(T, A, M, E) \
122  if (oct_z.is_ ## T ##_type ()) \
123  { \
124  const A ## NDArray z = oct_z.M ## array_value (); \
125  A ## NDArray psi_z (z.dims ()); \
126  \
127  const E* zv = z.data (); \
128  E* psi_zv = psi_z.fortran_vec (); \
129  const octave_idx_type n = z.numel (); \
130  for (octave_idx_type i = 0; i < n; i++) \
131  { \
132  if (*zv < 0) \
133  error ("psi: Z must be non-negative for polygamma (K > 0)"); \
134  \
135  *psi_zv++ = octave::math::psi (k, *zv++); \
136  } \
137  retval = psi_z; \
138  }
139 
140  FLOAT_BRANCH(double, , , double)
141  else FLOAT_BRANCH(single, Float, float_, float)
142  else
143  error ("psi: Z must be a floating point for polygamma (K > 0)");
144 
145 #undef FLOAT_BRANCH
146  }
147 
148  return retval;
149 }
150 
151 /*
152 %!shared em
153 %! em = 0.577215664901532860606512090082402431042; # Euler-Mascheroni Constant
154 
155 %!assert (psi (ones (7, 3, 5)), repmat (-em, [7 3 5]))
156 %!assert (psi ([0 1]), [-Inf -em])
157 %!assert (psi ([-20:1]), [repmat(-Inf, [1 21]) -em])
158 %!assert (psi (single ([0 1])), single ([-Inf -em]))
159 
160 ## Abramowitz and Stegun, page 258, eq 6.3.5
161 %!test
162 %! z = [-100:-1 1:200] ./ 10; # drop the 0
163 %! assert (psi (z + 1), psi (z) + 1 ./ z, eps*1000);
164 
165 ## Abramowitz and Stegun, page 258, eq 6.3.2
166 %!assert (psi (1), -em)
167 
168 ## Abramowitz and Stegun, page 258, eq 6.3.3
169 %!assert (psi (1/2), -em - 2 * log (2))
170 
171 ## The following tests are from Pascal Sebah and Xavier Gourdon (2002)
172 ## "Introduction to the Gamma Function"
173 
174 ## Interesting identities of the digamma function, in section of 5.1.3
175 %!assert (psi (1/3), - em - (3/2) * log(3) - ((sqrt (3) / 6) * pi), eps*10)
176 %!assert (psi (1/4), - em -3 * log (2) - pi/2, eps*10)
177 %!assert (psi (1/6), - em -2 * log (2) - (3/2) * log (3) - ((sqrt (3) / 2) * pi), eps*10)
178 
179 ## First 6 zeros of the digamma function, in section of 5.1.5 (and also on
180 ## Abramowitz and Stegun, page 258, eq 6.3.19)
181 %!assert (psi ( 1.46163214496836234126265954232572132846819620400644), 0, eps)
182 %!assert (psi (-0.504083008264455409258269304533302498955385182368579), 0, eps*2)
183 %!assert (psi (-1.573498473162390458778286043690434612655040859116846), 0, eps*2)
184 %!assert (psi (-2.610720868444144650001537715718724207951074010873480), 0, eps*10)
185 %!assert (psi (-3.635293366436901097839181566946017713948423861193530), 0, eps*10)
186 %!assert (psi (-4.653237761743142441714598151148207363719069416133868), 0, eps*100)
187 
188 ## Tests for complex values
189 %!shared z
190 %! z = [-100:-1 1:200] ./ 10; # drop the 0
191 
192 ## Abramowitz and Stegun, page 259 eq 6.3.10
193 %!assert (real (psi (i*z)), real (psi (1 - i*z)))
194 
195 ## Abramowitz and Stegun, page 259 eq 6.3.11
196 %!assert (imag (psi (i*z)), 1/2 .* 1./z + 1/2 * pi * coth (pi * z), eps *10)
197 
198 ## Abramowitz and Stegun, page 259 eq 6.3.12
199 %!assert (imag (psi (1/2 + i*z)), 1/2 * pi * tanh (pi * z), eps*10)
200 
201 ## Abramowitz and Stegun, page 259 eq 6.3.13
202 %!assert (imag (psi (1 + i*z)), - 1./(2*z) + 1/2 * pi * coth (pi * z), eps*10)
203 
204 ## Abramowitz and Stegun, page 260 eq 6.4.5
205 %!test
206 %! for z = 0:20
207 %! assert (psi (1, z + 0.5),
208 %! 0.5 * (pi^2) - 4 * sum ((2*(1:z) -1) .^(-2)),
209 %! eps*10);
210 %! endfor
211 
212 ## Abramowitz and Stegun, page 260 eq 6.4.6
213 %!test
214 %! z = 0.1:0.1:20;
215 %! for n = 0:8
216 %! ## our precision goes down really quick when computing n is too high.
217 %! assert (psi (n, z+1),
218 %! psi (n, z) + ((-1)^n) * factorial (n) * (z.^(-n-1)), 0.1);
219 %! endfor
220 
221 ## Test input validation
222 %!error psi ()
223 %!error psi (1, 2, 3)
224 %!error <Z must be> psi ("non numeric")
225 %!error <conversion of 5.3 to int.* value failed> psi (5.3, 1)
226 %!error <K must be non-negative> psi (-5, 1)
227 %!error <Z must be non-negative for polygamma> psi (5, -1)
228 %!error <Z must be a floating point> psi (5, uint8 (-1))
229 %!error <Z must be real value for polygamma> psi (5, 5i)
230 
231 */
double psi(double z)
Definition: lo-specfun.cc:3714
bool is_real_type(void) const
Definition: ov.h:667
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
for large enough k
Definition: lu.cc:606
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
void error(const char *fmt,...)
Definition: error.cc:570
#define FLOAT_BRANCH(T, A, M, E)
issues an error eealso single
Definition: ov-bool-mat.cc:594
JNIEnv void * args
Definition: ov-java.cc:67
int nargin
Definition: graphics.cc:10115
bool is_complex_type(void) const
Definition: ov.h:670
octave_value retval
Definition: data.cc:6294
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
std::complex< double > Complex
Definition: oct-cmplx.h:31