Mercurial > octave

annotate libinterp/corefcn/psi.cc @ 29359:7854d5752dd2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: merge stable to default.
author John W. Eaton <jwe@octave.org>
date Wed, 10 Feb 2021 10:10:40 -0500
parents c40a367a84c0 0a5b15007766
children 32c3a5805893
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
1 ////////////////////////////////////////////////////////////////////////
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
2 //
29358
0a5b15007766 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
3 // Copyright (C) 2016-2021 The Octave Project Developers
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
4 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
5 // See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
6 // distribution or <https://octave.org/copyright/>.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
7 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
8 // This file is part of Octave.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
9 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
10 // Octave is free software: you can redistribute it and/or modify it
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
11 // under the terms of the GNU General Public License as published by
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
12 // the Free Software Foundation, either version 3 of the License, or
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
13 // (at your option) any later version.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
14 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
15 // Octave is distributed in the hope that it will be useful, but
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
18 // GNU General Public License for more details.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
19 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
20 // You should have received a copy of the GNU General Public License
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
21 // along with Octave; see the file COPYING. If not, see
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
22 // <https://www.gnu.org/licenses/>.
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
23 //
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
24 ////////////////////////////////////////////////////////////////////////
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
25
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21580
diff changeset
26 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21200
diff changeset
27 # include "config.h"
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
28 #endif
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
29
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
30 #include "ov.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
31 #include "defun.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
32 #include "error.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
33 #include "dNDArray.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
34 #include "fNDArray.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
35
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
36 #include "lo-specfun.h"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
37
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
38 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
39 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
40 @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
41 @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
42 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
43
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
44 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
45 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
46 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
47 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
48
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
49 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
50
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
51 @tex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
52 $$
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
53 \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
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 @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
56 @ifnottex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
57
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
58 @example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
59 @group
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
60 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
61 @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
62 @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
63
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
64 @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
65
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
66 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
67 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
68 (@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
69
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21830
diff changeset
70 @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
71 @end deftypefn */)
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
72 {
20939
b17fda023ca6 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20853
diff changeset
73 int nargin = args.length ();
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 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
76 print_usage ();
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
77
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
78 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
79 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
80 if (k < 0)
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
81 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
82
b17fda023ca6 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20853
diff changeset
83 octave_value retval;
b17fda023ca6 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20853
diff changeset
84
b17fda023ca6 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20853
diff changeset
85 if (k == 0)
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
86 {
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
87 #define FLOAT_BRANCH(T, A, M, E) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
88 if (oct_z.is_ ## T ##_type ()) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
89 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
90 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
91 A ## NDArray psi_z (z.dims ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
92 \
23449
c763214a8260 maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents: 23220
diff changeset
93 const E *zv = z.data (); \
c763214a8260 maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents: 23220
diff changeset
94 E *psi_zv = psi_z.fortran_vec (); \
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
95 const octave_idx_type n = z.numel (); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
96 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
97 *psi_zv++ = octave::math::psi (*zv++); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
98 \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
99 retval = psi_z; \
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
100 }
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
101
23581
c3075ae020e1 maint: Deprecate is_complex_type and replace with iscomplex.
Rik <rik@octave.org>
parents: 23449
diff changeset
102 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
103 {
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
104 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
105 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
106 else
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
107 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
108 }
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
109 else
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
110 {
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
111 FLOAT_BRANCH(double, , , double)
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
112 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
113 else
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
114 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
115 }
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
116
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
117 #undef FLOAT_BRANCH
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
118 }
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
119 else
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
120 {
23582
0cc2011d800e maint: Deprecate is_real_type and replace with isreal.
Rik <rik@octave.org>
parents: 23581
diff changeset
121 if (! oct_z.isreal ())
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
122 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
123
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
124 #define FLOAT_BRANCH(T, A, M, E) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
125 if (oct_z.is_ ## T ##_type ()) \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
126 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
127 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
128 A ## NDArray psi_z (z.dims ()); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
129 \
23449
c763214a8260 maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents: 23220
diff changeset
130 const E *zv = z.data (); \
c763214a8260 maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents: 23220
diff changeset
131 E *psi_zv = psi_z.fortran_vec (); \
22197
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
132 const octave_idx_type n = z.numel (); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
133 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
134 { \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
135 if (*zv < 0) \
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
136 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
137 \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
138 *psi_zv++ = octave::math::psi (k, *zv++); \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
139 } \
e43d83253e28 refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
140 retval = psi_z; \
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
141 }
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
142
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
143 FLOAT_BRANCH(double, , , double)
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
144 else FLOAT_BRANCH(single, Float, float_, float)
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
145 else
20831
35241c4b696c eliminate return statements after calls to error
John W. Eaton <jwe@octave.org>
parents: 20801
diff changeset
146 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
147
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
148 #undef FLOAT_BRANCH
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
149 }
20154
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 return retval;
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
152 }
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
153
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 %!shared em
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
156 %! em = 0.577215664901532860606512090082402431042; # Euler-Mascheroni Constant
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
157
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
158 %!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
159 %!assert (psi ([0 1]), [-Inf -em])
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
160 %!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
161 %!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
162
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
163 ## 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
164 %!test
20424
099bdf98f724 psi: Fix failing unit tests on 32-bit systems
Mike Miller <mtmiller@octave.org>
parents: 20161
diff changeset
165 %! 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
166 %! 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
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.2
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
169 %!assert (psi (1), -em)
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 ## 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
172 %!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
173
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
174 ## 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
175 ## "Introduction to the Gamma Function"
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
176
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
177 ## Interesting identities of the digamma function, in section of 5.1.3
28915
c40a367a84c0 maint: Use Octave convention of space after function name in libinterp/.
Rik <rik@octave.org>
parents: 28888
diff changeset
178 %!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
179 %!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
180 %!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
181
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
182 ## 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
183 ## 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
184 %!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
185 %!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
186 %!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
187 %!assert (psi (-2.610720868444144650001537715718724207951074010873480), 0, eps*10)
20156
bd565f3e0ecb psi: refactor to reduce code duplication.
Carnë Draug <carandraug@octave.org>
parents: 20155
diff changeset
188 %!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
189 %!assert (psi (-4.653237761743142441714598151148207363719069416133868), 0, eps*100)
20155
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
190
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
191 ## Tests for complex values
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
192 %!shared z
20424
099bdf98f724 psi: Fix failing unit tests on 32-bit systems
Mike Miller <mtmiller@octave.org>
parents: 20161
diff changeset
193 %! 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
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.10
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
196 %!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
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.11
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
199 %!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
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.12
20424
099bdf98f724 psi: Fix failing unit tests on 32-bit systems
Mike Miller <mtmiller@octave.org>
parents: 20161
diff changeset
202 %!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
203
1fae49e34a1a psi: add support for complex numbers.
Carnë Draug <carandraug@octave.org>
parents: 20154
diff changeset
204 ## 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
205 %!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
206
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
207 ## 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
208 %!test
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
209 %! for z = 0:20
21580
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
210 %! assert (psi (1, z + 0.5),
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
211 %! 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
212 %! eps*10);
20161
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
213 %! endfor
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
214
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
215 ## 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
216 %!test
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
217 %! 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
218 %! for n = 0:8
21580
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
219 %! ## 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
220 %! assert (psi (n, z+1),
ecce63c99c3f maint: Add semicolons to terminate code in %! blocks.
Rik <rik@octave.org>
parents: 21316
diff changeset
221 %! 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
222 %! endfor
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
223
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
224 ## Test input validation
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
225 %!error psi ()
65e22ba879f0 psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents: 20156
diff changeset
226 %!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
227 %!error <Z must be> psi ("non numeric")
26359
cd44edea6a31 psi.cc: Fix BIST test.
Rik <rik@octave.org>
parents: 26355
diff changeset
228 %!error <K must be an integer> psi ({5.3}, 1)
20704
571508c1ed06 eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents: 20424
diff changeset
229 %!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
230 %!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
231 %!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
232 %!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
233
20154
45565ecec019 New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff changeset
234 */