Mercurial > octave

annotate libinterp/corefcn/psi.cc @ 26355:bb3ebaa9460e stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
psi.cc: Fix static analyzer detected issues (bug #55347). * psi.cc (Fpsi): Call args(2).xidx_type_value() rather than incorrect version without prepended 'x'.
author Rik <rik@octave.org>
date Tue, 01 Jan 2019 21:12:24 -0800
parents 6652d3823428
children cd44edea6a31
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
1 /*
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
2
25054
6652d3823428 maint: Update copyright dates in all source files.
John W. Eaton <jwe@octave.org>
parents: 24534
diff changeset
3 Copyright (C) 2016-2018 Carnë Draug
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
4
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
5 This file is part of Octave.
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
6
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23582
diff changeset
7 Octave is free software: you can redistribute it and/or modify it
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22407
diff changeset
8 under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23582
diff changeset
9 the Free Software Foundation, either version 3 of the License, or
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22407
diff changeset
10 (at your option) any later version.
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
11
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22407
diff changeset
12 Octave is distributed in the hope that it will be useful, but
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22407
diff changeset
13 WITHOUT ANY WARRANTY; without even the implied warranty of
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22407
diff changeset
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22407
diff changeset
15 GNU General Public License for more details.
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
16
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
17 You should have received a copy of the GNU General Public License
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
18 along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23582
diff changeset
19 <https://www.gnu.org/licenses/>.
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
20
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
21 */
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
22
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21580
diff changeset
23 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21200
diff changeset
24 # include "config.h"
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
25 #endif
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
26
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
27 #include "ov.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
28 #include "defun.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
29 #include "error.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
30 #include "dNDArray.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
31 #include "fNDArray.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
32
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
33 #include "lo-specfun.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
34
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
35 DEFUN (psi, args, ,
21966
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
36 doc: /* -*- texinfo -*-
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
37 @deftypefn {} {} psi (@var{z})
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
38 @deftypefnx {} {} psi (@var{k}, @var{z})
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
39 Compute the psi (polygamma) function.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
40
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
41 The polygamma functions are the @var{k}th derivative of the logarithm
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
42 of the gamma function. If unspecified, @var{k} defaults to zero. A value
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
43 of zero computes the digamma function, a value of 1, the trigamma function,
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
44 and so on.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
45
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
46 The digamma function is defined:
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
47
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
48 @tex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
49 $$
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
50 \Psi (z) = {d (log (\Gamma (z))) \over dx}
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
51 $$
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
52 @end tex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
53 @ifnottex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
54
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
55 @example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
56 @group
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
57 psi (z) = d (log (gamma (z))) / dx
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
58 @end group
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
59 @end example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
60
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
61 @end ifnottex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
62
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
63 When computing the digamma function (when @var{k} equals zero), @var{z}
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
64 can have any value real or complex value. However, for polygamma functions
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
65 (@var{k} higher than 0), @var{z} must be real and non-negative.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
66
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
67 @seealso{gamma, gammainc, gammaln}
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
68 @end deftypefn */)
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
69 {
20939
b17fda023ca6 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20853
diff changeset
70 int nargin = args.length ();
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
71
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
72 if (nargin < 1 || nargin > 2)
20801
a542a9bf177e eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20704
diff changeset
73 print_usage ();
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
74
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
75 const octave_value oct_z = (nargin == 1) ? args(0) : args(1);
26355
bb3ebaa9460e psi.cc: Fix static analyzer detected issues (bug #55347).
Rik <rik@octave.org>
parents: 25054
diff changeset
76 const octave_idx_type k = (nargin == 1) ? 0 : args(0).xidx_type_value ("psi: K must be an integer");
20704
571508c1ed06 eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents: 20424
diff changeset
77 if (k < 0)
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
78 error ("psi: K must be non-negative");
20939
b17fda023ca6 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20853
diff changeset
79
b17fda023ca6 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20853
diff changeset
80 octave_value retval;
b17fda023ca6 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20853
diff changeset
81
b17fda023ca6 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20853
diff changeset
82 if (k == 0)
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
83 {
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
84 #define FLOAT_BRANCH(T, A, M, E) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
85 if (oct_z.is_ ## T ##_type ()) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
86 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
87 const A ## NDArray z = oct_z.M ## array_value (); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
88 A ## NDArray psi_z (z.dims ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
89 \
23449
c763214a8260 maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents: 23220
diff changeset
90 const E *zv = z.data (); \
c763214a8260 maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents: 23220
diff changeset
91 E *psi_zv = psi_z.fortran_vec (); \
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
92 const octave_idx_type n = z.numel (); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
93 for (octave_idx_type i = 0; i < n; i++) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
94 *psi_zv++ = octave::math::psi (*zv++); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
95 \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
96 retval = psi_z; \
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
97 }
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
98
23581
c3075ae020e1 maint: Deprecate is_complex_type and replace with iscomplex.
Rik <rik@octave.org>
parents: 23449
diff changeset
99 if (oct_z.iscomplex ())
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
100 {
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
101 FLOAT_BRANCH(double, Complex, complex_, Complex)
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
102 else FLOAT_BRANCH(single, FloatComplex, float_complex_, FloatComplex)
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
103 else
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
104 error ("psi: Z must be a floating point");
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
105 }
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
106 else
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
107 {
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
108 FLOAT_BRANCH(double, , , double)
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
109 else FLOAT_BRANCH(single, Float, float_, float)
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
110 else
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
111 error ("psi: Z must be a floating point");
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
112 }
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
113
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
114 #undef FLOAT_BRANCH
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
115 }
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
116 else
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
117 {
23582
0cc2011d800e maint: Deprecate is_real_type and replace with isreal.
Rik <rik@octave.org>
parents: 23581
diff changeset
118 if (! oct_z.isreal ())
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
119 error ("psi: Z must be real value for polygamma (K > 0)");
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
120
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
121 #define FLOAT_BRANCH(T, A, M, E) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
122 if (oct_z.is_ ## T ##_type ()) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
123 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
124 const A ## NDArray z = oct_z.M ## array_value (); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
125 A ## NDArray psi_z (z.dims ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
126 \
23449
c763214a8260 maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents: 23220
diff changeset
127 const E *zv = z.data (); \
c763214a8260 maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents: 23220
diff changeset
128 E *psi_zv = psi_z.fortran_vec (); \
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
129 const octave_idx_type n = z.numel (); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
130 for (octave_idx_type i = 0; i < n; i++) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
131 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
132 if (*zv < 0) \
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
133 error ("psi: Z must be non-negative for polygamma (K > 0)"); \
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
134 \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
135 *psi_zv++ = octave::math::psi (k, *zv++); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
136 } \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
137 retval = psi_z; \
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
138 }
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
139
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
140 FLOAT_BRANCH(double, , , double)
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
141 else FLOAT_BRANCH(single, Float, float_, float)
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
142 else
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
143 error ("psi: Z must be a floating point for polygamma (K > 0)");
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
144
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
145 #undef FLOAT_BRANCH
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
146 }
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
147
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
148 return retval;
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
149 }
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
150
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
151 /*
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
152 %!shared em
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
153 %! em = 0.577215664901532860606512090082402431042; # Euler-Mascheroni Constant
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
155 %!assert (psi (ones (7, 3, 5)), repmat (-em, [7 3 5]))
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
156 %!assert (psi ([0 1]), [-Inf -em])
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
157 %!assert (psi ([-20:1]), [repmat(-Inf, [1 21]) -em])
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
158 %!assert (psi (single ([0 1])), single ([-Inf -em]))
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
159
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
160 ## Abramowitz and Stegun, page 258, eq 6.3.5
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
161 %!test
20424
099bdf98f724 psi: Fix failing unit tests on 32-bit systems
Mike Miller <mtmiller@octave.org>
parents: 20161
diff changeset
162 %! z = [-100:-1 1:200] ./ 10; # drop the 0
21580
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
163 %! assert (psi (z + 1), psi (z) + 1 ./ z, eps*1000);
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
164
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
165 ## Abramowitz and Stegun, page 258, eq 6.3.2
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
166 %!assert (psi (1), -em)
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
167
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
168 ## Abramowitz and Stegun, page 258, eq 6.3.3
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
169 %!assert (psi (1/2), -em - 2 * log (2))
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
170
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
171 ## The following tests are from Pascal Sebah and Xavier Gourdon (2002)
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
172 ## "Introduction to the Gamma Function"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
173
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
174 ## Interesting identities of the digamma function, in section of 5.1.3
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
175 %!assert (psi (1/3), - em - (3/2) * log(3) - ((sqrt (3) / 6) * pi), eps*10)
20424
099bdf98f724 psi: Fix failing unit tests on 32-bit systems
Mike Miller <mtmiller@octave.org>
parents: 20161
diff changeset
176 %!assert (psi (1/4), - em -3 * log (2) - pi/2, eps*10)
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
177 %!assert (psi (1/6), - em -2 * log (2) - (3/2) * log (3) - ((sqrt (3) / 2) * pi), eps*10)
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
178
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
179 ## First 6 zeros of the digamma function, in section of 5.1.5 (and also on
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
180 ## Abramowitz and Stegun, page 258, eq 6.3.19)
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
181 %!assert (psi ( 1.46163214496836234126265954232572132846819620400644), 0, eps)
21762
289ee97765ad psi: Adjust tolerances for failing tests on i386 systems (bug #47965)
Mike Miller <mtmiller@octave.org>
parents: 21724
diff changeset
182 %!assert (psi (-0.504083008264455409258269304533302498955385182368579), 0, eps*2)
289ee97765ad psi: Adjust tolerances for failing tests on i386 systems (bug #47965)
Mike Miller <mtmiller@octave.org>
parents: 21724
diff changeset
183 %!assert (psi (-1.573498473162390458778286043690434612655040859116846), 0, eps*2)
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
184 %!assert (psi (-2.610720868444144650001537715718724207951074010873480), 0, eps*10)
20156
bd565f3e0ecb psi: refactor to reduce code duplication.
Carnë Draug <carandraug@octave.org>
parents: 20155
diff changeset
185 %!assert (psi (-3.635293366436901097839181566946017713948423861193530), 0, eps*10)
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
186 %!assert (psi (-4.653237761743142441714598151148207363719069416133868), 0, eps*100)
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
187
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
188 ## Tests for complex values
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
189 %!shared z
20424
099bdf98f724 psi: Fix failing unit tests on 32-bit systems
Mike Miller <mtmiller@octave.org>
parents: 20161
diff changeset
190 %! z = [-100:-1 1:200] ./ 10; # drop the 0
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
191
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
192 ## Abramowitz and Stegun, page 259 eq 6.3.10
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
193 %!assert (real (psi (i*z)), real (psi (1 - i*z)))
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
194
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
195 ## Abramowitz and Stegun, page 259 eq 6.3.11
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
196 %!assert (imag (psi (i*z)), 1/2 .* 1./z + 1/2 * pi * coth (pi * z), eps *10)
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
197
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
198 ## Abramowitz and Stegun, page 259 eq 6.3.12
20424
099bdf98f724 psi: Fix failing unit tests on 32-bit systems
Mike Miller <mtmiller@octave.org>
parents: 20161
diff changeset
199 %!assert (imag (psi (1/2 + i*z)), 1/2 * pi * tanh (pi * z), eps*10)
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
200
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
201 ## Abramowitz and Stegun, page 259 eq 6.3.13
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
202 %!assert (imag (psi (1 + i*z)), - 1./(2*z) + 1/2 * pi * coth (pi * z), eps*10)
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
203
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
204 ## Abramowitz and Stegun, page 260 eq 6.4.5
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
205 %!test
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
206 %! for z = 0:20
21580
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
207 %! assert (psi (1, z + 0.5),
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
208 %! 0.5 * (pi^2) - 4 * sum ((2*(1:z) -1) .^(-2)),
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
209 %! eps*10);
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
210 %! endfor
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
211
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
212 ## Abramowitz and Stegun, page 260 eq 6.4.6
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
213 %!test
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
214 %! z = 0.1:0.1:20;
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
215 %! for n = 0:8
21580
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
216 %! ## our precision goes down really quick when computing n is too high.
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
217 %! assert (psi (n, z+1),
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
218 %! psi (n, z) + ((-1)^n) * factorial (n) * (z.^(-n-1)), 0.1);
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
219 %! endfor
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
220
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
221 ## Test input validation
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
222 %!error psi ()
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
223 %!error psi (1, 2, 3)
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
224 %!error <Z must be> psi ("non numeric")
21830
cbd8cf0a8a5c Fix psi input validation with --enable-64
Mike Miller <mtmiller@octave.org>
parents: 21786
diff changeset
225 %!error <conversion of 5.3 to int.* value failed> psi (5.3, 1)
20704
571508c1ed06 eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents: 20424
diff changeset
226 %!error <K must be non-negative> psi (-5, 1)
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
227 %!error <Z must be non-negative for polygamma> psi (5, -1)
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
228 %!error <Z must be a floating point> psi (5, uint8 (-1))
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
229 %!error <Z must be real value for polygamma> psi (5, 5i)
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
230
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
231 */